Chapter 2 Simulating random variables and vectors

Discrete distributions in R

Distributions Random numbers generation
Binomial, \({\rm B}(n,p)\) rbinom(size,prob)
Geometric, \({\mathcal G}(p)\) rgeom(prob)
Negative Binomial, \({\rm NB}(k,p)\) rnbinom(size,prob)
Hypergeometric, \({\rm H}(N_1,N_2,k)\) rhyper(m,n,k)
Poisson, \({\mathcal P}(\lambda)\) rpois(lambda)

Continuous distributions in R

Distributions Random number generation
Uniform, \({\rm U}(a,b)\) runif(min=0,max=1)
Exponential, \({\rm Exp}(\lambda)\) rexp(rate=1)
Normal, \({\rm N}(\mu,\sigma)\) rnorm(mean=0,sd=1)
Gamma, \({\rm Gamma}(k,\lambda)\) rgamma(shape,rate=1)
Beta, \({\rm Beta}(\alpha,\beta)\) rbeta(shape1,shape2)
Chi-square, \(\chi^2_n\) rchisq(df)
Student’s \(t\), \(t_n\) rt(df)
Fisher’s \(F\), \(F_{n_1,n_2}\) rf(df1,df2)

2.1 Inverse transform

Consider rv \(X\) with cdf \(F(x)=P(X\leq x)\) and having as quantile function the generalized inverse function of \(F\), \[F^{-1}(u) = \inf \{x : \, F(x) \geq u \}.\]

If \(U\sim {\rm U}(0,1)\) we have \(P(U\leq u)=u\) for any \(u\in [0,1]\) and given rv \(X=F^{-1}(U)\), its cdf is

that is, \(X\sim F\).

Inverse transform for discrete random variables

Consider \(X\) discrete rv that can assume the \(k\) values \(x_1,x_2, \ldots, x_k\) with probabilities \(p_1\geq p_2\geq \ldots\geq p_k\). The cdf of \(X\) is \(F(x)=\sum_{i:x_i\leq x} p_i\). The inverse transform algorithm to obtain values from \(X\) at random is

  1. Generate \(U\) random number in \((0,1)\).
  2. Set \(X=x_i\) if \(F(x_{i-1}) < U \leq F(x_i)\) and stop.

The values in the support of \(X\) have been ordered in terms of their probabilities in order to gain efficiency.

We are simulated from rv \(X\) with probability mass function \[P(X=0)= 0.05 \quad P(X=1)= 0.1 \quad P(X=2)=0.45 \quad P(X=3)=0.4. \] Which algorithm is more efficient?

  1. Generate \(U\) random number in \((0,1)\).
  2. If \(U<0.05\) set \(X=0\) and stop.
  3. If \(U<0.15\) set \(X=1\) and stop.
  4. If \(U<0.6\) set \(X=2\) and stop.
  5. Set \(X=3\) and stop.
  1. Generate \(U\) random number in \((0,1)\).
  2. If \(U<0.45\) set \(X=2\) and stop.
  3. If \(U<0.85\) set \(X=3\) and stop.
  4. If \(U<0.95\) set \(X=1\) and stop.
  5. Set \(X=0\) and stop.
