# RandomVariate returns values outside the support of a PDF

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

-
I get the same but the following works dist2 = ProbabilityDistribution[ 2/(Abs[x]*Log[Abs[x]]^2), {x, 0, Exp[-2]}]. –  b.gatessucks Dec 23 '12 at 17:59
Interestingly, it only overshoots at the upper boundary of the range but never at the lower one (v 8.0.4, Win7). –  István Zachar Dec 23 '12 at 18:46
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". –  Mike Honeychurch Dec 23 '12 at 21:51
Thanks for the suggestion Mike. I have added the 'bugs' tag. –  wolfies Dec 24 '12 at 19:34

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}

-
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? –  wolfies Dec 24 '12 at 19:19
@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. –  Sasha Dec 24 '12 at 20:18