In Cooperation with:

American Society for Quality Statistics Division

American Statistical Association

Bernoulli Society for Mathematical Statistics and Probability

Institute of Mathematical Statistics

International Biometric Society

International Chinese Statistical Association

International Society for Bayesian Analysis

International Statistical Institute

Royal Statistical Society

Statistical Society of Canada / Société statistique du Canada

Generating Random Variables

ref

Generating Random Variables

Christian Robert¹
Université Paris Dauphine and CREST, INSEE

George Casella²
University of Florida

Keywords and Phrases: Random Number Generator, Probability Integral Transform, Accept Reject Algorithm, Metroplis-Hastings Algorithm, Gibbs Sampling, Markov Chains

Introduction

The methods described in this article mostly rely on the possibility of producing (with a computer) a supposedly endless flow of random variables (usually iid) for well-known distributions. Such a simulation is, in turn, based on the production of uniform random variables. There are many ways in which uniform pseudorandom numbers can be generated. For example there is the Kiss algorithm of Marsaglia and Zaman (1993); details on other random number generators can be found in the books of Rubinstein (1981), Ripley (1987), Fishman (1996), and Knuth (1998).

Generating Nonuniform Random Variables

The generation of random variables that are uniform on the interval

, the Uniform

distribution, provides the basic probabilistic representation of randomness. The book by Devroye (1986) is a detailed discussion of methods for generating nonuniform variates, and the subject is one of the many covered in Knuth (1998). Formally, we can generate random variables with any distribution by means of the following.

Probability Integral Transform

For a function

on $\Re$ , the generalized inverse of

is $F^{-}(u) = \inf \; \{x:\; F(x) \geq u\} \;$ . If

is a cumulative distribution function (cdf), $F^{-}(U)$ has the cdf

and to generate $X \sim F$ ,

Generate according to Uniform
Make the transformation .

Generating Discrete Variables

For discrete distributions, if the random variable satisfies

$\displaystyle P(X=k)=p_k, \quad k=0,1,2,\ldots$

then if $U \sim$ uniform

can be generated by

$\displaystyle k=0:$		$\displaystyle U \le p_0 \Rightarrow X=0,$
$\displaystyle k=1, 2, \ldots,:$		$\displaystyle \sum_{j=0}^{k-1} p_j \le U < \sum_{j=0}^{k}p_{j} \Rightarrow X=k.$

Mixture Distributions

In a mixture representation a density has the form

$\displaystyle f(x)$	$\displaystyle =$	$\displaystyle \int_{\cal{Y}} \; f(x\vert y) g(y)\; dy$ (continuous)
$\displaystyle f(x)$	$\displaystyle =$	$\displaystyle \sum_{i \in {\cal Y}} \; p_{i} \; f_{i}(x)\;$ (discrete).

We then can simulate

1.: $Y \sim g(y)$ , $X \sim f(x\vert Y=y)$ , or
2.: , $X \sim f_i(x)$

Accept-Reject Methods

Another class of methods only requires the form of the density of interest - called the target density. We simulate from a density , called the candidate density.

Given a target density , we need a candidate density and a constant such that

$\displaystyle f(x) \leq M g(x)$

on the support of

Accept-Reject Algorithm:

To produce a variable

distributed according to

1.: Generate $Y \sim g$ , $U\sim$ Uniform ;
2.: Accept if $U \leq f(Y) / Mg(Y)$ ;
3.: Otherwise, return to 1.

Notes:

(a).: The densities and need be known only up to a multiplicative factor.
(b).: The probability of acceptance is , when evaluated for the normalized densities.

Envelope Accept-Reject Algorithm:

If the target density

is difficult to evaluate, the Envelope Accept-Reject Algorithm (called the squeeze principle by Marsaglia 1977) may be appropriate.

If there exist a density $g_{m}$ , a function $g_{l}$ and a constant such that

$\displaystyle g_{l}(x) \leq f(x) \leq Mg_{m}(x) \;,$

then the algorithm

1.: Generate $X \sim g_{m}(x)$ , $\; U \sim$ Uniform ;
2.: Accept if $U \leq g_{l}(X) / M g_{m}(X)$ ;
3.: Otherwise, accept if $U \leq f(X) /M g_{m}(X)$
4.: Otherwise, return to 1.

produces $X \sim f$ . The number of evaluations of

is potentially decreased by a factor ${1\over M}$ . If

is log-concave, Gilks and Wild (1992) construct a generic accept-reject algorithm that can be quite efficient.

Markov Chain Methods

Every simulation method discussed thus far has produced independent random variables whose distribution is exactly the target distribution. In contrast, Markov chain methods produce a sequence of dependent random variables whose distribution converges to the target. Their advantage is their applicability in complex situations.