sim.disc1 <- function(MC) {
  sim.disc <- rep(3,MC)
  sim.u <- runif(MC)
  for(i in 1:MC){
    if(sim.u[i]<0.05) sim.disc[i] <- 0
    else if(sim.u[i]<0.15) sim.disc[i] <- 1
    else if(sim.u[i]<0.6) sim.disc[i] <- 2
  }
  return(sim.disc)
} 
sim.disc2 <- function(MC) {
  sim.disc <- rep(0,MC)
  sim.u <- runif(MC)
  for(i in 1:MC){
    if(sim.u[i]<0.45) sim.disc[i] <- 2
    else if(sim.u[i]<0.85) sim.disc[i] <- 3
    else if(sim.u[i]<0.95) sim.disc[i] <- 1
  }
  return(sim.disc)
} 
sim.disc3 <- function(MC) {
  sim.disc <- rep(3,MC)
  sim.u <- runif(MC)
  zero <- which(sim.u<0.05) ; sim.disc[zero] <- 0
  one <- which(sim.u>0.05 & sim.u<0.15) ; sim.disc[one] <- 1
  two <- which(sim.u>0.15 & sim.u<0.6) ; sim.disc[two] <- 2
  return(sim.disc)
} 
require(microbenchmark)
## Loading required package: microbenchmark
## Warning: package 'microbenchmark' was built under R version 3.6.1
microbenchmark(sim.disc1(10),sim.disc2(10),sim.disc3(10))
## Unit: microseconds
##           expr   min    lq     mean median    uq      max neval
##  sim.disc1(10) 3.138 3.706 81.28680  3.992 4.562 7660.430   100
##  sim.disc2(10) 2.850 3.422 76.91332  3.706 4.278 7259.576   100
##  sim.disc3(10) 4.846 5.418 83.22548  5.989 7.128 7651.022   100
microbenchmark(sim.disc1(1000),sim.disc2(1000),sim.disc3(1000))
## Unit: microseconds
##             expr     min      lq      mean  median      uq     max neval
##  sim.disc1(1000) 136.280 138.275 147.73858 138.848 140.843 428.226   100
##  sim.disc2(1000) 112.616 114.043 120.12362 114.898 116.038 264.576   100
##  sim.disc3(1000)  60.158  61.014  63.08006  61.584  62.296  96.366   100

We can use the chi-square GoF test to check that any of the algorithms that we have written actually samples from the desired distribution.

MC <- 1000 ; set.seed(1)
sim.disc <- sim.disc3(MC)
chisq.test(x=table(sim.disc),p=c(0.05,0.1,0.45,0.4))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(sim.disc)
## X-squared = 1.5389, df = 3, p-value = 0.6733

Inverse transform for discrete uniform random variable

Consider \(X\) discrete rv with \(P(X=r)=1/k\) for \(x\in\{1,2,\ldots,k\}\).

