1
$\begingroup$

The following seems buggy to me, but perhaps I'm confused about something. Or maybe there's a workaround that makes it moot?

First, I ask Mathematica to find a distribution that describes a dataset that's roughly uniform between 3 and 4. Then I truncate the result to never be outside the interval [1,7].

d = TruncatedDistribution[{1, 7}, FindDistribution[{3, 3.4, 3.6, 4}]];

The result of that is:

TruncatedDistribution[{3., 7}, ParetoDistribution[7., 2.60877, 8.13978, 3.]]

But plotting the PDF or CDF of that shows nothing and if I ask for something like CDF[d, 5] I get a complex number. How does that happen? Can I stipulate that I only want distributions over real numbers?

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ It looks like you have not enough data points to safely define the distribution. $\endgroup$
    – Daniel Huber
    Aug 9, 2022 at 7:44

3 Answers 3

Reset to default
3
$\begingroup$

I agree with OP that this is some kind of bug, in Version 12.3 at least. Compare:

  • With exact parameters:
CDF[TruncatedDistribution[{1,7},ParetoDistribution[7,2,8,3]],4]
(* (1-(1+7^(-1/8))^(-2))/(1-(1+2^(1/4)/7^(1/8))^(-2))
   the //N of which is 0.9366484688088881 *)
  • With non-exact parameters
CDF[TruncatedDistribution[{1.,7.},ParetoDistribution[7.,2.,8.,3.]],4.]
(* 0.9351246575722288 - 0.000639203984437722*I *)

Clearly this should not give an imaginary part, let alone such a not-so-small imaginary part.

Possible cause. By Trace-ing the above calculation with non-exact parameters, one finds that one expression generated early in the course of the evaluation is

DistributionDomain[ParetoDistribution[7.,2.,8.,3.]]//FullForm
(* Interval[List[2.9999999999999996`,DirectedInfinity[1]]] *)

For some reason, the last argument 3. was transformed to a number very slightly below 3.. It is conceivable that somewhere along the line, this leads to an evaluation similar to

Quiet[NIntegrate[(x-3.)^(-7/8),{x,2.9999999999999996`,4}]]

which generates a not-so-small imaginary part.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Oh ho, thanks for the debugging here! $\endgroup$
    – dreeves
    Aug 9, 2022 at 21:31
2
$\begingroup$

I think decimals cause it. The following works well.

d = TruncatedDistribution[{1, 7}, Rationalize[FindDistribution[{3, 3.4, 3.6, 4}], 0]]
CDF[d, 5] // N

0.975879

Plot[CDF[d, x], {x, 1, 8}]

enter image description here

Addition. Another possible approach consists in the following.

d = TruncatedDistribution[{1, 7}, FindDistribution[{3, 3.4, 3.6, 4}, 
   TargetFunctions -> {UniformDistribution}]]

UniformDistribution[{3.07532,3.97558}]

Improve this answer
$\endgroup$
2
$\begingroup$

If we ignore that feeding just 3 data points to a black box function that estimates a distribution is rarely recommended, there is a simple workaround that avoids obtaining complex numbers.

Look at the "found" distribution:

d0 = FindDistribution[{3, 3.4, 3.6, 4}]
(* ParetoDistribution[7., 2.60877, 8.13978, 3.] *)

The pdf of that distribution is

PDF[d0, x] // FullSimplify

pdf of found distribution

One sees that the pdf is 0 for $x<3$. But when you request for that distribution be truncated, you truncate between 1 and 7. There's no point in truncating below 3 as the pdf is zero below 3. So changing the truncation to 3 and 7 gives the desired result.

d = TruncatedDistribution[{3, 7}, FindDistribution[{3, 3.4, 3.6, 4}]];
PDF[d, x]

pdf of truncated distribution

Plot[PDF[d, x], {x, 2, 8}, PlotRange -> {Automatic, {0, 1}}]

Plot of pdf of truncated distribution

Another alternative is @user64494 's solution: Rationalize the distribution prior to evaluating the pdf or cdf.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Ah, thank you! I would think that truncating a distribution below where the pdf is zero should be a no-op. Sounds like another general workaround for my problem would be to take the interval [a,b] you want to truncate to and then also find the support of the distribution, i.e., where it's non-zero. Call the support [sa, sb]. So then just truncate the distribution to be between max(a, sa) and min(b, sb). You've shown that that works in my example. I'm not sure how to do it in general though. So maybe the Rationalize trick is the way to go. $\endgroup$
    – dreeves
    Aug 9, 2022 at 18:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.