Recall that a sequence $X_0, X_1, \ldots, X_n,\allowbreak \ldots$ of random variables is a Markov chain if, for any , the conditional distribution of given $x_{t-1},x_{t-2},\ldots, x_0$ is the same as the distribution of given $x_{t-1}$ ; that is,

$\displaystyle P(X_{k+1} \in A \vert x_0, x_1, x_2, \ldots, x_k)= P(X_{k+1} \in A \vert x_k).$

Metropolis - Hastings Algorithm

The Metropolis-Hastings (M-H) algorithm (Metropolis et al. 1953, Hastings 1970) associated with the target density and the candidate density $q(y\vert x)$ produces a Markov chain $(X^{(t)})$ through

Metropolis -Hastings

Given $x^{(t)}$ ,

Generate $\; Y_t \sim q(y\vert x^{(t)})$ .
Take
$\displaystyle X^{(t+1)} =\left\{ \begin{array}{ll} Y_t & \mbox{ with prob. } \r... ...t) \ x^{(t)} & \mbox{ with prob. } \ 1 - \rho(x^{(t)},Y_t)\end{array}\right.$
where
$\displaystyle \rho(x,y) = \min \left\{ {f(y)\over f(x)} \; {q(x\vert y)\over q(y\vert x)} \;, 1 \right\} \;,$
Then $X^{(t)}$ converges in distribution to $X \sim f$ .

Notes:

(a).: For $q(\cdot\vert x)=q(\cdot)$ we have the independent M-H algorithm, and for $q(x\vert y) = q(y\vert x)$ we have a symmetric M-H algorithm, where $\rho$ does not depend on . Also, $q(x\vert y) = q(y-x)$ , symmetric around zero, is a random walk M-H algorithm.
(b).: Like the Accept-Reject method, the Metropolis - Hastings algorithm only requires knowing and up to normalizing constants.

The Gibbs Sampler

For , write the random variable ${\mathbf X} \in {\cal X}$ as ${\mathbf X}=(X_1,\ldots,X_p) \sim f$ , where the 's are either uni- or multidimensional, with conditional distributions

$\displaystyle X_i\vert x_1, x_2, \ldots, x_{i-1}, x_{i+1}, \ldots,x_p \sim f_i(x_i\vert x_1, x_2, \ldots, x_{i-1}, x_{i+1}, \ldots,x_p)$

for $i=1,2,\ldots,p$ . The associated Gibbs sampler is given by the following transition from $X^{(t)}$ to $X^{(t+1)}$ :

The Gibbs sampler

Given ${\bf x}^{(t)} = (x_1^{(t)},\ldots,x_p^{(t)})$ , generate

1.: $X_1^{(t+1)} \sim f_1(x_1\vert x_2^{(t)},\ldots,x_p^{(t)})$ ;
2.: $X_2^{(t+1)} \sim f_2(x_2\vert x_1^{(t+1)},x_3^{(t)},\ldots,x_p^{(t)})$ ,; $\qquad\vdots$
p.: $X_p^{(t+1)} \sim f_p(x_p\vert x_1^{(t+1)},\ldots,x_{p-1}^{(t+1)})$ ,

then $X^{(t)}$ converges in distribution to $X \sim f$ .

Notes:

(a).: The densities $f_1,\ldots,f_p$ are called the full univariate conditionals.
(b).: Even in a high dimensional problem, all of the simulations can be univariate.
(c).: The Gibbs sampler is, formally, a special case of the M-H algorithm (see Robert and Casella 2004, Section ) but with acceptance rate equal to .

The Slice Sampler

A particular version of the Gibbs sampler, called the slice sampler (Besag and Green 1993, Damien et al. 1999), can sometimes be useful. Write $f(x)= \prod_{i=1}^k f_i(x)$ where the

's are positive functions, not necessarily densities

Slice sampler

Simulate

1.: $\omega_1^{(t+1)} \sim$ Uniform $_{[0,f_1(x^{(t)})]}$ ;; $\qquad\ldots$
k.: $\omega_k^{(t+1)} \sim$ Uniform $_{[0,f_k(x^{(t)})]}$ ;
k+1.: $x^{(t+1)} \sim$ Uniform $_{A^{(t+1)}}$ , with
$\displaystyle A^{(t+1)} = \{ y;\; f_i(y) \ge \omega_i^{(t+1)},\ i=1,\ldots,k \}.$

then $X^{(t)}$ converges in distribution to $X \sim f$ .

Application

In this section we give some guidelines for simulating different distributions. See Robert and Casella (2004, 2010) for more detailed explanations, and Devroye (1986) for more algorithms. In what follows,

