$$\newcommand{\E}{\mathbb{E}}\newcommand{\Var}{\operatorname{Var}}\newcommand{\Cov}{\operatorname{Cov}}\newcommand{\iid}{\stackrel{\text{iid}}{\sim}}\newcommand{\ind}{\perp\!\!\!\perp}\newcommand{\Normal}{\mathcal{N}}\newcommand{\Uniform}{\mathcal{U}}\newcommand{\Exp}{\operatorname{Exp}}\newcommand{\Pois}{\operatorname{Pois}}\newcommand{\Ber}{\operatorname{Ber}}\newcommand{\Bin}{\operatorname{Bin}}\newcommand{\Gam}{\operatorname{Gamma}}\newcommand{\Betad}{\operatorname{Beta}}\newcommand{\argmax}{\mathop{\operatorname{arg\,max}}}\newcommand{\argmin}{\mathop{\operatorname{arg\,min}}}\newcommand{\bx}{\mathbf{x}}\newcommand{\by}{\mathbf{y}}\newcommand{\bz}{\mathbf{z}}\newcommand{\bu}{\mathbf{u}}\newcommand{\bX}{\mathbf{X}}\newcommand{\bL}{\mathbf{L}}\newcommand{\bI}{\mathbf{I}}\newcommand{\bmu}{\boldsymbol{\mu}}\newcommand{\bSigma}{\boldsymbol{\Sigma}}\newcommand{\btheta}{\boldsymbol{\theta}}\newcommand{\R}{\mathbb{R}}\newcommand{\one}{\mathbf{1}}$$

Lab 1: Monte Carlo Methods and Random Variable Generation

STAT 536 — Statistical Computing

Author

Soumojit Das

Published

Fall 2026

Setup and Introduction

Load Required Libraries

library(ggplot2)
library(patchwork)
library(tidyverse)
library(microbenchmark)
library(testthat)

set.seed(536)
theme_set(theme_bw(base_size = 14))

test_exercise <- function(exercise_num, test_func) {
  tryCatch({
    test_func()
    print(paste("Exercise", exercise_num, "passed all tests!"))
    TRUE
  }, error = function(e) {
    print(paste("Exercise", exercise_num, "failed:", e$message))
    FALSE
  })
}

Learning Objectives

By the end of this lab, you will be able to:

  1. Implement numerically stable algorithms for probability computations
  2. Apply Monte Carlo integration to estimate complex integrals
  3. Generate random variables using inverse transform and transformation methods
  4. Validate random number generators using statistical tests
  5. Visualize convergence of Monte Carlo estimators

Part 1: Numerical Stability and Foundations

Exercise 1.1: The Log-Sum-Exp Trick

When computing \(\log\left(\sum_{i=1}^n e^{x_i}\right)\) for very large or small \(x_i\), we face numerical issues.

# TODO: Implement the log-sum-exp function
log_sum_exp <- function(x) {
  # Your code here
  # Hint: Use c = max(x), ...

  # implement your solution
}

# Test your implementation
test_values <- c(-1000, -999, -998)  # These would underflow with naive method

# Naive approach (will fail)
naive_result <- tryCatch(
  log(sum(exp(test_values))),
  error = function(e) -Inf
)

# Your approach
stable_result <- log_sum_exp(test_values)

print(paste("Naive result:", naive_result))
[1] "Naive result: -Inf"
print(paste("Stable result:", stable_result))
[1] "Stable result: "
[1] "Exercise 1.1 failed: is.finite(log_sum_exp(c(-1000, -999, -998))) is not TRUE\n\n`actual`:       \n`expected`: TRUE"
[1] FALSE

Exercise 1.2: Stable Probability Computation

Compute the likelihood for 1000 Bernoulli trials with very small success probability.

# Scenario: 1000 coin flips, 3 heads, p = 0.001
n <- 1000
k <- 3
p <- 0.001

# TODO: Implement stable log-likelihood computation
compute_log_likelihood <- function(n, k, p) {
  # Your code here
  # Hint: Use log probabilities throughout
  # ...
  return(log_lik)
}

