3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Jolyon Aug 16 '18 at 15:37
2
$\begingroup$

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:

enter image description here

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

enter image description here

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:

enter image description here

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}

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – Jolyon Aug 16 '18 at 15:50
5
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Anjan Kumar Aug 16 '18 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.