is Uniform

unless otherwise specified.

Arcsine

$f(x ) = \frac{1}{\pi\sqrt{x(1-x)}},\quad 0 \leq x \leq 1, \quad \sin^2\left(\pi U/2 \right) \sim f$ .

Beta

$f(x) = {\Gamma(r+s) \over{\Gamma(r)\Gamma(s )}}x^{r-1}(1-x)^{s -1}$ , $0 \leq x \leq 1, r > 0, s > 0$ , $\frac{X_1}{X_1+X_2} \sim f$ , where $X_1 \sim$ Gamma

and $X_2 \sim$ Gamma

, independent.

Cauchy $(\mu,\sigma )$

$f(x\vert \mu,\sigma ) = {1\over{\sigma \pi }}{1\over{1+\frac{(x-\mu)^2}{\sigma^2}}}$ , $-\infty < x < \infty$ , $\sigma \tan\left(\frac{\pi}{2}(2U-1) \right) + \mu \sim$ Cauchy $(\mu,\sigma )$ .

If $U \sim$ uniform ${[-\pi/2,\pi/2]}$ , then $\tan(U) \sim$ Cauchy . Also, $X/Y \sim$ Cauchy , where $X, Y, \sim N(0,1)$ , independent.

Chi squared( )

$f(x\vert p) = {1\over{\Gamma (p/2)2^{p/2}}}x^{(p/2)-1}e^{-x/2}$ , $0 \leq x < \infty, p = 1, 2, \dots$ .

$\displaystyle - 2 \sum_{j=1}^\nu \log(U_j ) \sim \chi^2_ {2 \nu},$ or $\displaystyle \chi^2_ {2 \nu} + Z^2 \sim \chi^2_ {2 \nu+1},$

where

is an independent Normal

random variable.

Double exponential $(\mu,\sigma )$

$f(x\vert \mu,\sigma ) = {1\over{2\beta }}e^{-{\left\vert{x-\mu }\right\vert}/\beta }$ , $-\infty < x < \infty$ , $-\infty < \mu < \infty$ , $\beta > 0$ . If $Y \sim$ Exponential $(\beta )$ , then attaching a random sign (

otherwise) gives $X \sim$ Double exponential

, and $\sigma X +\mu \sim$ Double exponential $(\mu,\sigma)$ . This is also known as the Laplace distribution.

Exponential $(\beta )$

$f(x\vert \beta ) = {1\over{\beta }}e^{-x/\beta },\quad 0 \leq x < \infty,\quad \beta > 0$ , $-\beta \log U \sim$ Exponential $(\beta)$ .

Extreme Value $(\alpha, \gamma )$

$f(x\vert \alpha, \gamma ) = e^{-\frac{x-\alpha}{\gamma}-e^{-\frac{x-\alpha}{\gamma}}},\quad \alpha \leq x < \infty, \quad \gamma >0$ . If $X \sim$ Exponential

, $\alpha - \gamma \log (X) \sim$ Extreme Value $(\alpha, \gamma )$ . This is also known as the Gumbel distribution.

F Distribution

$f(x\vert {\nu }_1,\nu_2) = {{\Gamma \left( {{{\nu }_1+\nu_2} \over 2}\right) }... ...t(1 + \left( {{\nu }_1 \over \nu_2 } \right)x\right) ^{({\nu }_1+{ \nu }_2)/2}}$ , $0 \leq x < \infty.$ If $X_1 \sim {\chi }^2_{{\nu }_1}$ and $X_2 \sim {\chi }^2_{{\nu }_2}$ , independent, $\frac{X_1/\nu_1}{X_2/\nu_2} \sim f(x\vert {\nu }_1,\nu_2).$

Gamma $(\alpha,\beta )$

$f(x\vert \alpha,\beta ) = {1\over{\Gamma (\alpha ){\beta }^ {\alpha }}}x^{\alpha -1}e^{-x/\beta }$ , $0 \leq x < \infty$ , $\alpha,\ \beta > 0$ . Then $-\beta \sum^a_{j=1} \log(U_j ) \sim$ Gamma $(a, \beta)$ ,

an integer

If $\alpha$ is not an integer, indirect methods can be used. For example, to generate a Gamma $(\alpha,\beta )$ use Algorithm 3.1 or 4.1 with candidate distribution Gamma , with $a = [\alpha]$ and $b=\beta \alpha/a$ , where $[\alpha]$ is the greatest integer less than $\alpha$ . For the Accept-Reject algorithm the bound on the normalized is $M = \frac{\Gamma(a)}{\Gamma(\alpha)}\frac{\alpha^\alpha}{a^a}e^{-(\alpha-a)}$ . There are many other efficient algorithms.

