Updated 2017-11-05
Given an uniform random variable, how do we generate random numbers with a specific probability distribution?
To generate discrete random variables with distribution \(\mathbb P(X=x_i)=p_i\), for \(i=1,\dotsc,n\), we follow the following algorithm.
To generate continuous random variables with distribution \(F_X(x)\), follow these steps:
“Pseudo” Proof: Suppose we have the uniform random variable \(Y\sim U(0,1)\) and we use the inverse distribution function \(F_X^{-1}\) to obtain random \(X\). Note that the distribution function for \(Y\) is, by definition, \(P_Y(Y<y)=y\).
\[\begin{align} \mathbb P(X\leq x)&=\mathbb P(F_X^{-1}(Y)\leq x)\\ &=\mathbb P(Y\leq F_X(x))\\ &=F_X(x) \end{align}\]Which is what we wanted.
Note: This method does not work directly on Gaussian random variables, use Polar Algorithm for Gaussian random variables instead.
As noted above, we can’t use inverse transformation method for Gaussian random variables. In order to simulate continuous random variable to have the distribution \(F_X(x)=N(0,1)\), we need to follow these steps:
To obtain a random variable with non-standard normal distribution \(Z\sim N(\mu, \sigma^2)\), we use the linear properties of the Gaussian random variables.
Use this method if the distribution function \(F_X\) is composed of a sum of other distribution functions \(F_{X_i}\):
\[F_X(x)=\sum_{i=1}^{n}p_i F_{X_i}(x)\qquad\sum_{i=1}^n p_i=1,\qquad p_i\geq 0\]Where \(F_{X_i}\) are CDFs.
To get a random variable with distribution \(F_X\), follow these steps:
We can use this method if it’s hard to use previous methods to generate a random variable \(X\) with distribution \(F_X\).
So we find a proposal distribution \(F_Y\) where it is easy to sample. Note that we must ensure the range of the density functions, \(f_X\) and \(f_Y\) are the same.
Next, we need to know the upper bound \(a\), where \(a\geq \frac{f_X(x)}{f_Y(x)}\).
With all that in mind, we may begin the algorithm:
Note: The probability of \(Y\) being accepted is \(\frac{1}{a}\), therefore we should always choose \(a\) to be as small as possible given \(a\geq 1\).
Example: half normal distribution
The half normal distribution function has the density \(f_X=\frac{2}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\). The inverse CDF is very difficult to find, so we use the A-R method.
The range of the density of this distribution spans \([0,\infty)\), so it would be appropriate to choose the exponential random variable distribution as our proposal distribution. The density of the proposal distribution is \(f_Y(x)=e^{-x}\).
Thus
\[\frac{f_X(x)}{f_Y(x)}=\frac{2}{\sqrt{2\pi}}e^{-\frac{x^2}{2}+x}\]Taking the derivative and the second derivative, we find the maximum is at \(x=1\). Thus, our upper bound is
\[\frac{f_X(1)}{f_Y(1)}=\frac{2}{\sqrt{2\pi}}e^\frac12\approx1.3=a\]Now we apply the A-R algorithm. First we generate an exponential random variable, \(Y\sim f_Y\).
Next, we generate \(Z\sim U(0,1)\).
Lastly, we set \(X=Y\) if \(Z\leq \frac{\frac{2}{\sqrt{2\pi}}e^{y-\frac{y^2}{2}}}{1.3}\)
Suppose we have a vector consists of two random variables where \((x, y)\sim F_{X,Y}(x,y)\).
If \(X\) and \(Y\) are independent. Great! Then we can just generate each of them independently.
If not, then the algorithm is as follows:
Express the joint PDF of the two random variables as a conditional PDF
\[f_{X,Y}(x,y)=f_{y\vert x}(y\vert x)f_X(x)\]For \(i=1,2,\dotsc,n\) do:
Example:
Suppose \(f_{X,Y}(x,y)=x+y\) and \(0\leq x,y \leq 1\), then the marginal PDF for \(X\) is
\[f_X(x)=\int_0^1 f_{X,Y}(x,y)\mathrm dy=x+\frac12\]Next, using the marginal PDF for \(X\), we can find the conditional PDF:
\[f_{Y\vert X}(y\vert x)=\frac{f_{X,Y}(x,y)}{f_X(x)}=\frac{x_i+y}{x_i+\frac{1}{2}}\]Following the rest of the algorithm, we sample \(x_i\sim x+\frac12\) and sample \(y_i\sim \frac{x_i + y}{x_i + \frac12}\) for all \(i\)’s.