log_lik <- compute_log_likelihood(n, k, p)
print(paste("Log-likelihood:", log_lik))
print(paste("Likelihood:", exp(log_lik)))

Part 2: Monte Carlo Integration

Exercise 2.1: Estimating \(\pi\) Using Monte Carlo

Estimate \(\pi\) by computing the area of a unit circle using Monte Carlo integration.

set.seed(536)

# TODO: Implement MC estimation of pi
estimate_pi_mc <- function(n_samples) {
  # Your code here
  # Hint: Generate points in [0,1] x [0,1], check if inside quarter circle
  # Area of quarter circle = pi/4, so pi = 4 * proportion inside
  # ...
}

# Test with different sample sizes
sample_sizes <- c(100, 1000, 10000, 100000)
pi_estimates <- sapply(sample_sizes, estimate_pi_mc)

results_df <- data.frame(
  n = sample_sizes,
  estimate = pi_estimates,
  error = abs(pi_estimates - pi)
)

print(results_df)

Exercise 2.2: Visualizing Monte Carlo Convergence

Create a visualization showing how the estimate converges to \(\pi\).

set.seed(536)

# TODO: Create convergence plot
visualize_pi_convergence <- function(n_max = 10000) {
  # Your code here
  # Generate cumulative estimates and create a line plot
  # Add confidence bands using CLT
  # ...

  n_seq <- c(1:n_max)

  # Simple plot
  plot(n_seq[10:n_max], pi_estimates[10:n_max],
       type = "l", col = "blue",
       xlab = "Number of samples", ylab = expression(paste("Estimate of ", pi)),
       main = expression(paste("Monte Carlo Convergence to ", pi)))
  abline(h = pi, col = "red", lty = 2)
}

visualize_pi_convergence()

Interpretation. In one or two sentences, interpret your plot by quoting the specific numbers you obtained (e.g., the estimate and its standard error at \(n = 10^4\)), not a generic description.

Exercise 2.3: High-Dimensional Integration

Compute the volume of a unit sphere in d dimensions.

set.seed(536)

# TODO: Implement d-dimensional sphere volume estimation
estimate_sphere_volume <- function(d, n_samples = 10000) {
  # Your code here
  # Hint: Generate points in [-1,1]^d hypercube
  # Check if sum(x^2) <= 1
  # Volume = 2^d * proportion inside
  # ...
  # ...

}

# Test for different dimensions
dimensions <- c(2, 3, 5, 10)
for (d in dimensions) {
  vol_est <- estimate_sphere_volume(d)
  # True volume: pi^(d/2) / gamma(d/2 + 1)
  vol_true <- pi^(d/2) / gamma(d/2 + 1)
  print(sprintf("d=%d: Estimate=%.4f, True=%.4f, Error=%.4f",
                d, vol_est, vol_true, abs(vol_est - vol_true)))
}

Part 3: Law of Large Numbers Visualization

Exercise 3.1: LLN Visualization

Create a visualization of the Law of Large Numbers.

set.seed(536)

# TODO: Implement LLN visualization
visualize_lln <- function(n = 10000, p = 0.3) {
  # Your code here
  # Generate Bernoulli trials
  # Compute running average
  # Create plot with confidence bands

  # PLACEHOLDER - Remove this {...} and implement your solution
  samples <- {...}
  running_avg <- cumsum(samples) / (1:n)

  plot(1:n, running_avg, type = "l", col = "blue",
       xlab = "Sample size", ylab = "Running average",
       main = "Law of Large Numbers")
  abline(h = p, col = "red", lty = 2)
}

visualize_lln()

Exercise 3.2: Convergence Rate Analysis

Empirically verify the \(\sqrt{n}\) convergence rate of Monte Carlo.

set.seed(536)

# TODO: Analyze convergence rate
analyze_convergence_rate <- function() {
  # Your code here
  # For different n, estimate variance of MC estimator
  # Plot log(variance) vs log(n), should have slope -1

  # PLACEHOLDER - Remove this and implement your solution
  sample_sizes <- round(10^seq(1, 3.5, length.out = 10))
  variances <- numeric(length(sample_sizes))

  for (i in seq_along(sample_sizes)) {
    n <- sample_sizes[i]
    # Estimate integral of x^2 from 0 to 1 multiple times
    estimates <- replicate(100, mean(runif(n)^2))
    variances[i] <- var(estimates)
  }

  plot(log(sample_sizes), log(variances),
       xlab = "log(n)", ylab = "log(variance)",
       main = "Convergence Rate Analysis")
  abline(lm(log(variances) ~ log(sample_sizes)), col = "red")
}

