RandomVariate returns values outside the support of a PDF

Question

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:

Max[data]
0.15312

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

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}]}], 
    10^6]]]

{-0.135334, 0.135335}

The PDF is in good agreement with the histogram:

Show[
  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]]

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

FullSimplify[Quantile[ProbabilityDistribution[1/(Abs[x]*Log[Abs[x]]^2), 
       {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}