Note:

If $X \sim$ Gamma $(\alpha,1 )$ then $\beta X \sim$ Gamma $(\alpha,\beta )$ . Some special cases are Exponential

= gamma

, and Chi squared

= gamma

. Also,

has the inverted (or inverse) gamma distribution.

Logistic $(\mu,\beta )$

$f(x\vert \mu,\beta ) = {1\over{\beta }} {{e^{-(x-\mu )/\beta }}\over{{{\left[{1 + e^{-(x-\mu )/\beta }}\right]}}^2}}$ , $-\infty < x < \infty$ , $-\infty < \mu < \infty$ , $\beta >0$ .

$\displaystyle -\beta \log \left(\frac{1-U}{U}\right) + \mu \sim$ Logistic $\displaystyle (\mu,\beta )$

Lognormal $(\mu,{\sigma }^2)$

$f(x\vert \mu,{\sigma }^2) = {1\over{\sqrt{2\pi }\sigma }} {{e^{-(\log x-\mu )^2/(2{\sigma }^2)}}\over x}$ , $0 \leq x < \infty$ , $-\infty < \mu < \infty$ . If $X \sim$ Normal $(\mu,{\sigma }^2)$ .

$\displaystyle e^X \sim$ Lognormal $\displaystyle (\mu,{\sigma }^2)$

Noncentral chi squared $(\lambda,p)$ , $\lambda \ge 0$

$f_p(x\vert\lambda)= \sum_{k=0}^\infty \frac{x^{p/2+k-1}e^{-x/2}}{\Gamma(p/2+k)2^{p/2+k}}\frac{\lambda^k e^{-\lambda}}{k!}$ , $0<x<\infty$ .

$\displaystyle K \sim$ Poisson $\displaystyle (\lambda/2), \quad X \sim \chi_{p+2K}^{2} \Rightarrow X \sim f_p(x\vert\lambda)$

where

is the degrees of freedom and $\lambda$ is the noncentrality parameter. A more efficient algorithm is

$\displaystyle Z\sim\chi_{p-1}^{2}$ and $\displaystyle Y\sim$ N $\displaystyle (\sqrt{\lambda}, 1) \Rightarrow Z+Y^2\sim f_p(x\vert\lambda).$

Normal $(\mu,{\sigma }^2)$

$f(x\vert \mu,{\sigma }^2) = {1\over{\sqrt{2\pi }\sigma }}e^{-(x-\mu )^2/(2{\sigma }^2)}$ , $-\infty < x < \infty$ , $-\infty < \mu < \infty$ , $\sigma > 0$ . The Box-Muller algorithm simulates two normals from two uniforms:

$\displaystyle X_{1} = \sqrt{- 2 \log (U_{1}}) \; \cos (2\pi U_{2})$ and $\displaystyle X_{2} = \sqrt{- 2 \log (U_{1}}) \; \sin (2\pi U_{2}) \;,$

then $X_1, X_2 \sim$ Normal $(\mu,{\sigma }^2)$ .

There are many other ways to generate normal random variables.

(a).: Accept Reject using Cauchy When $f(x) = (1/\sqrt{2\pi})\exp (-x^{2}/2)$ and $g(x) = (1/\pi)1/(1 + x^{2})$ , densities of the normal and Cauchy distributions, respectively, then $f(x)/g(x) \le \sqrt{2\pi/e} = 1.520$ .
(b).: Accept Reject using double exponentialWhen $f(x) = \frac{1}{\sqrt{2\pi}}\exp (-x^{2}/2)$ and $g(x) = (1/2)\allowbreak \exp(-\vert x\vert)$ , $f(x)/g(x) \le \sqrt{2 e/\pi}=1.315$ .
(c).: Slice Sampler
$\displaystyle W\vert x \sim$ uniform $\displaystyle {[0,\exp(-x^2/2)]}\;, \qquad X\vert w \sim$ uniform $\displaystyle {[-\sqrt{-2\log(w)},\sqrt{-2\log(w)}]} \;,$
yields $X \sim N(0,1)$ .

Pareto $(\alpha,\beta )$ , $\alpha > 0,\quad \beta > 0$

$f(x\vert \alpha, \beta ) = {{\beta {\alpha }^{\beta }}\over{x^{\beta +1}}}$ , $\alpha < x < \infty$ ,

$\displaystyle \frac{\alpha}{(1-U)^{1/\beta}} \sim$ Pareto $\displaystyle (\alpha,\beta)$

Student's $t_\nu$

