STAT0549-001: Statistical Computing and Lab

Note 2 - Methods for Generating Random Variables

Gyuhyeong Goh
Department of Statistics
Kyungpook National University

1
Introduction

In this chapter, our objective is to simulate random variables from probability

distributions using R.

2
 Sampling from a finite population

1 > #toss some coins

2 > sample(0:1, size = 10, replace = TRUE)
3 [1] 0 1 1 1 0 1 1 1 1 0
4 > #choose some lottery numbers
5 > sample(1:100, size = 6, replace = FALSE)
6 [1] 51 89 26 99 74 73
7 > #permutation of letters a-z
8 > sample(letters,size=5)
9 [1] "q" "a" "j" "x" "f"

3
 Random Generators of Common Probability Distributions

1 > n<-3
2 > rbeta(n,shape1=2,shape2=2)
3 [1] 0.9840114 0.4127489 0.2268932
4 > rbinom(n,size=5,prob=0.5)
5 [1] 1 1 3
6 > rchisq(n,df=2)
7 [1] 0.1987044 1.1532773 0.3874872
8 > rgamma(n,shape=3,rate=1)
9 [1] 2.7883450 0.6982289 2.3801886
0 > rnorm(n,mean=1,sd=2)
1 [1] 0.3269119 4.1151027 -0.1809909
2 > rpois(n,lambda=5)
3 [1] 3 8 3
4 > runif(n,min=0,max=2)
5 [1] 0.4537955 0.9852103 1.3223165

4
The Inverse Transform Method

Theorem (Probability Integral Transformation)

If X is a continuous random variable with cdf FX (x), then FX (X ) ∼ Uniform(0, 1).

This implies that if U ∼ Uniform(0, 1), then FX−1 (U) has the same
distribution as X .

5
Inverse Transform Method, Continuous Case

Let X ∼ FX (x). Then, a random sample of X can be generated as follows:

1. Generate u from Uniform(0, 1).
2. Compute x = FX−1 (u).

6
 Example
Inverse transform method, continuous case

Our goal is to generate a random sample with pdf fX (x) = 3x 2 , 0 < x < 1.
Rx
Note that FX (x) = 0 3t 2 dt = x 3 and FX−1 (u) = u 1/3 .
Then, the R code can be written as follows:

1 > n <- 1000

2 > u <- runif(n)
3 > x <- u^(1/3)

7
Example (cont.)
Inverse transform method, continuous case

f(x) = 3x2

2.5
2.0
1.5
Density

1.0
0.5
0.0

0.0 0.2 0.4 0.6 0.8 1.0

8
Inverse Transform Method, Discrete Case

Let X be a discrete random variable on the support {x1 , x2 , . . .} such that

x1 < x2 < · · · .

The inverse transform method can be implemented as follows:

1. Generate u from Uniform(0, 1).
2. Select x = xi such that F (xi−1 ) < u ≤ F (xi ), where F (x0 ) = 0.

9
 Example: Bernoulli distribution
Inverse Transform Method, Discrete Case

Let X ∼ Bernoulli(p = 0.4). Then, x1 = 0, x2 = 1, F (x1 ) = 0.6, and

F (x2 ) = 1.
Recall that
1. Generate u from Uniform(0, 1).
2. Select x = xi such that F (xi−1 ) < u ≤ F (xi ).
Hence, we generate u from Uniform(0, 1) and then set x = 0 if u ≤ 0.6 and
x = 1 if 0.6 < u ≤ 1.

1 > n <- 1000

2 > p <- 0.4
3 > u <- runif(n)
4 > x <- [Link](u > 0.6) #(u > 0.6) is a logical vector

10
 Example: Geometric distribution
Inverse Transform Method, Discrete Case

Let X ∼ Geometric(p = 0.25). Then, the cdf is F (x) = 1 − q x+1 ,

x = 0, 1, 2, . . . , where q = 1 − p.
Recall that
1. Generate u from Uniform(0, 1).
2. Select x = xi such that F (xi−1 ) < u ≤ F (xi ).
Hence, we generate u from Uniform(0, 1) and then set x = xi such that

log(1 − u)
1 − q xi < u ≤ 1 − q xi +1 ⇔ xi < ≤ xi + 1.
log(q)

1 > n <- 1000

