# 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

1. Draw mixture component \(m\) according to weights \(\alpha_i\)
2. 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

1. Draw \(Y\) according to \(g_Y\)

2. Draw \(U \sim \mathcal{U}({0,1})\)

3. 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

1. Draw Bernoulli random variable \(I\) with \(\Pr(I = 1) = \alpha\)
2. 1. If \(I = 1\), sample random variable \(Z\)
3. 2. 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:

1. sample \(i\) from \(1,\ldots, 4\) according to \(\alpha\)
2. draw \(U \sim \mathcal{U}(0,1)\)
3. \(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. \(1\)
> 2. \(2(1 − x)\)
> 3. \(2(1+3x) / 5\)

> Which one of the three methods is best?

## Question 5

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

1. \(\sup_x (1 - x) (1 + 3 x) = 4 / 3\)
2. \(\sup_x (1 + 3x) / 2 = 2\)
3. \(\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

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

### Inverse CDF

1. \(F_1^{-1}(p) = p\)
2. \(F_2^{-1}(p) = 1 - \sqrt{1 - p}\)
3. \(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))