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

Question

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?

Answer 1

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}

Answer 2

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".)