GENERATING RANDOM VARIATES
REJECTION SAMPLING
Remember that there are some techniques to generate random sample (numbers) from a distribution which are;
– Inverse transform method
– Acceptance-rejection method
– Convolution method
– Composition method
We cannot use inverse transform method distributions which do not have closed from of CDF like normal, gamma, beta, etc. In such a cases we can use Acceptance-rejection method.
Acceptance-rejection method
For instance, if one wants to generate 100 random numbers that belong to a normal distribution in R, it is as simple as executing:
rnorm(100)However, how does this process actually work under the algorithm? How can an algorithm know whether a random number belongs to a particular distribution or not?
The answer is through rejection sampling.
Acceptance-Rejection sampling is a way to simulate random samples from an unknown or difficult to sample from distribution by using random samples from a similar probability distribution. A random subset of the generated samples are rejected; the rest are accepted. The goal is for the accepted samples to be distributed as if they were from the unknown probability distribution.
How Rejection Sampling Works
Assume that we have a random variable x, which has density function f(x). We want to generate data from this density function f(x).
To do this find an alternative probability distribution G, with density function g(x), from which we already have an efficient algorithm for generating from, but also such that the function g(x) is “close” to f(x). Then, we assume that the ratio f(x)/g(x) is bounded by a constant c > 0.
We can find a maximum x point that maximize the f(x), by taking derivative and equalize it to 0.
Continuous Case
Example 1)
Create an algorithm (function) to generate data from the density function below.
The function is:
f(x)=100x³(1−x)², 0≤x≤1
g(x)=1, 0≤x≤1.
After creating sample from that beta distribution, estimate P(0 ≤ X ≤ 0.8) with k = 1000 repetitions (n=1000).
Let applied the procedure,
1. We choose density g(x)=1, 0 ≤ x ≤ 1
2. Find constant c;
c= max(f(x)/g(x))
That is why c=3.5 is chosen.
3. Generate a random variable Y distributed as G(Y). g(x)=1 which is standard uniform.
4. Generate a uniform random number U. u= runif(1)
5. If U ≤ 100x³(1-x)²/(1*3.5) then accept set Y=X, else reject and go to step 3.
beta.rejection <- function(f, c, g, n) {
set.seed(2022)
k <- 0 # counter for accepted
i <- 0 # number of iteration
s <- c() #newly created random sample
while (k < n) {
x <- runif(1)
i <- i+1
u <- runif(1)
if (u <= f(x) / (c*g(x))) {
k <- k + 1
s[k] = x
}
}
print(summary(s))
return(s)
}
f <- function(x) 100*(x^3)*(1-x)^2
g <- function(x) 1
c <- 3.5 # an upper bound for f(x)/g(x)
prob.result <- beta.rejection(f, c, g, n=1000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09818 0.44099 0.58905 0.57391 0.70535 0.98020print(length(prob.result[prob.result <= 0.8])/1000)## [1] 0.9hist(prob.result,freq = FALSE)
points(seq(0,1,0.01), dbeta(seq(0,1,0.01),4,3),type = "l")
Example 2)
Write a function which is generated a random variable from
f(x)= x²/18, -3 ≤x≤ 3
using acceptance — rejection method and generate n=10000 random samples.
Assume that g(x) = 1, 0 ≤ x ≤ 1
Let applied the procedure,
1. We choose density g(x)=1, 0 ≤ x ≤ 1
2. Find constant c;
When x=0, f(x) is also 0. In this way we cannot find maximum point. However, by putting x=3 which is maximum value that x can take, we can find maximum value of f(x). Thus, when x=3, f(x)/g(x)=1/2
That is why c=0.5 is chosen.
3. Generate a random variable Y distributed as G(Y). g(x)=1 which is standard uniform
4. Generate a uniform random number U. u= runif(1)
5. If U ≤ x²/(18*1*0.5) then accept set Y=X, else reject and go to step 3
set.seed(2022) # for reproducibility
n = 10^4 # sample size
i <- 0 # number of iteration
k <- 0 # counter for accepted
y <- numeric(n) # random sample
while (k < n) {
u = runif(1)
i <- i+1
x = runif(1, -3,3)
if (u <= x^2/(18*1*0.5)) {
k <- k + 1
y[k] <- x
}
}
i;k;length(y)## [1] 30244
## [1] 10000
## [1] 10000hist(y, prob = TRUE, main ="f(x)=x^2/18")
a = seq(-3,3,0.01)
lines(a,a^2/18,col="red",lwd = 2.5)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.99960 -2.38161 0.60966 0.01376 2.37344 2.99960
We can write the code also without while loop;
set.seed(1234)
n = 100000
x = runif(n, -3, 3)
u = runif(n)
# acceptance condition
s = x[u <= x^2/18]
length(s)## [1] 16659hist(s, prob = TRUE, main ="f(x)=x^2/18")
a = seq(-3,3,0.01)
lines(a,a^2/18,col="red",lwd = 2.5)
plot(x,u, pch=".", col="maroon")
points(s,u[u <= x^2/18], pch=".", col="green3")
Example 3)
Write an algorithm to generate n=10000 random variables from Gamma(α = 4, β = 3) by using the relation between Exponential and Gamma distributions. Graph the histogram of the sample with the theoretical Gamma(α = 4, β = 3) density.
Reminding!
x ∼ Exp(λ) = Gamma(1, λ)
set.seed(361)
gamma <- function(n,t,lambda) {
U = matrix(runif(n*t),t, n)
logU = -log(U) / lambda # exponential rv by using inverse transform method
X = colSums(logU) # column sums of matrix logU
hist(X, prob=TRUE, main ="Gamma Distribution (3,2)")
y = seq(0,6,0.01)
lines(y,dgamma(y,3,2),col="red") # theoretical curve for gamma dist.
}
n = 10
t = 3
lambda = 2
gamma(n,t,lambda)
Discrete Case
Now, suppose we want to generate an X that is a discrete random variable with probability mass function(PMF) p(k)=P(X=k). Further suppose that we can easily generate a discrete random variable Y with PMF q(k)=P(Y=k) such that p(k) /q(k) ≤ c.
- Step 1: Generate a random variable Y which is distributed as q(k).
- Step 2: Generate U (independent from Y).
- Step 3: If u ≤ p(k)/(q(k)*c), then set X = Y; otherwise go back to step 1.
Example 4)
Simulate a discrete random variable X that has PMF given by;
P(X = 1) = 0.12
P(X = 2) = 0.23
P(X = 3) = 0.18
P(X = 4) = 0.37
P(X = 5) = 0.07
P(X = 6) = 0.03
Lets apply the procedure above.
Step 1: Generate a random variable Y which is distributed as q(k).
Since the domain of X is 1, 2, 3, 4, 5, 6 we can assume that random variable Y has discrete uniform distribution between (1,6). The PMF of discrete uniform distribution is
q(k)=P(Y=k)= 1/6, k=1,2,3,4,5,6
Step 2: Generate U (independent from Y). = runif(1)
Step 3: If u ≤ p(k)/(q(k)*c), then set X = Y; otherwise go back to step 1.
Firstly find minimum value of c value by finding maximum value of pmf p(k).
p(k)/q(k) ≤ c
Since the biggest pmf for X is 0.37, we can choose minimum c by the equation:
n <- 10000 # sample size
probs <- c(0.12, 0.23, 0.18, 0.37, 0.07, 0.03) # probabilities for desired density
c <- 2.219 #finded "c" value above
cq <- c * (1/6)
cq## [1] 0.3698333# for accepted variates
accepted_x <- numeric(n)
accepted_y <- numeric(n)
# for rejected variates
rejected_x <- numeric(n)
rejected_y <- numeric(n)
i <- 1 # initialization of index for accepted r.v’s
j <- 1 # initialization of index for rejected r.v’s
set.seed(2022)
while(i <= n){
k <- sample(1:6, size = 1) # random number for q(k)
u <- runif(1) # random number for comparison
cq <- 0.3698333
if (u <= probs[k]/cq){
accepted_x[i] <- k
accepted_y[i] <- u * cq
i <- i + 1
}else{
rejected_x[j] <- k
rejected_y[j] <- u * cq
j <- j + 1
}
}
plot(accepted_x, accepted_y, type="p", col="blue",ylim=c(0,1))
points(rejected_x, rejected_y,col="red")
table(accepted_x)/n## accepted_x
## 1 2 3 4 5 6
## 0.1219 0.2254 0.1859 0.3667 0.0716 0.0285x<- c(1,2,3,4,5,6)
theoretical_probabilities <- c(0.12, 0.23, 0.18, 0.37, 0.07, 0.03)
estimated_probabilities <- c(0.1219, 0.2254, 0.1859, 0.3667, 0.0716, 0.0285)
d<- data.frame(x,theoretical_probabilities,estimated_probabilities)d## x theoretical_probabilities estimated_probabilities
## 1 1 0.12 0.1219
## 2 2 0.23 0.2254
## 3 3 0.18 0.1859
## 4 4 0.37 0.3667
## 5 5 0.07 0.0716
## 6 6 0.03 0.0285