2 > p <- 0.25
3 > u <- runif(n)
4 > x <- ceiling(log(1-u)/log(1-p))-1

11
The Acceptance-Rejection (AR) Method

Let X be a random variable with pdf (or pmf) f (x).

Our goal is to generate a sample of X , but it is impossible to draw a sample

from f (x) directly.

Let Y be a random variable with pdf (or pmf) g (y ) such that f (t)/g (t) ≤ c
for some constant c.
The AR method proceeds as follows:
1. Generate y from the distribution with g .
2. Generate u from from Uniform(0, 1).
3. If u < f (y )/{cg (y )}, then accept y . i.e., x = y ;
otherwise reject y go to Step 1.

Note that g and c are determined by us.

12
 Example: Beta(2,2)
AR method
Let X ∼ Beta(2, 2). Then, the pdf is f (x) = 6x(1 − x), 0 < x < 1.

Let g (x) be the Uniform(0, 1) density, i.e., g (x) = 1, 0 < x < 1.

Note that
f (x) 6x(1 − x)
= ≤ 6 ≡ c.
g (x) 1
f (y )
Hence, we accept if u < cg (y ) = y (1 − y ).
1 n <- 1000
2 x <- rep(0,n) #storage for generated x-values
3 k <- 0 #count of accepted samples
4 while (k < n){
5 u <- runif(1)
6 y <- runif(1) #random variate from g
7 if (y * (1-y) > u) { #we accept y if the condition is true
8 k <- k + 1
9 x[k] <- y
0 }
1 }
13
Example: Beta(2,2) (cont.)
AR method

1.0
●
●●●
●●●
●
●
●
●●
●●●●●
●
●●
●
●
●
●●
●●
●
●●
●
●●
●
●
●●
●
●
●●
●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●●
●
●●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
0.8
●
●
●●
●
●
●●
●
●
●●
●
●
●
●●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●●
●
●
●●
●
●
●●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●●
●
●●
●●
●●
●
●
●●
●
●●
●
●
●●
●
●
●
●●
●
●
●
●●
●
●●
●●
●
●●
●
●●
●●
●
●
●
●●
●
●●
●
●
●
●●
●
●●
●
●
●
●●
●●
●
●
●
●
●●
●
●●
●
●
0.6

●
●●
●
●
●
●●
●●
●
●●
●
●●
●
●●
●
●●
●
●
●
●●
●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●
Sample

●●
●
●●
●
●●
●
●●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●●
●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●
●●
●
●
●●
●
●●
●
●
●●
●
●
●●
●
●●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●●
●
●●
●
●
●●
●
●●
●
●●
0.4

●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●
●●
●●
●
●●
●●
●
●
●●
●
●●
●
●
●●
●
●
●
●●
●
●
●●
●
●●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●●
●
●
●
●●
●●
●
●●
●
●●
●
●●
●
●●
●●
●
●●
●●
●●
●
●●
●
●
●
●●
●●
●
●
●
●●
●
●
●
●
●
●●
0.2

●●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●●
●●
●
●
●●
●
●●
●
●
●
●
●●
●
●●
●
●
●
●
●
●●
●
●●
●
●
●
●●
●
●●
●
●●
●
●
●●
●
●
●
●●
●●
●
●
●
●●
●
●●
●
●
●●
●
●
●
●
●
●
●●
●●
●
●●
●●
●
●●●
●●●
●●
●●
●
0.0

0.0 0.2 0.4 0.6 0.8 1.0

Beta(2, 2)

14
Remark on AR

Let T be the number of iterations.

It can be shown that P(Accept) = 1/c.

Hence, the average number of the accepted samples is T /c.

This implies that choosing c as small as possible improves the efficiency of

AR method.

15
Exercise 1

Suppose that we are interested in sampling from U(0, k).

Using the inverse transform method, generate 1, 000 samples from

U(0, k = 10).

16
Exercise 2

Suppose that we are interested in sampling from Poisson(λ).

Using the inverse transform method, generate 1, 000 samples from

Poisson(λ = 2).

17
Exercise 3

Suppose that we are interested in sampling from a half-normal distribution

with pdf √
2
f (x) = √ exp(−x 2 /2), x > 0.
π
Using the AR method, generate 1, 000 samples from the half-normal
distribution.

18