Sidney Richard Percy (1822-1886) was an English landscape painter.
We have sections
\begin{equation}{\label{a}}\tag{A}\mbox{}\end{equation}
The Inverse Transform Method.
Suppose we want to generate the value of a discrete random variable \(X\) having probability mass function
\[\mathbb{P}\{X=x_{j}\}=p_{j}, j=0,1,2,\cdots ,\mbox{ where }\sum_{j} p_{j}=1.\]
To accomplish this, we generate a random number \(u\), and set
\[x=\left\{\begin{array}{ll}
x_{0} & \mbox{if \(u<p_{0}\)}\\
x_{1} & \mbox{if \(p_{0}\leq u<p_{0}+p_{1}\)}\\
\vdots &\\
x_{j} & \mbox{if \(\sum_{i=1}^{j-1} p_{i}\leq u\leq\sum_{i=1}^{j}p_{i}\)}\\
\vdots &
\end{array}\right .\]
For \(0<a<b<1\), since \(\mathbb{P}\{a\leq U<b\}=b-a\), we have
\[\mathbb{P}\{X=x_{j}\}=\mathbb{P}\left\{\sum_{i=1}^{j-1}p_{i}\leq U<\sum_{i=1}^{j}p_{i}\right\}=p_{j}\]
and so \(X\) has the desired distribution.
The preceding discussions can be written algorithmically as
\begin{equation}{\label{ch2eq1}}\tag{1}
\begin{array}{l}
\mbox{Generate a random number \(u\)}\\
\mbox{If \(u<p_{0}\) set \(x=x_{0}\) and stop}\\
\mbox{If \(u<p_{0}+p_{1}\) set \(x=x_{1}\) and stop}\\
\mbox{If \(u<p_{0}+p_{1}+p_{2}\) set \(x=x_{2}\) and stop}\\
\vdots
\end{array}
\end{equation}
Let the values \(x_{i}\) be ordered satisfying \(x_{0}<x_{1}<x_{2}<\cdots\), and let \(F\) denote the distribution function of \(X\). Then, we have
\begin{align*} F(x_{k}) & =\mathbb{P}\{X\leq x_{k}\}\\ & =\sum_{i=0}^{k}\mathbb{P}\{X=x_{i}\}\\ & =\sum_{i=1}^{k}p_{i},\end{align*}
which says that \(x=x_{j}\) when
\[\sum_{i=1}^{j-1}p_{i}=F(x_{j-1})\leq u<F(x_{j})=\sum_{i=1}^{j}p_{i}.\]
In other words, after generating a random number \(u\) we determine the value of \(x\) by finding the interval \([F(x_{j-1}),F(x_{j})]\) in which \(u\) lies or equivalently, by finding the inverse of \(F(u)\). For this reason, the above method is called the discrete inverse transform method for generating \(X\). Furthermore, we have the following algorithm
\begin{equation}{\label{sch2eq2}}\tag{2}
\begin{array}{l}
\mbox{Generate a random number \(u\)}\\
\mbox{If \(u<F(x_{0})\) set \(x=x_{0}\) and stop}\\
\mbox{If \(u<F(x_{1})\) set \(x=x_{1}\) and stop}\\
\mbox{If \(u<F(x_{2})\) set \(x=x_{2}\) and stop}\\
\vdots
\end{array}
\end{equation}
The amount of time it takes to generate a discrete random variable by the above method is proportional to the number of intervals one must search. For this reason it is sometimes worthwhile to consider the possible values \(x_{j}\) of \(X\) in decreasing order of the \(p_{j}\).
Example. If we want to simulate a random variable \(X\) satisfying
\[p_{1}=0.2, p_{2}=0.15, p_{3}=0.25, p_{4}=0.4,\mbox{ where }p_{j}=\mathbb{P}\{X=j\}\]
then we could generate \(u\) and do the following
\[\begin{array}{l}
\mbox{If \(u<0.2\) set \(x=1\) and stop}\\
\mbox{If \(u<0.35\) set \(x=2\) and stop}\\
\mbox{If \(u<0.6\) set \(x=3\) and stop}\\
\mbox{otherwise set \(x=4\)}
\end{array}\]
However, a more efficient procedure is the following
\[\begin{array}{l}
\mbox{If \(u<0.4\) set \(x=4\) and stop}\\
\mbox{If \(u<0.65\) set \(x=3\) and stop}\\
\mbox{If \(u<0.85\) set \(x=1\) and stop}\\
\mbox{otherwise set \(x=2\)}
\end{array}\]
The following two examples show that it can be generated without needing to search for the relevant interval.
\begin{equation}{\label{sch2ex2}}\tag{3}\mbox{}\end{equation}
Example \ref{sch2ex2}. Suppose that we want to generate the value of \(X\) which is equally likely to take on any of the values \(1,2,\cdots ,n\). That is, we have \(\mathbb{P}\{X=j\}=1/n\) for \(j=1,2,\cdots ,n\). Using the above results, it follows that we can accomplish this by generating \(u\) and setting
\[x=j\mbox{ if }\frac{j-1}{n}\leq u<\frac{j}{n}.\]
Therefore, \(x\) will equal \(j\) when \(j-1\leq nu<j\). Equivalently, we have \(x=[nu]+1\), where \([x]\) is the integral part of \(x\). \(\sharp\)
Example. Recall that \(X\) is said to be a geometric random variable with parameter \(p\) when we have
\[\mathbb{P}\{X=i\}=pq^{i-1}\]
for \(i\geq 1\), where \(q=1-p\). The random variable \(X\) can be thought of as representing the time of the first success when independent trials, each of which is a success with probability \(p\), are performed. Since
\begin{align*}
\sum_{i=1}^{j-1} P\{X=i\} & =\mathbb{P}\{X\leq j-1\}=1-P\{X>j-1\}\\
& =1-\mathbb{P}\{\mbox{first \(j-1\) trials are all failures}\}=1-q^{j-1}\mbox{ for }j\geq 1
\end{align*}
we can generate the value of \(x\) by generating a random number \(u\) and setting \(x\) equal to that value \(j\) for which
\[1-q^{j-1}\leq u<1-q^{j}\]
or equivalently, for which
\[q^{j}<1-u\leq q^{j-1}\]
In other words, we can define \(x\) by
\[x=\min\{j:q^{j}<1-u\}\mbox{ since }0<q<1.\]
Equivalently, the value \(x\) can be expressed as
\begin{align*} x & =\min\{j:j\ln q<\ln (1-u)\}\\ & =\min\left\{j:j>\frac{\ln (1-u)}{\ln q}\right\}\end{align*}
since \(\ln q\) is negative ($0<q<1$). Therefore, we can express \(x\) as
\[x=\left [\frac{\ln (1-u)}{\ln q}\right ]+1.\]
Finally, by noting that \(1-U\) is also uniformly distributed on \((0,1)\), it follows
\[X=\left [\frac{\ln U}{\ln q}\right ]+1\]
is also a geometric with parameter \(p\). \(\sharp\)
Example. A random variable \(X\) that takes on one of the values \(0,1,2,\cdots\) is said to be a Poisson random variable with parameter \(\lambda\) when its probability mass function is given by
\[p_{i}=\mathbb{P}\{X=i\}=e^{-\lambda}\frac{\lambda^{i}}{i!}\]
for \(i=0,1,2,\cdots\). To compute the Poisson probabilities, we make use of the following recursive formula
\begin{align*} \frac{p_{i+1}}{p_{i}} & =\frac{e^{-\lambda}\lambda^{i+1}/(i+1)!}{e^{-\lambda}\lambda^{i}/i!}\\ & =\frac{\lambda}{i+1}\end{align*}
Equivalently, we have
\[p_{i+1}=\frac{\lambda}{i+1}p_{i}.\]
The inverse transform algorithm for generating a Poisson random variable with mean \(\lambda\) can be expressed as follows by referring to (\ref{ch2eq1}) and (\ref{sch2eq2}). The quantity \(i\) refers to the value presently under consideration; \(p=p_{i}\) is the probability that \(X\) equals \(i\), and \(F=F(i)\) is the probability that \(X\) is less than or equal to \(i\)
\[\begin{array}{l}
\mbox{Step 1. Generate a random number \(u\).}\\
\mbox{Step 2. \(i=0\), \(p=e^{-\lambda}=p_{0}\), \(F=p\).}\\
\mbox{Step 3. If \(u<F\) set \(X=i\) and stop}\\
\mbox{Step 4. \(p=\lambda p/(i+1)\), \(F\leftarrow F+p\), \(i\leftarrow i+1\).}\\
\mbox{Step 5. Go to Step 3.}
\end{array}\]
Example. Suppose that we want to generate the value of a binomial \(B(n,p)\) random variable \(X\), i.e.,
\[\mathbb{P}\{X=i\}=\frac{n!}{i!(n-i)!}p^{i}(1-p)^{n-i}\]
for \(i=0,1,2,\cdots ,n\). The following recursive formula expressing \(p_{i+1}\) in terms of \(p_{i}\) is useful
\begin{align*}
p_{i+1} & =\frac{n!}{(n-i-1)!(i+1)!}p^{i+1}(1-p)^{n-i-1}\\
& =\frac{n!(n-i)}{(n-i)!i!(i+1)}p^{i}(1-p)^{n-i}\frac{p}{1-p}\\
& =\frac{n-i}{i+1}\cdot\frac{p}{1-p}\cdot p_{i}.
\end{align*}
With \(i\) denoting the value currently under consideration, \(pr=\mathbb{P}\{X=i\}\) the probability that \(X\) is equal to \(i\), and \(F=F(i)\) the probability that \(X\) is less than or equal to \(i\), the algorithm can be expressed as follows
\[\begin{array}{l}
\mbox{Step 1. Generate a random number \(u\)}\\
\mbox{Step 2. \(c=p/(1-p)\), \(i=0\), \(pr=(1-p)^{n}\), \(F=pr\).}\\
\mbox{Step 3. If \(u<F\), set \(X=i\) and stop.}\\
\mbox{Step 4. \(pr\leftarrow [c(n-i)/(i+1)]pr\), \(F\leftarrow F+pr\), \(i\leftarrow i+1\).}\\
\mbox{Step 5. Go to Step 3.}
\end{array}\]
\begin{equation}{\label{b}}\tag{B}\mbox{}\end{equation}
The Acceptance-Rejection Technique.
Suppose that we have an efficient method for simulating a random variable having probability mass function \(\{q_{j}\}\). We can use this as the basis for simulating from the distribution having mass function \(\{p_{j}\}\) by first simulating a random variable \(Y\) having mass function \(\{q_{j}\}\) and then accepting this simulated value with a probability proportional to \(p_{j}/q_{j}\). Specifically, let \(c\) be a constant satisfying
\[\frac{p_{j}}{q_{j}}\leq c\mbox{ for all \(j\) such that \(p_{j}>0\)}\]
(in order to make \(p_{j}/cq_{j}\leq 1\) for all \(j\)). We shall present the acceptance-rejection algorithm for simulating a random variable \(X\) having mass function \(p_{j}=P\{X=j\}\).
\[\begin{array}{l}
\mbox{Step 1. Simulate the value of \(Y=j\) having probability mass function \(q_{j}\).}\\
\mbox{Step 2. Generate a random number \(u\)}.\\
\mbox{Step 3. If \(u<p_{j}/cq_{j}\), set \(X=j\) and stop. Otherwise go to Step 1.}
\end{array}\]
Remark A. Note that the way in which we “accept the value \(Y=j\) with probability \(p_{j}/cq_{j}\)” is by generating a random number \(u\) and then accepting \(Y=j\) and setting \(X=j\) when \(u\leq p_{j}/cq_{j}\). \(\sharp\)
Theorem. The acceptance-rejection algorithm generates a random variable \(X\) satisfying
\[\mathbb{P}\{X=j\}=p_{j}\mbox{ for }j=0,1,2,\cdots .\]
In addition, the number of iterations of the algorithm needed to obtain \(X\) is a geometric random variable with mean \(c\).
Proof. Let us determine the probability that a single iteration produces the accepted value \(j\).
\begin{align*}
\mathbb{P}\{Y=j,\mbox{it is accepted}\} & =\mathbb{P}\{Y=j\}\cdot \mathbb{P}\{\mbox{accept}|Y=j\}\\
& =q_{j}\cdot\frac{p_{j}}{cq_{j}}\mbox{ (by Step 3 and Remark A)}=\frac{p_{j}}{c}.
\end{align*}
Summing over \(j\) yields the probability that a generated random variable is accepted
\[\mathbb{P}\{\mbox{accepted}\}=\sum_{j}\frac{p_{j}}{c}=\frac{1}{c}.\]
As each iteration independently results in an accepted value with probability \(1/c\), we see that the number of iterations needed is geometric with mean \(c\). Also, we have
\begin{align*} \mathbb{P}\{X=j\} & =\sum_{n} \mathbb{P}\{\mbox{accepted on iteration \(n\)}\}\\ & =\sum_{n}
\left (1-\frac{1}{c}\right )^{n-1}\cdot\frac{p_{j}}{c}=p_{j}.\end{align*}
This completes the proof. \(\blacksquare\)
Example. Suppose that we want to simulate the value of a random variable \(X\) that takes one of the values \(1,2,\cdots ,10\) with respective probability \(0.11, 0.12,0.09, 0.08, 0.12, 0.1, 0.09, 0.09, 0.1, 0.1\), Whereas one possibility is to use the inverse transform algorithm, a better approach is to use the acceptance-rejection method with \(q\) being the discrete uniform density on \(1,2,\cdots ,10\). That is, \(q_{j}=1/10\), \(j=1,2,\cdots ,10\). For this choice of \(\{q_{j}\}\), we can choose \(c\) by
\[c=\max\frac{p_{j}}{q_{j}}=1.2\]
Therefore, the algorithm would be as follows
\[\begin{array}{l}
\mbox{Step 1. Generate a random number \(u_{1}\) and set \(Y=j=[10u_{1}]+1\) according Example 3}\\
\mbox{Step 2. Generate a second random number \(u_{2}\).}\\
\mbox{Step 3. If \(u_{2}\leq p_{j}/cq_{j}=p_{j}/1.2\cdot 0.1=p_{y}/0.12\), set \(X=j\) and stop. Otherwise go to Step 1.}
\end{array}\]
On average, this algorithm requires only \(1.2\) iterations to obtain the generated value of \(X=j\). \(\sharp\)
\begin{equation}{\label{c}}\tag{C}\mbox{}\end{equation}
The Composition Approach.
Suppose that we have an efficient method to simulate the value of a random variable having either of the two probability mass functions \(\{p_{1j}\}\) and \(\{p_{2j}\}\), and that we want to simulate the value of the random variable \(X\) having mass function
\begin{equation}{\label {sch2eq12}}\tag{4}
\mathbb{P}\{X=j\}=\alpha p_{1j}+(1-\alpha )p_{2j}
\end{equation}
where \(0<\alpha <1\). One way to simulate such a random variable \(X\) is to note that if \(X_{1}\) and \(X_{2}\) are random variables having respective mass functions \(\{p_{1j}\}\) and \(\{p_{2j}\}\), then the random variable \(X\) defined by
\[X=\left\{\begin{array}{ll}
X_{1} & \mbox{with probability \(\alpha\)}\\
X_{2} & \mbox{with probability \(1-\alpha\)}
\end{array}\right .\]
will have its mass function given by (\ref{sch2eq12}). According to the inverse transform method, it follows that we can generate the value of such a random variable by first generating a random number \(u\) and then generating a value of \(X=X_{1}=x_{1}\) if \(u<\alpha\) and of \(X=X_{2}=x_{2}\) if \(u>\alpha\).
Example. Suppose that we want to generate the value of a random variable \(X\) satisfying
\[p_{j}=\mathbb{P}\{X=j\}=\left\{\begin{array}{ll}
0.05 & \mbox{for \(j=1,2,3,4,5\)}\\
0.15 & \mbox{for \(j=6,7,8,9,10\)}
\end{array}\right .\]
Since \(p_{j}=0.5p_{1j}+0.5p_{2j}\), where \(p_{1j}=0.1\) for \(j=1,2,\cdots ,10\) and
\[p_{2j}=\left\{\begin{array}{ll}
0 & \mbox{for \(j=1,2,3,4,5\)}\\
0.2 & \mbox{for \(j=6,7,8,9,10\)}
\end{array}\right .\]
we can accomplish this by first generating a random number \(u\) and then generating from the discrete uniform over \(1,2,\cdots ,10\) if \(u<0.5\) and from the discrete uniform over \(6,7,8,9,10\) if \(u>0.5\). That is, we can simulate \(X\) as follows
\[\begin{array}{l}
\mbox{Step 1. Generate a random number \(u_{1}\)}\\
\mbox{Step 2. Generate a random number \(u_{2}\)}\\
\mbox{Step 3. If \(u_{1}<0.5\), set \(x=[10u_{2}]+1\). Otherwise, set \(x=[5u_{2}]+6\) by referring to Example 3.}
\end{array}\]
If \(F_{i}\) for \(i=1,2,\cdots ,n\) are distribution functions and \(\alpha_{i}\) for \(i=1,2,\cdots ,n\) are nonnegative numbers summing to \(1\), then the distribution function given by
\[F(x)=\sum_{i=1}^{n} \alpha_{i}F_{i}(x)\]
is said to be a mixture or a composition of the distribution functions \(F_{i}\) for \(i=1,2,\cdots ,n\). One way to simulate from \(F\) is to
simulate a random variable \(I=i\) with probability \(\alpha_{i}\) for \(i=1,2,\cdots ,n\), and then to simulate from the distribution \(F_{I}\). In other words, if the simulated value of \(I\) is \(I=j\), then the second simulation is from \(F_{j}\). This approach to simulating from \(F\) is often referred to as the composition method.
