ST303 Stochastic Simumlation

Class 5

Philipp Sterzinger

p.sterzinger@lse.ac.uk

Office hours: Thursday, COL 8.05

13 November 2025

Recap

Mixing

Assume we have a random variable \(X\) with pdf \(f_X\) such that

\[ f_X(x) = \sum_{i = 1}^N \alpha_i f_i(x) \] where \(f_i\) are known pdf functions (that we can sample from), and \(\alpha_i \in [0,1]\), \(\sum_{i = 1}^N \alpha_i = 1\), \(N \in \mathbb{N} \cup \infty\).

Sampling algorithm

Draw mixture component \(m\) according to weights \(\alpha_i\)
For \(m = i\), sample Y from pdf \(f_i\)

Mixing

Claim: \(Y\) has the same distribution as \(X\).

\[ \begin{aligned} F_X(z) &= \int_{- \infty}^z f(x) dx \\ &= \int_{- \infty}^z \sum_{i = 1}^N \alpha_i f_i(x) dx && {\Huge|} \textrm{Fubini's Theorem} \\ &= \sum_{i = 1}^N \alpha_i \int_{- \infty}^z f_i(x) dx \\ &= \sum_{i = 1}^N \alpha_i F_i(z) \\ &= \sum_{i = 1}^N \alpha_i \Pr \left(Y \leq z \mid m = i\right) \\ &= \sum_{i = 1}^N \Pr \left(Y \leq z \mid m = i \right) \Pr \left(m = i \right) \\ &= \Pr \left(Y \leq z \right) \\ % &= F_Y(z) \end{aligned} \]

Rejection sampling

Have:

pdf \(f_X\) of a random variable \(X\) that is easy to evaluate but difficult to sample from
another pdf \(g_Y\) of a random variable \(Y\) with same support as \(X\), that we can sample from easily and for which \(\sup_x \, \frac{f_X(x)}{g_Y(x)} \leq c\), \(c \geq 1\)

Want: A sampling routine for \(X\)

Algorithm

Draw \(Y\) according to \(g_Y\)
Draw \(U \sim \mathcal{U}({0,1})\)
Check if \(U \leq \frac{f_X(Y)}{cg_Y(Y)}\)

3.a. If yes, \(X = Y\)

3.b. Otherwise, go to 1.

Rejection sampling

Claim: \(X\) from the rejection sampling scheme has distribution according to \(f_X\).

\[ \begin{aligned} \Pr \left(X \leq z \right) &= \Pr \left(Y \leq z {\huge \mid} U \leq \frac{f_X(Y)}{c g_Y(Y)} \right)\\ % &= \frac{\Pr \left(Y \leq z \cap \textrm{accept} \right) }{\Pr(\textrm{accept})} \\ &= \frac{\Pr \left(Y \leq z \cap U \leq \frac{f_X(Y)}{g_Y(Y)}\right) }{\Pr\left(U \leq \frac{f_X(Y)}{c g_Y(Y)}\right)} \\ &= \frac{\int_{-\infty}^z \int_{0}^{\frac{f_X(y)}{c g_Y(y)}} f_{Y, U}(y,u) du dy }{\int_{-\infty}^{\infty} \int_{0}^{\frac{f_X(y)}{c g_Y(y)}} f_{Y, U}(y,u) du dy} \\ &= \frac{\int_{-\infty}^z \int_{0}^{\frac{f_X(y)}{c g_Y(y)}} g_Y(y) du dy }{\int_{-\infty}^{\infty} \int_{0}^{\frac{f_X(y)}{c g_Y(y)}} g_Y(y) du dy} && {\Huge|} f_{Y,U}(y,u) = g_Y(y) \mathbb{1}\{ u \in (0,1) \} \\ &= \int_{-\infty}^{z} f_X(y) dy = F_X(z) && {\Huge|} \int_0^{\frac{f_X(y)}{c g_Y(y)}} g_Y(y) du = \frac{f_X(y)}{c} \end{aligned} \]

Rejection sampling

Code

library(ggplot2) 
library(patchwork) 
xs <- seq(-3, 3, length.out = 1000)

p <- function(x){exp(-2 * abs(x))} 

q <- function(x){3.5 * dcauchy(x)} 

set.seed(123) 
q_sample <- rcauchy(length(xs)) 
us <- q(q_sample) * runif(length(xs)) 
threshold <- p(q_sample) / q(q_sample)

x_point <- 0.5 
y_point <- q(x_point) 
mid_point <- p(x_point)


df <- data.frame(x = xs, y = c(p(xs), q(xs)), group = rep(c(1,2), each = length(xs)), u = us, q = q_sample, fill = factor(us < p(q_sample))) 
rej_plot <- ggplot(df, aes(x = x, y = y, group = group)) + 
                geom_point(aes(x = q, y = u, fill = fill), 
                    shape = 21, 
                    color = "black", 
                    stroke = NA, 
                    size = 1.5, 
                    alpha = .3) +
                scale_color_manual(
                    name = "Distribution",
                    values = c("1" = "steelblue", "2" = "grey"),
                    labels = c("target", "proposal")
                ) + 
                scale_fill_manual(
                    name = "Proposal sample",
                    values = c("TRUE" = "steelblue", "FALSE" = "tomato"), 
                    labels = c("TRUE" = "accept", "FALSE" = "reject")
                ) + 
                geom_segment(data = data.frame(x = x_point, y = mid_point),
                                aes(x = x, xend = x, y = 0, yend = y, group = 1), 
                                color = "steelblue", lty = 1) + 
                geom_segment(data = data.frame(x = x_point, y = mid_point, y_end = y_point),
                                aes(x = x, xend = x, y = y, yend = y_end, group = 1), 
                                color = "tomato", lty = 1) + 
                geom_line(aes(colour = factor(group))) + 
                theme_minimal() + 
                coord_cartesian(xlim = c(-3, 3)) +
                theme(plot.margin = margin(t = 0, b = 0), 
                    axis.title.x = element_blank(),  
                    axis.text.x = element_blank(), 
                    axis.title.y = element_blank(),  
                    axis.ticks.x = element_blank())  
                

rej_sample <- q_sample[us < threshold] 
bottom_hist <- ggplot(data.frame(x = rej_sample), aes(x = x)) + 
    geom_histogram(aes(y = after_stat(density)), bins = 30, color = NA, fill = "steelblue", alpha = .3) +
    geom_line(data = data.frame(x = xs, y = p(xs)), aes(x = x, y = y), color = "steelblue", alpha = .3) + 
    scale_x_continuous(limits = c(-3, 3)) +  
    scale_y_reverse(limits = c(1, 0), breaks = rev(c(0, .3, .6, .9))) + 
    theme_minimal() + 
    theme(plot.margin = margin(t = 0, b = 0), 
    axis.title.x = element_blank(),  
    axis.title.y = element_blank(),  
    axis.ticks.y = element_blank()) 

rej_plot / bottom_hist + plot_layout(guides = "collect")

Expected acceptance rate

The expected acceptance rate of the rejection sampler is \(\Pr(\text{accept}) = n/c\).

\[ \mathrm{E}[N_{\text{accept}}] = \mathrm{E} \left[\sum_{i = 1}^n \mathbb{1} \{\text{accept } X_i \}\right] = \sum_{i = 1}^n \mathrm{E}[\mathbb{1} \{\text{accept } X_i \}] = n \Pr(\text{accept } X_i) \,. \]

Acceptance Probability

\[ \begin{aligned} \Pr(\text{accept}) &= \Pr\left(U \leq \frac{f_X(Y)}{c g_Y(Y)}\right) \\ &= \int_{-\infty}^{\infty} \int_{0}^{1} \mathbb{1}\left\{u \leq \frac{f_X(y)}{c g_Y(y)}\right\} f_{Y,U}(y,u)\, du\, dy \\ &= \int_{-\infty}^{\infty} \int_{0}^{\frac{f_X(y)}{c g_Y(y)}} g_Y(y)\, du\, dy, \quad {\Huge |}\ f_{Y,U}(y,u) = g_Y(y)\,\mathbb{1}\{u \in (0,1)\} \\ &= \int_{-\infty}^{\infty} g_Y(y)\,\frac{f_X(y)}{c g_Y(y)}\, dy \\ &= \frac{1}{c} \int_{-\infty}^{\infty} f_X(y)\, dy \\ &= \frac{1}{c}. \end{aligned} \]

Questions

Question 1

Explain how you would generate random variables with density \[ f(x) = (1 - \alpha) + \alpha \frac{2}{(1 + x)^2}, \quad x \in (0, 1) \] If you want to write the code, set \(\alpha = 0.6\).

We see that the random variable \(X\) is a discrete mixture of a uniform random variable \(Y\) and another random variable \(Z\) with density \(2 / (1 + x^2)\).

Sampling scheme

Draw Bernoulli random variable \(I\) with \(\Pr(I = 1) = \alpha\)
1. If \(I = 1\), sample random variable \(Z\)
1. Else, sample uniform random variable \(Y\)

Question 1

Sampling \(Z\)

\[ F_Z(z) = \int_0^z \frac{2}{(1 + x)^2} dx = \left[ - \frac{2}{1 + x}\right]_0^z = \frac{2z}{1 + z}\]

Thus

\[ F^{-1}_Z(p) = \frac{p}{2 - p} \]

Code

alpha <- 0.6 
N <- 10^5

inv_cdf <- function(x){
   x / (2 - x)
}

set.seed(123) 
Is <- rbinom(N, 1, alpha) 
Xs <- runif(N)
Xs[Is == 1] <- inv_cdf(Xs[Is == 1])

xs <- seq(from = 0, to = 1, length.out = 1000)

pdf_fun <- function(x){
    (1 - alpha)  + alpha * 2 / (1 + x)^2
}

hist(Xs, freq = FALSE) 
lines(xs, pdf_fun(xs), col = "blue", lwd = 2)

Question 2

Generate a sample from the density \[ f(x) = \sum_{i = 1}^5 \alpha_i \frac{i^3 x^2 e^{-ix}}{2} \] for \(\alpha = (0.1, 0.2, 0.2, 0.4, 0.1)^\top\). Use the sample function for the weights. Plot a histogram and calculate the mean (with a \(95\%\) confidence interval) comparing it with the actual mean of the distribution.

Note that the densities of the mixture components are from a gamma distribution with shape \(\alpha = 3\), and rate \(\beta = i\), i.e.

\[f_i(x) = \frac{i^3 x^2 e^{-ix}}{\Gamma(3)} = f_{\Gamma}(x; \alpha = 3, \beta = i) \]

Question 2

Code

N <- 10^5 
alphas <- c(0.1,0.2,0.2,0.4,0.1) 

set.seed(123)
betas <- sample(1:5, N, prob = alphas, replace = TRUE)
Xs <- rgamma(N, 3, rate = betas) 

mixture_pdf <- function(x, alphas){
    0.5 * x^2 * (vapply(1:length(alphas), function(i){i^3 * exp(-i * x)}, rep(0., length(x))) %*% alphas)
}

hist(Xs, 
    freq = FALSE, 
    ylim = c(0, .8), 
    main = bquote("Histogram of discrete mixture of " ~ Gamma(alpha, beta) ~ " distributions"))  

xs <- seq(from = 0, to = 14, length.out = 1000) 
dens <- mixture_pdf(xs, alphas) 
lines(xs, dens, col = "red", lty = 1, lwd = 2)
lines(density(Xs), col = "blue", lty = 2, lwd = 2)
legend("topright", 
    legend = c("Emp pdf", "pdf", "Est pdf"), 
    pch = c(15, NA, NA), 
    col = c("grey", "red", "blue"), 
    lty = c(NA, 1, 2), 
    lwd = c(NA, 2, 2))

Code

m <- mean(Xs)  
sdx <- sd(Xs) / sqrt(N)
alpha <- 0.05 
z <- qnorm((1 - alpha/ 2)) 
ci <- m + c(-1, 1) * z * sdx
true_mean <- (0.1 + 0.2 / 2 + 0.2 / 3 + 0.4 / 4 + 0.1 / 5) * 3
cat("Estimated mean: ", m, ", 95% CI: ", ci, ", true mean: ", true_mean)

Estimated mean:  1.16124 , 95% CI:  1.154816 1.167664 , true mean:  1.16

Question 3

Let \(U \sim \Gamma(3,1)\), \(V \sim \Gamma(2, 1)\) independent of each other. For \(X = \frac{U}{U + V}\), give the qqplot of a sample \(X\) of size \(1000\) against \(Y \sim B(3,2)\). What do you see?

Code

N <- 1000
set.seed(123) 

U <- rgamma(N, shape = 3, rate = 1)
V <- rgamma(N, shape = 2, rate = 1)
X <- U / (U + V) 
Y <- rbeta(N, 3, 2) 

par(mfrow = c(1, 2))
qqplot(Y, X, xlab = "B(3, 2)", ylab = "X=U/(U+V)", col = "red")
abline(0, 1, col = "blue")

plot(density(X), 
    lwd = 2, 
    col = "blue", 
    xlab = "x",
    main = "Probability Density Function")
lines(density(Y), lwd = 2, col = "red")

Question 4

Generate a sample from the density \[ f(x) = \frac{1}{6} + \frac{x}{2} + x^2 + x^3, \quad 0 < x < 1 \,. \] Plot a histogram and calculate the mean (with a \(95\%\) confidence interval) comparing it with the actual mean of the distribution.

Rejection sampling

proposal distribution: \(\mathcal{U}(0,1)\)
scaling: \(c = \sup_x f(x) / q(x) = \sup_x f(x) = 8 / 3\)
Algorithm:
1. Draw \(U_1, U_2\)
2. If \(U_2 \leq f(U_1) / c\), accept, else reject

Question 4

Rejection sampling

Code

pdfs <- function(x) {
  1 / 6 + x / 2 + x^2 + x^3
}

c <- 8 / 3 
N <- 10^5

set.seed(123) 
U_1 <- runif(N)  
U_2 <- runif(N)
threshold <- pdfs(U_1) / c
X <- U_1[U_2 <= threshold]

hist(X, freq = FALSE, breaks = 50, main = "Rejection Sampling", xlab = "x")
curve(pdfs, from = 0, to = 1, col = "red", add = TRUE, lwd = 2)

Question 4

Mixture sampling

\[ f(x) = \sum_{i = 1}^{4} \alpha_i f_i(x) = \sum_{i = 1}^{4} \alpha_i (c_i x_i^{i-1}) \]

\[ F_i(z) = \int_0^z x^{i-1} dx = \frac{1}{i}[x^i]^{z}_0 = \frac{z^i}{i} \]

To sum to one, we need to choose \(c_i\) such that

\[ F(1) = 1 \implies c_i = i \]

Hence

\[ f(x) = \sum_{i = 1}^{4} \alpha_i (i x_i^{i-1}) \]

Question 4

Mixture sampling

Matching mixture weights

\[ [1 \cdot \alpha_1, 2 \cdot \alpha_2,3 \cdot \alpha_3,4 \cdot \alpha_4]^\top = [1/6, 1/2, 1, 1]^\top \implies \alpha = [1/6, 1/4, 1/3, 1/4]^\top \,. \]

Question 4

Inverse transform sampling

\[F_i(z) = \int_0^z i x^{i-1} dx = [x^i]^{z}_0 = z^i\]

and

\[ F_i^{-1}(p) = p^{1 / i} \]

Alogrithm:

sample \(i\) from \(1,\ldots, 4\) according to \(\alpha\)
draw \(U \sim \mathcal{U}(0,1)\)
\(X = U^{1 / i}\)

Question 4

Mixture sampling

Code

N <- 10^5 
set.seed(123)
alphas <- c(1 / 6, 1 / 4, 1 / 3, 1 / 4) 
Is <- sample(1:4, N, replace = TRUE, prob = alphas)
Us <- runif(N) 
Xs <- Us^(1 / Is)

hist(Xs, freq = FALSE, breaks = 50, main = "Mixture Sampling", xlab = "x")
curve(pdfs, from = 0, to = 1, col = "red", add = TRUE, lwd = 2)

Question 4

Mean

\[ E_{f_i}(x) = \int_0^1 x f_i(x) dx = \int_0^1 i x^i dx = \frac{i}{i + 1} \]

Code

m <- mean(Xs)
sdx <- sd(Xs) / sqrt(N)
alpha <- 0.05 
z <- qnorm((1 - alpha/ 2)) 
ci <- m + c(-1, 1) * z * sdx
true_mean <- alphas %*% ((1 : 4) / (2 : 5))
cat("Estimated mean: ", m, ", 95% CI: ", ci, ", true mean: ", true_mean)

Estimated mean:  0.700031 , 95% CI:  0.6985485 0.7015136 , true mean:  0.7

Question 5

Using acceptance-rejection methods to generate random variable \(X\) with density \[ (1 − x)(1 + 3x), \quad 0 < x < 1 \,. \] You should try with the following envelopes:

\(1\)

\(2(1 − x)\)

\(2(1+3x) / 5\)

Which one of the three methods is best?

Question 5

Calculating upper bounds on \(f(x) / g(x)\)

\(\sup_x (1 - x) (1 + 3 x) = 4 / 3\)
\(\sup_x (1 + 3x) / 2 = 2\)
\(\sup_x 5(1 - x) / 2 = 5 / 2\)

Question 5

Code

c_1 <- 4/3
c_2 <- 2 
c_3 <- 5/2

pdf <- function(x){(1 - x) * (1 + 3 * x)} 
r_1 <- function(x){rep(c_1, length(x))} 
r_2 <- function(x){c_2 * (1 + 3 * x) / 2} 
r_3 <- function(x){c_3 * (1 - x) * 5 / 2}

xs <- seq(from = 0, to = 1, length.out = 1000)


plot(xs, pdf(xs), type = "l", xlim = c(0, 1), ylim = c(0, 2.5), main = "Acceptance region", ylab = "", xlab = "", col = "blue")

polygon(c(xs, rev(xs)),
        c(pdf(xs), rep(0, length(xs))),
        col = adjustcolor("lightblue", alpha.f = 0.5),
        border = NA)

lines(xs, pdf(xs), col = "blue", lwd = 2)
curve(r_1, from = 0, to = 1, col = "red", add = TRUE, lwd = 2)
curve(r_2, from = 0, to = 1, col = "purple", add = TRUE, lwd = 2)
curve(r_3, from = 0, to = 1, col = "green", add = TRUE, lwd = 2)

legends <- c(
  "Target",
  expression((1 - x) * (1 + 3 * x)),
  expression(frac(1 + 3 * x, 2)),
  expression(frac(5 * (1 - x), 2))
)

par(xpd = NA)
legend("topright",
       #inset = c(-0.1, 0),
       legend = legends,
       lwd = rep(2, 4),
       col = c("blue", "red", "purple", "green"),
       cex = 0.5,
       bty = "n")

Code

par(xpd = FALSE)

Question 5

Calculating CDF of proposal distribution

\(F_1(z) = z\)
\(F_2(z) = 2 \int_0^z (1 - x) dx = 2[x + 1/2 x^2]^z_0 = 2 (z - z^2 / 2)\)
\(F_3(z) = 2 / 5 \int_0^z (1 + 3x) dx = 2 [x + 3/2 x^2]^z_0 / 5\)

Inverse CDF

\(F_1^{-1}(p) = p\)
\(F_2^{-1}(p) = 1 - \sqrt{1 - p}\)
\(F_3^{-1}(p) = \{-1 + \sqrt{1 + 15p} \} / 3\)

Question 5

Code

f <- function(x){
    (1 - x) * (1 + 3 * x) 
}

f_1 <- function(x){
    1
}

F_1_inv <- function(p){
    p
} 

f_2 <- function(x){
    2 * (1 - x) 
}

F_2_inv <- function(p){
    1 - sqrt(1 - p) 
}

f_3 <- function(x){
    2 * (1 + 3 * x) / 5
}

F_3_inv <- function(p){
    (sqrt(1 + 15 * p) - 1) / 3 
}

N <- 10^5 
set.seed(123)

U_1 <- runif(N) 
U_2 <- runif(N) 

X_1 <- F_1_inv(U_1)
X_2 <- F_2_inv(U_1)
X_3 <- F_3_inv(U_1)

c_1 <- 4 / 3
c_2 <- 2 
c_3 <- 5 / 2

Y_1 <- X_1[U_2 <= f(X_1) / (f_1(X_1) * c_1)] 
Y_2 <- X_2[U_2 <= f(X_2) / (f_2(X_2) * c_2)] 
Y_3 <- X_3[U_2 <= f(X_3) / (f_3(X_3) * c_3)] 

ar_1 <- length(Y_1) / N 
ar_2 <- length(Y_2) / N 
ar_3 <- length(Y_3) / N 

th_acc_rates <- 1 / c(c_1, c_2, c_3)
emp_acc_rates <- c(ar_1, ar_2, ar_3) 

print(rbind(th_acc_rates, emp_acc_rates))

                 [,1]   [,2]    [,3]
th_acc_rates  0.75000 0.5000 0.40000
emp_acc_rates 0.75158 0.5002 0.40085

Code

curve(f, from = 0, to = 1,  col = "black", lwd = 2, lty = 2)  
lines(density(Y_1), col = "red", lwd = 2, xlab = "", main = "") 
lines(density(Y_2), col = "blue", lwd = 2) 
lines(density(Y_3), col = "green", lwd = 2) 

legend("topright", 
    legend = c("pdf", "Envelope 1", "Envelope 2", "Envelope 3"), 
    col = c("black", "red", "blue", "green"), 
    lty = c(2, 1, 1, 1), 
    lwd = rep(2, 4))

ST303 Stochastic Simumlation Class 5

Recap

Mixing

Sampling algorithm

Mixing

Rejection sampling

Algorithm

Rejection sampling

Rejection sampling

Expected acceptance rate

Acceptance Probability

Questions

Question 1

Sampling scheme

Question 1

Sampling \(Z\)

Question 2

Question 2

Question 3

Question 4

Rejection sampling

Question 4

Rejection sampling

Question 4

Mixture sampling

Question 4

Mixture sampling

Question 4

Inverse transform sampling

Question 4

Mixture sampling

Question 4

Mean

Question 5

Question 5

Calculating upper bounds on \(f(x) / g(x)\)

Question 5

Question 5

Calculating CDF of proposal distribution

Inverse CDF

Question 5

ST303 Stochastic Simumlation

Class 5