\[F_X(x)=\left\{\begin{array}{ll}0&\textrm{ if }x<1\\r/k&\textrm{ if }r\leq x<r+1,\,i\in\{1,2,\ldots,k-1\}\\1&\textrm{ if }x\geq k\end{array}\right..\]

\[F^{-1}_X(u)=r \textrm{ if }\frac{r-1}{k}<u\leq\frac{r}{k},\quad\textrm{equiv.}\quad F^{-1}_X(u)=\lfloor{uk}\rfloor+1\]

  1. Generate \(U\) random number in \((0,1)\).
  2. Return \(X=\lfloor{U*k}\rfloor+1\). (or equiv. \(X=\lceil{U*k}\rceil\))
MC <- 1000 ; set.seed(2)
k <- 5
sim.dunif <- ceiling(runif(MC)*k)
chisq.test(x=table(sim.dunif),p=rep(1,k)/k)
## 
##  Chi-squared test for given probabilities
## 
## data:  table(sim.dunif)
## X-squared = 5.71, df = 4, p-value = 0.2219

Inverse transform for geometric random variable

Consider \(X\sim{\mathcal G}(p)\) with \(P(X=r)=p(1-p)^r\) for \(r\in\{0,1,2,\ldots\}.\)

\[F_X(x)=\left\{\begin{array}{ll}0&\textrm{ if }x<0\\1-(1-p)^{r+1}&\textrm{ if }r\leq x<r+1,\,r\in\{0,1,\ldots,k-1\}\end{array}\right..\]

\[F^{-1}_X(u)=\lfloor\ln(1-u)/\ln(1-p)\rfloor\]

  1. Generate \(U\) random number in \((0,1)\).
  2. Set \(X=\lfloor\ln(U)/\ln(1-p)\rfloor\) and stop.
MC <- 1000 ; set.seed(2)
p <- 0.1
sim.geom1 <- floor(log(runif(MC))/log(1-p))
sim.geom2 <- rgeom(MC,prob=p)
plot(sort(sim.geom1),sort(sim.geom2))
abline(0,1)

Inverse transform for continuous random variables

Consider \(X\) continuous rv with cdf \(F\) and quantile function \(F^{-1}\). The inverse transform algorithm to obtain values from \(X\) at random is

  1. Generate \(U\) random number in \((0,1)\).
  2. Set \(X=F^{-1}(U)\) and stop.

Example (Exponential)

Let \(X\sim{\rm Exp}(\lambda)\), its cdf is \(F(x)=1- e^{-\lambda x}\) for \(x>0\), while its quantile function is \[\begin{multline*} u\in(0,1),\,\, F^{-1}(u)=x \Leftrightarrow u=F(x) \Leftrightarrow u=1- e^{-\lambda x} \Leftrightarrow 1-u=e^{-\lambda x} \\ \Leftrightarrow -\mbox{ln}(1-u)=\lambda x \Leftrightarrow x=-\mbox{ln}(1-u)/\lambda\,. \end{multline*}\]

Inverse transform for exponential random variable

  1. Generate \(U\) random number in \((0,1)\).
  2. Set \(X=-\mbox{ln}(1-U)/\lambda\) (or \(X=-\mbox{ln}(U)/\lambda\)) and stop.
MC <- 1000 ; set.seed(1)
sim.exp <- -log(runif(MC))
ks.test(sim.exp,"pexp",rate=1)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.exp
## D = 0.024366, p-value = 0.5928
## alternative hypothesis: two-sided

2.2 Aceptance-rejection

Our goal is to simulate rvs with dmf (pmf) \(f\). We need an instrumental dmf (pmf) \(g\) from which we can simulate. Some values simulated from \(g\) will be rejected. Conditions on \(g\):

  • If \(f(x)>0\), then \(g(x)>0\).
  • There exists \(M>0\) such that \(f(x)/g(x) \leq M\) for every \(x\).
  1. Generate \(Y\) with density (mass) \(g\).
  2. Generate \(U\) random number in \((0,1)\).
  3. If \(U<\frac{f(Y)}{Mg(Y)}\) set \(X=Y\) and stop. Otherwise, return to Step 1. \[\begin{multline*} P\left(Y\leq x | U< \frac{f(Y)}{Mg(Y)}\right) = \dfrac{P\left(Y\leq x , U< \frac{f(Y)}{Mg(Y)}\right)}{P\left(U< \frac{f(Y)}{Mg(Y)}\right)}\\ =\dfrac{\int_{-\infty}^x F_U\left(\frac{f(y)}{Mg(y)}\right)\,g(y)dy}{\int_{-\infty}^{\infty}F_U\left(\frac{f(y)}{Mg(y)}\right)\, g(y)dy} =\dfrac{\int_{-\infty}^x\frac{f(y)}{Mg(y)} g(y)dy}{\int_{-\infty}^{\infty}\frac{f(y)}{Mg(y)} g(y)dy}=P(X\leq x).\\ \end{multline*}\]

Rejection for \(|Z|\) with \(Z\sim{\rm N}(0,1)\)

The dmf of \(|Z|\) with \(Z\sim{\rm N}(0,1)\) is \(f(x)=\frac{2}{\sqrt{2\pi}}e^{-x^2/2}\) for \(x>0\), while the dmf of an Exponential random variable with \(\lambda=1\) is \(g(x)=e^{-x}\) for \(x>0\). \[\frac{f(x)}{g(x)}=\frac{\frac{2}{\sqrt{2\pi}}e^{-x^2/2}}{e^{-x}}=\sqrt{2/\pi}e^{x-x^2/2}\leq\sqrt{2e/\pi}.\] Set \(M=\sqrt{2e/\pi}\) and \[\frac{f(x)}{Mg(x)}=e^{x-x^2/2-1/2}=e^{-(x-1)^2/2}\]

  1. Generate \(Y\) with density (mass) \(g\).
  2. Generate \(U\) random number in \((0,1)\).
  3. If \(U<\exp\{-(Y-1)^2/2\}\) set \(X=Y\) and stop. Otherwise, return to Step 1.
s.absnorm <- function(n) {
  y <- vector(length=n)
  for(i in 1:n){
    y[i] <- -log(runif(1))
    while(runif(1)>exp(-(y[i]-1)^2/2)) y[i] <- -log(runif(1))
  }
  return(y)
}
MC <- 1000 ; set.seed(1)
hist(s.absnorm(MC),probability=T)

Rejection for truncated models

A Zero Trucated Poisson rv \(X\) with parameter \(\lambda\) has pmf \[p(r)=P(X=r)=\frac{e^{-\lambda}}{1-e^{-\lambda}}\frac{\lambda^r}{r!}\,,\quad r\in\{1,2,\ldots\}\] Set \(q\) equal to the pmf of a \({\mathcal P}(\lambda)\) rv and \(M=(1-e^{-\lambda})^{-1}\), then \(p(r)/(Mq(r))=1\) (and accepted) if \(r\geq 1\), and \(p(0)/(Mq(0))=0\) (and rejected).

  1. Generate \(Y\) a Poisson rv.
  2. If \(Y\geq 1\) set \(X=Y\) and stop. Otherwise, return to Step 1.

2.3 Composition approach

If we must simulate a random variable \(X\) with a mixture distribution either finite (or countable) or continuous \[F(x)=\sum_{i\in I}p_iF_i\quad\textrm{or}\quad F(x)=\int_A\omega(a)F_a(x)\,{\rm d}a\] and we are able to simulate from \(F_i\) (\(F_a\)) and a rv with pmf \((p_i)_{i\in I}\) or dmf \(\omega\), then:

  1. Generate \(Y\) with pmf \((p_i)_{i\in I}\) or dmf \(\omega\).
  2. Generate \(X\) with cdf \(F_Y\) and stop.

Generating a Laplace random variable

The Laplace (double exponential) distribution model with location parameter \(m\in{\mathbb R}\) and rate parameter \(\lambda>0\) has density mass function

\[f(x)=\frac{\lambda}{2}e^{-\lambda|x-m|}\,.\]

It is thus the mixture of an Exponential rv with parameter \(\lambda\) and minus an Exponential rv with rate parameter \(\lambda\) with equal weights and shifted by \(m\) units.

  1. Generate \(U_1,U_2\) indep. random numbers in \((0,1)\).
  2. Return \(X={\rm sign}(U_1-0.5)\ln(U_2/\lambda)+m\).
MC <- 1000 ; set.seed(2)
m <- 2 ; lmbd <- 1
s.lap <- sign(runif(MC)-0.5)*log(runif(MC)/lmbd)+m
hist(s.lap,prob=T)

Normal rv (rejection + composition)

MC <- 1000 ; set.seed(1)
mu <- 2 ; sigma <- 0.5
s.norm <- sigma*(sign(runif(MC)-0.5)*s.absnorm(MC))+mu
shapiro.test(s.norm)$p.value
## [1] 0.6076152
qqnorm(s.norm) ; qqline(s.norm)

2.4 Multivariate distributions

Generate random vector with known conditional cdfs

\[F_{X_1,X_2,X_3}(x_1,x_2,x_3)=F_{X_1}(x_1)F_{X_2|X_1}(x_2|x_1)F_{X_3|X_1,X_2}(x_3|x_1,x_2)\]

If we can simulate from the conditional distributions of a random vector, we can do as follows

  1. Generate \(X_1\) with cdf \(F_{X_1}\).
  2. Generate \(X_2\) with cdf \(F_{X_2|X_1}\).
  3. Generate \(X_3\) with cdf \(F_{X_3|X_1,X_2}\).
  4. Set \({\bf X}=(X_1,X_2,X_3)^t\) and stop.

2.4.1 Multinomial distribution

We use the property that the univariate marginal and conditional distributions of a Multinomial random vector are Binomial in order to simulate the Multinomial model.

\[{\bf X}=(X_1,X_2,X_3,X_4)^t\sim{\rm M}(n,(p_1,p_2,p_3,1-p_1-p_2-p_3))\]

  1. Generate \(X_1\sim{\rm B}(n,p_1)\).
  2. Generate \(X_2\sim{\rm B}(n-X_1,p2/(1-p_1))\).
  3. Generate \(X_3\sim{\rm B}(n-X_1,p3/(1-p_1-p_2))\).
  4. Set \(X_4=n-X_1-X_2-X_3\).
  5. Set \({\bf X}=(X_1,X_2,X_3,X_4)^t\) and stop.

Alternatively use rmultinom(n, size, prob).

Random permutations

A permutation \(\pi\) of \(S=\{1,2,\ldots,n\}\) is any of the \(n!\) sortings of \(S\). That is \(\pi(i)\in S\) for \(i\in S\) and \(\pi(i)\neq\pi(j)\) for \(i,j\in S\) with \(i\neq j\).

We can use the condtional cdfs argument to generate a permutation at random from \(S=\{1,2,3,4,5\}\) as follows.

  1. Set \(\pi(1)\) random number from \(S\).
  2. Set \(\pi(2)\) random number from \(S\setminus\{\pi(1)\}\)
  3. Set \(\pi(3)\) random number from \(S\setminus\{\pi(1),\pi(2)\}\)
  4. Set \(\pi(4)\) random number from \(S\setminus\{\pi(1),\pi(2),\pi(3)\}\)
  5. Set \(\pi(5)\in S\setminus\{\pi(1),\pi(2),\pi(3),\pi(4)\}\) and stop.

If \(A\) and \(B\) are two sets, \(A\setminus B\) is the set formed by the elements in \(A\) that are not in \(B\).

k <- 10 ; perm <- (1:k)
for(i in 1:(k-1)){
  interchange <- ceiling(runif(1)*(k-i+1))
  aux <- perm[i]
  perm[i] <- perm[interchange+(i-1)]
  perm[interchange+(i-1)] <- aux 
}
perm
##  [1]  2 10  4  6  3  5  7  8  1  9

Alternatively

rank(runif(10)) ; sample(10,x=1:10,replace=F)
##  [1]  5  1  3  7  6 10  9  2  8  4
##  [1]  8  5  7  6  3  4  9  2  1 10

2.5 Multivariate normal distribution

Bivariate normal generation (Box-Muller’s polar method)

The transformation algorithm below can be used to simulate a bivariate standard normal rv.

  1. Generate \(U_1,U_2\) independent random numbers in \((0,1)\).
  2. Set \((X_1,X_2)=\sqrt{-2\log U_1}(\cos(2\pi U_2),\sin(2\pi U_2))\) and stop.
  • If \((X_1,X_2)\) represents a random point in the plane with \(X_1,X_2\) independent standard normal rvs.
  • The polar coordinates of \((X_1,X_2)\) are \(\sqrt{X_1^2+X_2^2}\) and \(\arctan2(X_2,X_1)\), where \(\arctan2(\cdot,\cdot)\) is the \(2\)-argument arctangent.
  • The distance to the origin \(\sqrt{X_1^2+X_2^2}\) is the square root of a chi-square rv with 2 degrees of freedom (equivalently the square root of an exponential with \(\lambda=1/2\)).
  • The angle \(\arctan2(X_2,X_1)\) is a random angle in \((0,2\pi)\).

If we transform a bivariate standard normal rv \((Z_1,Z_2)\sim{\rm N}_2({\bf 0},{\bf I}_2)\) as \(X_1=\mu_1 +\sigma_1 Z_1\) and \(X_2=\mu_2 +\sigma_2 (Z_1\rho + Z_2 \sqrt{1-\rho^2})\) then

  • \({\mathbb E}[X_1]={\mathbb E}[\mu_1 +\sigma_1 Z_1]=\mu_1\);
  • \({\rm Var}[X_1]={\rm Var}[\mu_1 +\sigma_1 Z_1]=\sigma_1^2\);
  • \({\mathbb E}[X_2]={\mathbb E}[\mu_2 +\sigma_2 (Z_1\rho + Z_2 \sqrt{1-\rho^2}]=\mu_2\);
  • \({\rm Var}[X_2]={\rm Var}[\mu_2 +\sigma_2 (Z_1\rho + Z_2 \sqrt{1-\rho^2}]=\sigma_2^2(\rho^2+1-\rho^2)=\sigma_2^2\);
  • \(Cov(X_1,X_2)={\mathbb E}[X_1X_2]-{\mathbb E}[X_1]{\mathbb E}[X_2]=\mu_1\mu_2+\rho\sigma_1\sigma_2-\mu_1\mu_2=\rho\sigma_1\sigma_2\);
  • The joint distribution of \((X_1,X_2)\) is bivariate normal.
  1. Generate \(Z_1\) and \(Z_2\) independent standard normal rvs.
  2. Return \(X_1=\mu_1 +\sigma_1 Z_1\) and \(X_2=\mu_2 +\sigma_2 (Z_1\rho + Z_2 \sqrt{1-\rho^2})\).
mu1 <- 1 ; mu2 <- 3 ; sigma1 <- 1 ; sigma2 <- 2; rho <-0.8 
MC <- 1000 ; set.seed(1)
s.bnorm <- matrix(rnorm(2*MC),nrow=2)
m.trans <- matrix(c(sigma1,0,sigma2*rho,sigma2*sqrt(1-rho^2)),
                ncol=2,byrow=T)
s.norm <- t(m.trans%*%s.bnorm+c(mu1,mu2))
plot(s.norm)

Multivariate normal

Consider a \(d\)-variate normal distribution with covariance matrix \(\Sigma\) and mean vector \(\mu\).

  1. Determine matrix \(C\) such that \(C C^t=\Sigma\).
  2. Generate \(Z_1,\ldots,Z_d\) independent standard normal rvs.
  3. Set \({\bf X}=C (Z_1,\ldots,Z_d)^t+\mu\) and stop.

The Choleski decomposition of the covariance matrix \(\Sigma\) is a triangular matrix \(L\) such that \(\Sigma=L L^t\), where \(L\) is computed in R as chol(\(\Sigma\)) and the ouput is an upper triangular matrix. In case we need a lower triangular matrix, we simply take the Choleski decomposition of \(\Sigma\) and traspose it.

The matrix that we computed in the bivariate example equals the one we obtain with the Choleski decomposition.

m.trans
##      [,1] [,2]
## [1,]  1.0  0.0
## [2,]  1.6  1.2
Sigma<-matrix(c(sigma1^2,rho*sigma1*sigma2,
                rho*sigma1*sigma2,sigma2^2),ncol=2,byrow=T)
t(chol(Sigma)) # Use C=t(chol(Sigma))
##      [,1] [,2]
## [1,]  1.0  0.0
## [2,]  1.6  1.2

Alternatively use rmvnorm(n,mean,sigma) from package mvtnorm.

require(mvtnorm)
mu1 <- 1 ; mu2 <- 3 ; sigma1 <- 1 ; sigma2 <- 2; rho <-0.8 
Sigma<-matrix(c(sigma1^2,rho*sigma1*sigma2,
                rho*sigma1*sigma2,sigma2^2),ncol=2,byrow=T)
MC <- 1000 ; set.seed(1)
plot(rmvnorm(MC,mean=c(mu1,mu2),sigma=Sigma))