Bug introduced in 8.0.1 and fixed in 10.2.0

Let $X$ be a random variable with pdf:

dist = ProbabilityDistribution[1/(Abs[x]*Log[Abs[x]]^2), {x, -E^-2, E^-2}]

Here are some pseudo-random drawings from it:

data = RandomVariate[dist, 500000];

The domain of support for this random variable is:

{x, -E^-2, E^-2} // N
{x, -0.135335, 0.135335}

and all the generated data should lie within this domain of support. But, if I try something like:


Mathematica 9 returns a number which lies outside the domain of support. So, the Mathematica random number generator is failing — it is generating values that lie outside the domain of support: i.e. not from this pdf.

I reported this under version 8 to WRI, and again under version 9 pre-testing and since it is still here, I am beginning to wonder whether it is just me? Do others get this problem as well? I am running a Mac Pro with OS X 10.6.8

  • 1
    $\begingroup$ I get the same but the following works dist2 = ProbabilityDistribution[ 2/(Abs[x]*Log[Abs[x]]^2), {x, 0, Exp[-2]}]. $\endgroup$ Dec 23, 2012 at 17:59
  • $\begingroup$ Interestingly, it only overshoots at the upper boundary of the range but never at the lower one (v 8.0.4, Win7). $\endgroup$ Dec 23, 2012 at 18:46
  • 1
    $\begingroup$ Colin when you reported this for V8 did they say it was a bug (they never use the word "bug" but you know what I mean). If so worth tagging the question as "bug". $\endgroup$ Dec 23, 2012 at 21:51
  • $\begingroup$ Thanks for the suggestion Mike. I have added the 'bugs' tag. $\endgroup$
    – wolfies
    Dec 24, 2012 at 19:34

1 Answer 1


The integrable singularity of the PDF at the origin is not gracefully dealt with by the underlying solvers.

Presently one can work around the issue by exploiting the symmetry of the PDF:

Through[{Min, Max}[
  sample = RandomVariate[
    TransformedDistribution[(-1)^x y,
      {x \[Distributed] BernoulliDistribution[1/2], 
       y \[Distributed] ProbabilityDistribution[2/(x Log[ x ]^2),
       {x, 0, E^-2}]}], 

{-0.135334, 0.135335}

The PDF is in good agreement with the histogram:

  Histogram[sample, Automatic, "PDF"],
  Plot[1/(Abs[x] Log[ Abs[x] ]^2), {x, -1/E^2, 1/E^2}, 
    PlotStyle -> Directive[Thin, Orange], PlotRange -> 25]]

enter image description here

Alternatively, one could use the closed form expression for Quantile and apply the direct inversion method:

       {x, -E^(-2), E^(-2)}], q], 0 <= q <= 1]

Piecewise[{{-E^(2/(-1 + 2*q)), 2*q < 1}, { E^(2/( 1 - 2*q)), 2*q > 1}}, 0]

Through[{Min, Max}[
  sample2 = Function[q, Sign[2 q - 1] E^(-Abs[2/(1 - 2*q)])][RandomReal[1, 10^6]]]]

{-0.135335, 0.135335}

  • $\begingroup$ The above seems like a complicated way of stating what b.gatessucks posted above (Comment the First). Is there a reason this bug cannot be fixed? And how pervasive are the problems with the Wolfram random number generator? $\endgroup$
    – wolfies
    Dec 24, 2012 at 19:19
  • $\begingroup$ @ColinRose Yes, the above workaround is based on the symmetry and uses the observation noted in the first comment. The issue is caused by inaccurate result of numerical integration of the pdf. There is no reason this particular issue can not be fixed. It is always a matter of resource allocation. The above workaround (especially the direct inversion) is a sound way to get correct samples. $\endgroup$
    – Sasha
    Dec 24, 2012 at 20: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.