analyze_convergence_rate()

Interpretation. In one or two sentences, interpret your plot by quoting the specific numbers you obtained (e.g., the fitted slope and the empirical variance at the largest \(n\)), not a generic description.

Part 4: Option Pricing Application

Exercise 4.1: European Call Option

Price a European call option using Monte Carlo simulation.

set.seed(536)

# Parameters
S0 <- 100      # Initial stock price
K <- 105       # Strike price
r <- 0.08      # Risk-free rate
sigma <- 0.2   # Volatility
T <- 0.25      # Time to maturity

# TODO: Implement option pricing
price_option_mc <- function(S0, K, r, sigma, T, n_paths = 10000) {
  # Your code here
  # Generate final stock prices: S(T) = S0 * exp((r - sigma^2/2)*T + sigma*sqrt(T)*Z)
  # Payoff = max(S(T) - K, 0)
  # Price = exp(-r*T) * mean(payoffs)

  # PLACEHOLDER - Remove this and implement your solution
  Z <- rnorm(n_paths)
  ST <- S0 * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)
  payoffs <- pmax(ST - K, 0)
  price <- exp(-r * T) * mean(payoffs)
  return(price)
}

option_price <- price_option_mc(S0, K, r, sigma, T)
print(paste("Option price estimate:", option_price))
[1] "Option price estimate: 2.74982950836044"

Part 5: Random Variable Generation - Inverse Transform

Exercise 5.1: Weibull Distribution

Generate samples from a Weibull distribution using inverse transform.

set.seed(536)

# TODO: Implement Weibull generation
# CDF: F(x) = 1 - exp(-(x/lambda)^k)
rweibull_inverse <- function(n, shape = 2, scale = 1) {
  # ...
  return(x)
}

# Generate samples and validate
samples <- rweibull_inverse(10000, shape = 2, scale = 1)

# Validation plot
par(mfrow = c(1, 2))
hist(samples, breaks = 50, freq = FALSE, main = "Generated Samples")
curve(dweibull(x, shape = 2, scale = 1), add = TRUE, col = "red", lwd = 2)

# Q-Q plot
qqplot(rweibull(10000, shape = 2, scale = 1), samples,
       main = "Q-Q Plot", xlab = "R's rweibull", ylab = "My implementation")
abline(0, 1, col = "red")

Exercise 5.2: Pareto Distribution

Implement Pareto distribution generator.

set.seed(536)

# TODO: Implement Pareto generation
# CDF: F(x) = 1 - (x_m/x)^alpha for x >= x_m
rpareto_inverse <- function(n, x_m = 1, alpha = 3) {
  # Your code here
  # ...

  return(x)
}

# Test and visualize
samples <- rpareto_inverse(10000, x_m = 1, alpha = 3)

# Create validation plots
hist(samples[samples < quantile(samples, 0.95)], breaks = 50,
     freq = FALSE, main = "Pareto Distribution (95% of data)")

Exercise 5.3: Custom Discrete Distribution

Generate from a custom discrete distribution efficiently.

set.seed(536)

# Weather distribution
weather_probs <- c(Sunny = 0.5, Cloudy = 0.3, Rainy = 0.2)

# TODO: Implement discrete RV generator with binary search
rweather <- function(n, probs = weather_probs) {
  # Your code here
  # Implement both linear and binary search versions

  # PLACEHOLDER - Remove this and implement your solution
  states <- names(probs)
  cum_probs <- cumsum(probs)
  u <- runif(n)
  idx <- findInterval(u, c(0, cum_probs))  # 2026 fix: the 2025 version added 1 and shifted every category
  idx <- pmin(idx, length(states))
  return(states[idx])
}

# Generate and tabulate
weather_samples <- rweather(10000)
table(weather_samples) / 10000
weather_samples
Cloudy  Rainy  Sunny 
0.2993 0.2056 0.4951

