Optimize RandomVariate for a $1/x^4$ distribution
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".)
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}