# Optimize RandomVariate for a $1/x^4$ distribution

I'm working with the following code, where I define a probability distribution proportional to $1/x^4$.

pdf = ProbabilityDistribution[
a^3 b^3/(b^3 - a^3)/x^4, {x, a, b}] /. {a -> 0.01, b -> 0.5};
RandomVariate[pdf, 10000]; // Timing
RandomVariate[pdf]; // Timing


To generate 10,000 random variables according to this distribution or to generate 1 usually takes about the same amount of time (about 2 seconds); about half the time, 10,000 samples is faster. I understand that there's an amazing speedup that happens because of vectorized manipulations, but this is a little ridiculous. I believe that the expensive piece here is computing the inverse CDF when I call RandomVariate; the actual sampling part is quite quick (hence why it doesn't care about the 10,000 number). Is there some way to store the result of the inverse CDF computation, so I can generate appropriately sampled numbers from the same PDF quickly later on?

• My apologies on the normalization - must've had a copy-paste error from when I computed it or something... Simply fixing the normalization by inserting a factor of 3 definitely speeds up the process! – Jolyon Aug 16 '18 at 15:37

This is covered in Chapter 2 of our book, "Mathematical Statistics with Mathematica" (section 2.6). A free copy of the whole text (or chapter) can be downloaded from:

http://www.mathstatica.com/book/bookcontents.html

Note that your pdf is not well-defined as it does not integrate to 1. Making that correction, your pdf has form:

Find the cdf $P(X<x)$:

where I am using the Prob function from the mathStatica package for Mathematica (or use Mathematica's CDF function, or just Integrate it yourself). Then, the inverse cdf is:

The first solution (the real one) is the one that interests us. For your case, the inverse cdf is thus:

blah = (a b)/(b^3 + a^3 u - b^3 u)^(1/3) /. {a -> 0.01, b -> 0.5};


Then, a single pseudorandom number can be generated with:

 blah /. u -> RandomReal[]; // AbsoluteTiming


{0.000035, Null}

The latter is some 200,000 times faster on my Mac than the OP's version.

Similarly, to generate 100,000 pseudorandom values:

 blah /. u -> RandomReal[{}, 100000]; // AbsoluteTiming


{0.002089, Null}

• This is definitely the fastest option. As Anjan Kumar points out, after n~10, the vectorized nature of RandomVariate beats @mef's solution. Computing the appropriate root explicitly as you've done wins hands down -- I find it to be about a factor of 20 faster than RandomVariate (with the correct normalization). – Jolyon Aug 16 '18 at 15:50

First of all your "pdf" does not integrate to one. (It integrates to 1/3.) Let's define a distribution that does integrate to one:

myDist[a_, b_] = ProbabilityDistribution[3 a^3 b^3/(b^3 - a^3)/x^4, {x, a, b}]


We can generate a random variate from this distribution:

RandomVariate[myDist[.01, .5]]


On my computer this is about 35 times faster than when using your "pdf". But that can be speeded up dramatically by computing the inverse CDF:

Simplify[Quantile[myDist[a, b], q], 0 < q < 1]


Root[-a^3 b^3 + (b^3 + a^3 q - b^3 q) #1^3 &, 1]

We can use this to construct our own random number generator:

randomDraw[a_, b_] :=
With[{q = RandomReal[]},
Root[-a^3 b^3 + (b^3 + a^3 q - b^3 q) #1^3 &, 1]
]


Now

randomDraw[.01,.5]


produces an additional speed up of a factor of about 150. (In other words it's about 5000 times faster than using "pdf".)

• Comparing the timings of RandomVariate[myDist[.01, .5], 10000]; // AbsoluteTiming and Table[randomDraw[0.01, 0.5], 10000]; // AbsoluteTiming, RandomVariate is at-least 5 times faster. Looks like for a length greater than 10, RandomVariate is superior. Perhaps, there is a better way to generate points without Table I'm not aware of. – Anjan Kumar Aug 16 '18 at 14:12