How to generate random numbers from a custom distribution?

[John Shannon]

I'm writing a maximum likelihood evaluator and I want to test that it works by using it on data drawn randomly from a distribution with known parameters. But I'm having trouble finding a way to generate a random variable from the distribution I want (I'm working with a variation of the Weibull distribution). Is there a way to use the existing random number functions in Stata to build a random number generator for a different distribution?

[Nick Cox]

Yes. the most general way is to pass uniform randoms to an expression for the quantile function (some people say inverse cumulative distribution function).

Example at http://stackoverflow.com/questions/4...ution-sampling

[Christophe Kolodziejczyk]

Kenneth Train in his book Discrete Choice Methods with Simulation has a very useful chapter on how to draw from densities. Here is a link to the chapter: http://eml.berkeley.edu/choice2/ch9.pdf
Other chapters are also available if you are interested .

See also these entries on the stata-blog
http://blog.stata.com/2016/03/10/how...bers-in-stata/
http://blog.stata.com/2012/07/18/usi...rators-part-1/

[John Shannon]

That's very helpful, thank you. To make sure I understand, I tested out the concept by writing some code. The idea (using Trainer's notation) is to exploit the invertibility of a monotone univariate cdf:
\[ F(\varepsilon)=\mu\Leftrightarrow\varepsilon=F^{-1}(\mu) \]
You can simulate a distribution by randomly drawing mu from the uniform distribution on [0,1]. Using the example of a weibull distribution, this becomes:
\[ \mu = 1 - \exp\left[-\left(\frac{\varepsilon}{b}\right )^a\right]
\Leftrightarrow
\varepsilon = b\left[-\ln(1-\mu)\right]^{1/a} \]
To simulate a weibull variable in Stata, I coded:

clear
set obs 10
gen p = runiform()
gen weibull1 = invweibull(1,2,p)
gen weibull2 = 2*(-ln(1-p))^1
list weibull1 weibull2

This produces identical variables. Interestingly, rweibull(1,2) should produce the same results as invweibull(1,2,p) because that's how the manual says rweibull() works. However, the random numbers in rweibull() must be generated differently than in runiform() since I get different estimates when I type:

clear
set obs 10
local state = c(rngstate)
gen p = runiform()
set rngstate `state'
gen weibull1 = rweibull(1,2)
gen weibull2 = invweibull(1,2,p)

[Kreshna Gopal (StataCorp)]

John Shannon makes a good point: in his code above, rweibull(1,2) is expected to produce the same result as invweibull(1,2,p), but it does not.

However, the reason for the difference is not because runiform() and rweibull() generate different random numbers. The random numbers are identical. The reason for the difference is that we use invweibulltail() rather than invweibull(). We chose invweibulltail() because it is numerically more stable.

In the following code, we can see that rweilbull(1,2) produces the same result as invweibulltail(1,2,p), which is the same as invweibull(1,2,1-p).

clear
set obs 10
local state = c(rngstate)
gen p = runiform()
set rngstate `state'
gen weibull1 = rweibull(1,2)
gen weibull2 = invweibull(1,2,p)
gen weibull3 = invweibulltail(1,2,p)
gen weibull4 = invweibull(1,2,1-p)

With either invweibull() or invweibulltail(), we correctly generate random variates from the Weibull distribution.

Thank you John Shannon for bringing this up. Our documentation will be changed in a future update.

-- Kreshna