Part 6: Transformation Methods

Exercise 6.1: Student’s t-distribution

Generate Student’s t-distribution using transformation method.

set.seed(536)

# TODO: Generate t-distribution
# Use: T = Z / sqrt(V/nu) where Z ~ N(0,1), V ~ chi^2(nu)
rt_transform <- function(n, df = 5) {
  # Your code here
  # ...

  return(T)
}

# Validate
samples <- rt_transform(10000, df = 5)

# Comparison plots
par(mfrow = c(1, 2))
hist(samples, breaks = 50, freq = FALSE, main = "Generated t(5)")
curve(dt(x, df = 5), add = TRUE, col = "red", lwd = 2)

qqplot(rt(10000, df = 5), samples, main = "Q-Q Plot")
abline(0, 1, col = "red")

Exercise 6.2: Cauchy from Uniform

Generate Cauchy distribution from Uniform.

set.seed(536)

# TODO: Implement Cauchy generator
# Use: X = tan(pi * (U - 0.5))
rcauchy_transform <- function(n) {
  # Your code here
  # ...

  return(x)
}

# Test
samples <- rcauchy_transform(10000)

# Visualization (be careful with outliers!)
samples_trimmed <- samples[abs(samples) < 10]
hist(samples_trimmed, breaks = 50, freq = FALSE,
     main = "Cauchy Distribution (trimmed)")

Part 7: Box-Muller Implementation

Exercise 7.1: Step-by-Step Box-Muller

Implement and visualize Box-Muller transformation.

set.seed(536)

# TODO: Implement Box-Muller
box_muller <- function(n) {
  # Your code here
  # Remember: generates pairs of normal variates


  normals <- c(x1, x2)[1:n]

  return(list(
    normals = normals,
    U1 = u1, U2 = u2,
    X1 = x1, X2 = x2
  ))
}

# Visualize the transformation process
visualize_box_muller <- function(n = 1000) {
  # Your visualization code here
  # Show: uniform inputs -> polar transform -> normal outputs

  # PLACEHOLDER - Remove this and implement your solution
  result <- box_muller(n)

  par(mfrow = c(1, 3))
  plot(result$U1, result$U2, main = "Uniform Inputs",
       xlab = "U1", ylab = "U2", pch = 20, col = rgb(0,0,1,0.3))
  plot(result$X1, result$X2, main = "Normal Outputs",
       xlab = "X1", ylab = "X2", pch = 20, col = rgb(1,0,0,0.3))
  hist(result$normals, breaks = 30, freq = FALSE, main = "Distribution")
  curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
}

visualize_box_muller()

Exercise 7.2: Efficiency Comparison

Compare different methods for generating normal random variables.

set.seed(536)

# TODO: Benchmark different methods
benchmark_normal_methods <- function(n = 10000) {
  # Your code here
  # Compare: rnorm, Box-Muller, Marsaglia's polar

  # PLACEHOLDER - Remove this and implement your solution
  bench <- microbenchmark(
    rnorm = rnorm(n),
    box_muller = {
      n_pairs <- ceiling(n/2)
      u1 <- runif(n_pairs)
      u2 <- runif(n_pairs)
      x1 <- sqrt(-2*log(u1)) * cos(2*pi*u2)
      x2 <- sqrt(-2*log(u1)) * sin(2*pi*u2)
      c(x1, x2)[1:n]
    },
    times = 100
  )

  print(bench)
  return(bench)
}

benchmark_normal_methods(10000)
Unit: microseconds
       expr     min       lq     mean   median       uq      max neval
      rnorm 260.309 268.6115 285.6376 278.4310 293.5395  448.868   100
 box_muller 252.601 267.3815 381.5407 280.0095 299.8330 9429.590   100
Unit: microseconds
       expr     min       lq     mean   median       uq      max neval
      rnorm 260.309 268.6115 285.6376 278.4310 293.5395  448.868   100
 box_muller 252.601 267.3815 381.5407 280.0095 299.8330 9429.590   100

Part 8: Multivariate Generation

Exercise 8.1: Correlated Normal Variables