$f(x\vert \nu ) = {{\Gamma \left( {{\nu +1}\over 2}\right) } \over{\Gamma \left... ...pi }}} {1\over{\left( 1 + \left({x^2\over \nu }\right)\right) ^{(\nu + 1)/2} }}$ , $-\infty < x < \infty$ , $\nu = 1,2,\dots$ .

$\displaystyle Y\sim\chi^2_\nu$ and $\displaystyle X\vert y \sim$ N $\displaystyle (0, \nu/y) \Rightarrow X \sim t_{\nu }.$

Also, if $X_1 \sim$ N

and $X_2 \sim \chi^2_\nu$ , then $X_1/\sqrt{X_2/\nu} \sim t_\nu$ .

Uniform

$f(x\vert a,b) = {1\over{b-a}}$ , $a \leq x \leq b$ ,

$\displaystyle (b-a)U + a \sim$ Uniform $\displaystyle (a,b).$

Weibull $(\gamma,\beta )$

$f(x\vert \gamma,\beta ) = {{\gamma }\over{\beta }}x^ {\gamma -1}e^{-x^{\gamma }/\beta }$ , $0 \leq x < \infty$ , $\gamma >0$ , $\beta > 0$ .

$\displaystyle X \sim {\rm Exponential}(\beta) \Rightarrow X^{1/\gamma } \sim {\rm Weibull}(\gamma,\beta ).$

References

Besag, J. and Green, P.J. (1993) Spatial statistics and Bayesian computation (with discussion). J. Roy. Statist. Soc. B 55, 25-38.
Damien, P., Wakefield, J. and Walker, S. (1999) Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables. J. Roy. Statist. Soc. B 61(2), 331-344.
Devroye, L. (1986) Non-Uniform Random Variate Generation. New York: Spring-Verlag.
Fishman, G.S. (1996) Monte Carlo. New York: Springer-Verlag
Gilks, W.R. and Wild, P. (1992) Adaptive rejection sampling for Gibbs sampling. Appl. Statist. 41, 337-348.
Hastings, W. (1970). Monte Carlo Sampling Methods Using Markov Chains and their Application. Biometrika 57 97-109.
Knuth, D. E. (1998). The Art of Computer Programming, Volume 2/ Seminumerical Algorithms, Third Edition. Addison Wesley, Reading.
Marsaglia, G. (1977) The squeeze method for generating gamma variables. Computers and Mathematics with Applications, 3, 321-325.
Marsaglia, G. and Zaman, A. (1991). A new class of random number generators. Annals Applied Probability 1, pp. 462-480.
Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E. (1953) Equations of state calculations by fast computing machines. J. Chem. Phys. 21, 1087-1092.
Ripley, B.D. (1987) Stochastic Simulation. New York: John Wiley
Robert, C. P. and Casella, G. (2004). Monte Carlo Statistical Methods, Second Edition. New York: Spring-Verlag.
Robert, C. P. and Casella, G. (2010). An Introduction to Monte Carlo Methods with R. New York: Spring-Verlag.
Rubinstein, R.Y. (1981) Simulation and the Monte Carlo Method. New York: John Wiley

Footnotes

... Robert ¹: Université Paris Dauphine, 75785 Paris cedex 16, France, and CREST, INSEE. Email: xian@ceremade.dauphine.fr
... Casella ²: Department of Statistics and Genetics Institute, University of Florida, Gainesville, FL 32611. Email: casella@ufl.edu

Generating Random Variables is owned by Christian P. Robert, George Casella, .

How to Cite This Entry

Christian P. Robert, George Casella, . "Generating Random Variables" (version 13). StatProb: The Encyclopedia Sponsored by Statistics and Probability Societies. Freely available at http://statprob.com/encyclopedia/GeneratingRandomVariables.html

Classification

AMS MSC:	60-08 (Probability theory and stochastic processes :: Computational methods )
	62-04 (Statistics :: Explicit machine computation and programs )
	62E99 (Statistics :: Distribution theory :: Miscellaneous)

Pending Errata and Addenda

None.

Username:
Password:

Register
I've forgotten my login details

Login

In Cooperation with:

Generating Random Variables

Generating Random Variables

Introduction

Generating Nonuniform Random Variables

Probability Integral Transform

Generating Discrete Variables

Accept-Reject Methods

Accept-Reject Algorithm:

Metropolis - Hastings Algorithm

The Gibbs Sampler

Noncentral chi squared ,

Normal

Footnotes

How to Cite This Entry

Classification

Pending Errata and Addenda

Discussion

Noncentral chi squared $(\lambda,p)$ , $\lambda \ge 0$

Normal $(\mu,{\sigma }^2)$