Generate bivariate normal with specified correlation.

set.seed(536)

# TODO: Generate correlated normals
rmvnorm_cholesky <- function(n, mu, Sigma) {
  # Your code here
  # Use Cholesky decomposition

  # PLACEHOLDER - Remove this and implement your solution
  d <- length(mu)
  L <- chol(Sigma)
  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
  X <- sweep(Z %*% L, 2, mu, "+")
  return(X)
}

# Test with different correlations
mu <- c(0, 0)
correlations <- c(-0.8, 0, 0.8)

par(mfrow = c(1, 3))
for (rho in correlations) {
  Sigma <- matrix(c(1, rho, rho, 1), 2, 2)
  samples <- rmvnorm_cholesky(1000, mu, Sigma)
  plot(samples, main = substitute(paste(rho, " = ", r), list(r = rho)),
       xlab = "X1", ylab = "X2", pch = 20, col = rgb(0, 0, 1, 0.3))
}

Part 9: Integration Challenge

Exercise 9.1: Joint Probability Estimation

Estimate P(X + Y > 3) where X ~ Gamma(2, 1), Y ~ Beta(3, 2).

set.seed(536)

# TODO: Implement the complete solution
estimate_joint_probability <- function(n_samples = 10000) {
  # Your code here
  # 1. Generate X using sum of exponentials
  # 2. Generate Y using gamma ratio method
  # 3. Estimate probability with confidence interval

  # PLACEHOLDER - Remove this and implement your solution
  # Generate X ~ Gamma(2, 1) as sum of two Exp(1)
  X <- rexp(n_samples) + rexp(n_samples)

  # Generate Y ~ Beta(3, 2) using gamma ratio
  gamma1 <- rgamma(n_samples, shape = 3, rate = 1)
  gamma2 <- rgamma(n_samples, shape = 2, rate = 1)
  Y <- gamma1 / (gamma1 + gamma2)

  # Estimate probability
  prob <- mean((X + Y) > 3)
  se <- sqrt(prob * (1 - prob) / n_samples)
  ci <- prob + c(-1.96, 1.96) * se

  return(list(
    probability = prob,
    se = se,
    ci = ci
  ))
}

result <- estimate_joint_probability()
print(result)
$probability
[1] 0.3154

$se
[1] 0.00464675

$ci
[1] 0.3062924 0.3245076

Part 10: Comprehensive Test Suite

Run all tests to check your implementations.

# Automated testing suite
run_all_tests <- function() {
  print("Running all tests...")

  results <- list()

  # Test numerical stability
  results$ex1.1 <- test_exercise("1.1", function() {
    test_that("Log-sum-exp is stable", {
      expect_true(is.finite(log_sum_exp(c(-1000, -999))))
    })
  })

  # Test Monte Carlo integration
  results$ex2.1 <- test_exercise("2.1", function() {
    test_that("Pi estimation is accurate", {
      pi_est <- estimate_pi_mc(100000)
      expect_equal(pi_est, pi, tolerance = 0.01)
    })
  })

  # Test random variable generation
  results$ex5.1 <- test_exercise("5.1", function() {
    test_that("Weibull generation is correct", {
      samples <- rweibull_inverse(10000, shape = 2, scale = 1)
      ks_test <- ks.test(samples, pweibull, shape = 2, scale = 1)
      expect_true(ks_test$p.value > 0.01)
    })
  })

  # Summary
  passed <- sum(unlist(results))
  total <- length(results)

  print(sprintf("%d/%d tests passed", passed, total))

  return(results)
}

# Uncomment to run tests
# run_all_tests()

Summary and Reflection

Takeaways

  1. Numerical Stability: Always use log-scale for probability computations
  2. Monte Carlo: Error decreases as \(O(1/\sqrt{n})\) regardless of dimension
  3. Validation: Always validate your generators with Q-Q plots and statistical tests
  4. Efficiency: Different methods have different computational costs

Self-Assessment Questions

  1. Why does the log-sum-exp trick work?
  2. What is the theoretical convergence rate of Monte Carlo integration?
  3. When would you use inverse transform vs. transformation methods?
  4. How do you validate a random number generator?