library(ggplot2)
library(patchwork)
set.seed(536)
theme_set(theme_bw(base_size = 14))
pal <- c(red = "#D62828", blue = "#2E86AB", green = "#06A77D", dark = "#023047", grey = "#6a6a6a")
find_pdflatex <- function() {
cand <- c(Sys.which("pdflatex"),
path.expand("~/Library/TinyTeX/bin/universal-darwin/pdflatex"),
file.path(Sys.getenv("APPDATA"), "TinyTeX/bin/windows/pdflatex.exe"),
path.expand("~/.TinyTeX/bin/x86_64-linux/pdflatex"))
cand <- cand[nzchar(cand) & file.exists(cand)]
if (!length(cand)) stop("pdflatex not found; install TinyTeX with tinytex::install_tinytex()")
cand[[1]]
}
tikz_png <- function(p, png_file, width = 9, height = 5, dpi = 150) {
for (pkg in c("tikzDevice", "pdftools", "png"))
if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
options(tikzLatex = find_pdflatex(),
tikzDocumentDeclaration = r"(\documentclass[10pt]{article})",
tikzLatexPackages = c(r"(\usepackage{tikz})",
r"(\usepackage[active,tightpage,psfixbb]{preview})",
r"(\PreviewEnvironment{pgfpicture})",
r"(\setlength\PreviewBorder{0pt})",
r"(\usepackage[T1]{fontenc})", r"(\usepackage{lmodern})",
r"(\usepackage{amsmath})"))
work <- file.path(tempdir(), "tikz_png"); dir.create(work, showWarnings = FALSE)
tex <- file.path(work, "fig.tex"); pdf <- file.path(work, "fig.pdf")
tikzDevice::tikz(tex, width = width, height = height, standAlone = TRUE)
print(p); dev.off()
ok <- system2(getOption("tikzLatex"),
c("-interaction=batchmode", "-halt-on-error",
paste0("-output-directory=", shQuote(work)), shQuote(tex)),
stdout = FALSE, stderr = FALSE)
if (ok != 0 || !file.exists(pdf)) stop("pdflatex failed; see ", file.path(work, "fig.log"))
dir.create(dirname(png_file), recursive = TRUE, showWarnings = FALSE)
png::writePNG(pdftools::pdf_render_page(pdf, page = 1, dpi = dpi), png_file)
png_file
}
# tick labels in math mode, so minus signs and decimals are set by LaTeX
tex_num <- function(v) ifelse(is.na(v), NA,
paste0("$", format(v, scientific = FALSE, drop0trailing = TRUE, trim = TRUE), "$"))
tex_theme <- theme(text = element_text(size = 17), plot.title = element_text(hjust = 0.5),
legend.position = "top", legend.key.width = unit(1.2, "cm"))Lab 4: Importance Sampling and Variance Reduction
STAT 536 — Statistical Computing
Learning Objectives
By the end of today’s lab, you will:
- Estimate one integral with several sampling densities and explain why their standard errors differ
- Write the same integrand as \(h \cdot f\) in more than one way, and check in code that the estimator does not change
- Recognize an importance sampler whose reported standard error cannot be trusted
- Compare antithetic pairs, a control variate, and stratified sampling at an equal number of function evaluations
- Predict when antithetic pairs and a control variate stop working
Setup
Figure labels are typeset by LaTeX. tikz_png() draws a ggplot through tikzDevice, compiles it with pdflatex, and returns a PNG. It needs a LaTeX installation (TinyTeX is enough), and it installs the three R packages it uses.
The Identity
Recall how we defined Monte Carlo integration in week 1? We wrote an integral as an expectation under the uniform, and averaged. The same move works with any density. For a function \(u\) and a density \(g\) with \(g(x) > 0\) wherever \(u(x) \neq 0\), we have
\[\theta = \int u(x)\,dx = \int \frac{u(x)}{g(x)}\,g(x)\,dx = \E_g\!\left[\frac{u(X)}{g(X)}\right], \qquad \hat\theta_n = \frac{1}{n}\sum_{i=1}^n \frac{u(X_i)}{g(X_i)}, \quad X_i \iid g.\]
Essentially, this says that every integral is an expectation, and we choose the distribution it is taken under.
The estimator is unbiased for every such \(g\). Its variance is
\[\Var[\hat\theta_n] = \frac{1}{n}\left(\int \frac{u(x)^2}{g(x)}\,dx - \theta^2\right).\]
Thus \(g\) does not change what we estimate. It changes how precisely we estimate it.
| Where it appeared | Integrand \(u\) | Sampling density \(g\) | Averaged quantity \(u/g\) |
|---|---|---|---|
| Week 1 | \(u\) on \((a,b)\) | \(\Uniform(a,b)\) | \((b-a)\,u\) |
| Week 4, basic Monte Carlo | \(h \cdot f\) | \(f\) | \(h\) |
| Week 4, importance sampling | \(h \cdot f\) | a chosen \(g\) | \(h \cdot f/g\) |
In the last two rows \(f\) is part of the question. It is a model now, and a posterior later in the course. The density \(g\) is part of the method.
Part 1: One Integral, Four Sampling Densities
The target is
\[\theta = \int_0^1 \frac{e^{-x}}{1+x^2}\,dx,\]
an integrand that Rizzo, Statistical Computing with R, uses to compare importance functions. It has no elementary antiderivative, so we take a reference value from integrate().
[1] 0.5247971Exercise 1.1: Uniform Draws
Take \(g = \Uniform(0,1)\), which gives the week 1 estimator. The function mc_ratio() takes draws x from a density g and returns the estimate, its standard error, and the effective sample size. We use it in every exercise of this part.
Recall the effective sample size from the week 4 slides? We have \(n_{\text{eff}} = (\sum_i w_i)^2 / \sum_i w_i^2\), with \(w_i = u(X_i)/g(X_i)\). Essentially, it says how many equally weighted draws the sample is worth. It equals \(n\) when every \(w_i\) is the same, and it falls as a few of them dominate. We report it for every density in this part.
set.seed(536)
U <- runif(n)
# TODO: effective sample size of the weights w
ess <- function(w) ___________
mc_ratio <- function(x, g) {
# TODO: the quantity being averaged
w <- ___________
c(estimate = mean(w), se = sd(w) / sqrt(length(w)), ess = ess(w))
}
g0 <- function(x) rep(1, length(x))
res0 <- mc_ratio(U, g0)
print(round(res0, 5))Report. Your estimate, standard error, and effective sample size, and how many standard errors the estimate is from theta.
Exercise 1.2: Two Better Densities
Recall the inverse transform method from week 2? Two densities on \((0,1)\) follow the shape of \(u\) more closely than the uniform does:
\[g_1(x) = \frac{e^{-x}}{1 - e^{-1}}, \qquad g_2(x) = \frac{4}{\pi(1+x^2)}.\]
Find each CDF, invert it, and sample from the uniforms U you already have.
g1 <- function(x) exp(-x) / (1 - exp(-1))
g2 <- function(x) 4 / (pi * (1 + x^2))
# TODO: inverse-transform draws from g1 and from g2, both built from U
x1 <- ___________
x2 <- ___________
res1 <- mc_ratio(x1, g1)
res2 <- mc_ratio(x2, g2)
tab12 <- rbind(uniform = res0, g1 = res1, g2 = res2)
print(round(tab12, 5))
# draws the uniform needs to match each standard error
print(round((tab12[, "se"] / tab12["uniform", "se"])^2, 3))Check. Confirm your samplers before you trust the table. range(x1) and range(x2) must lie inside \((0,1)\), and hist(x1, freq = FALSE); curve(g1, add = TRUE) must agree.
Report. The table. Then, quoting your own numbers: which density gives the smallest standard error, and what fraction of the uniform’s draws does it need for the same precision? Do the effective sample sizes rank the three densities the same way?
The next chunk draws the three densities next to \(u/\theta\). It is complete. Run it as is.
xs <- seq(0, 1, length.out = 400)
dens_df <- rbind(
data.frame(x = xs, y = u(xs) / theta, g = "gstar"),
data.frame(x = xs, y = g1(xs), g = "g1"),
data.frame(x = xs, y = g2(xs), g = "g2"),
data.frame(x = xs, y = g0(xs), g = "unif")
)
dens_df$g <- factor(dens_df$g, levels = c("gstar", "g1", "g2", "unif"))
dens_lab <- c(r"($g^* = u/\theta$)", "$g_1$", "$g_2$", "uniform")
p_dens <- ggplot(dens_df, aes(x, y, colour = g, linetype = g)) +
geom_line(linewidth = 1.6) +
scale_colour_manual(values = unname(pal[c("dark", "green", "blue", "red")]), labels = dens_lab) +
scale_linetype_manual(values = c("solid", "dashed", "dashed", "dotted"), labels = dens_lab) +
scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
labs(x = "$x$", y = "density", colour = NULL, linetype = NULL,
title = r"(The closer $g$ sits to $u/\theta$, the smaller the variance)") +
tex_theme
knitr::include_graphics(tikz_png(p_dens, "figures/lab04/p1-densities.png"))Question. On \((0,1)\) the integrand \(u\) falls from \(1\) to \(0.184\). Work out the range of \(u/g_1\) and of \(u/g_2\) by hand. How does the ranking of those ranges compare with the ranking of your standard errors? Which density would give a standard error of zero, and why can we not use it?
Exercise 1.3: The Same Integrand as \(h \cdot f\)
The integrand factors as \(h \cdot f\), with \(f\) a density, in more than one way. Fill in the two missing functions.
| Reading | Density \(f\) | Function \(h\) |
|---|---|---|
| A | \(f_0(x) = 1\) | \(h_0(x) = e^{-x}/(1+x^2)\) |
| B | \(f_1(x) = e^{-x}/(1 - e^{-1})\) | \(h_1(x) = \;?\) |
| C | \(f_2(x) = 4/\{\pi(1+x^2)\}\) | \(h_2(x) = \;?\) |
Each row says \(\theta = \E_f[h(X)]\) for a different pair \((f, h)\). Basic Monte Carlo under reading B draws from \(f_1\) and averages \(h_1\).
# TODO: h1 and h2 from the table
h1 <- function(x) ___________
h2 <- function(x) ___________
# basic Monte Carlo under readings B and C, on the draws from Exercise 1.2
print(c(reading_B = mean(h1(x1)), from_1.2 = unname(res1["estimate"])))
print(c(reading_C = mean(h2(x2)), from_1.2 = unname(res2["estimate"])))Now take reading B as the question, \(\theta = \E_{f_1}[h_1(X)]\), and suppose we prefer to draw from \(g_2\). The week 4 importance sampling estimator averages \(h_1(X_i)\,f_1(X_i)/g_2(X_i)\) over \(X_i \sim g_2\).
Report. The three pairs of numbers. In one sentence, say why each pair agrees to every printed digit.
Question. The estimator uses only \(u\) and the density the draws came from. Why, then, does a statistics problem usually write its target as \(\E_f[h(X)]\), with \(f\) named separately?
Exercise 1.4: A Density That Fails
Take \(g_3(x) = 3x^2\) on \((0,1)\). Its CDF is \(x^3\), thus \(X = U^{1/3}\). We have \(g_3 > 0\) on \((0,1)\), so the identity holds and the estimator is unbiased.
idx <- seq(10, n, by = 10)
run_df <- rbind(
data.frame(n = idx, est = (cumsum(u(U)) / seq_len(n))[idx], g = "unif"),
data.frame(n = idx, est = (cumsum(w3) / seq_len(n))[idx], g = "g3")
)
p_run <- ggplot(run_df, aes(n, est, colour = g)) +
geom_hline(yintercept = theta, linetype = "dashed") +
geom_line(linewidth = 1.2) +
scale_colour_manual(values = c(g3 = pal[["red"]], unif = pal[["blue"]]),
labels = c(g3 = "$g_3(x) = 3x^2$", unif = "uniform")) +
scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
coord_cartesian(ylim = theta + c(-0.12, 0.12)) +
labs(x = "$n$", y = r"($\hat\theta_n$)", colour = NULL,
title = "Same uniforms, two sampling densities") +
tex_theme
knitr::include_graphics(tikz_png(p_run, "figures/lab04/p1-running.png"))One run cannot show whether the reported standard error can be trusted. The next chunk repeats the whole experiment 500 times under the uniform and under \(g_3\). It records how often the nominal 95% interval contains theta. It is complete. Run it as is.
set.seed(536)
R <- 500
one_run <- function(sampler, g) {
r <- mc_ratio(sampler(n), g)
c(r, covered = unname(abs(r["estimate"] - theta) <= 1.96 * r["se"]))
}
rep_unif <- replicate(R, one_run(function(n) runif(n), g0))
rep_g3 <- replicate(R, one_run(function(n) runif(n)^(1/3), g3))
summarize_reps <- function(r) c(
sd_of_estimates = sd(r["estimate", ]),
median_reported_se = median(r["se", ]),
coverage_95 = mean(r["covered", ]),
frac_below_theta = mean(r["estimate", ] < theta)
)
cov_tab <- rbind(uniform = summarize_reps(rep_unif), g3 = summarize_reps(rep_g3))
print(round(cov_tab, 4))Report. From the single run: the estimate, its reported standard error, the effective sample size next to those of the other three densities, the largest value of \(u/g_3\), and the running-estimate plot. From the repeated runs: the coverage of the 95% interval under each density, and the typical reported standard error next to the actual spread of the estimates.
Question. Near \(x = 0\) the integrand \(u(x)\) is close to \(1\). Use that to decide whether \(\int_0^1 u(x)^2/g_3(x)\,dx\) is finite. What does your answer say about the variance of the estimator, about the central limit theorem, and about the number sd(w3) / sqrt(n)? Recall the strong law from week 1? The estimator is still consistent. What does the strong law need, and does it hold here?
Question. Most of the repeated runs land on one side of theta. Which side, and why? State the rule for choosing \(g\) that \(g_3\) breaks.
Part 2: Same Budget, Four Estimators
The target is \(\theta = \int_0^1 e^x\,dx = e - 1\), the example on the antithetic slide. Every estimator below uses exactly \(n = 1000\) evaluations of the integrand. We measure the variance of an estimator directly: we run it \(R = 10000\) times and take the variance of the results. The function rep_var() does that.
The sampling density is the uniform throughout this part. Thus the averaged quantity \(u/g\) is the integrand itself, and we call it \(h\) as the variance reduction slides do.
A simulated variance from \(R\) replications has relative standard error about \(\sqrt{2/R}\). A ratio of two of them has relative standard error about \(2/\sqrt{R}\), which is 2% here. Thus a simulated ratio within 4% of its analytic target agrees with it.
For the analytic targets, we have \(\Var(e^U) = 0.2420\), \(\Cov(e^U, e^{1-U}) = -0.2342\), \(\Cov(e^U, U) = 0.1409\), and \(\Var(U) = 1/12\).
Exercise 2.1: Plain Monte Carlo
Exercise 2.2: Antithetic Pairs
Draw \(n/2\) uniforms and evaluate the integrand at \(U\) and at \(1 - U\). That is \(n\) evaluations in total.
Report. Your variance ratio. Compute \(\rho = \Cov(e^U, e^{1-U})/\Var(e^U)\) from the targets above, and compare your ratio with \(1 + \rho\).
Question. The slide gives a reduction factor of \(\frac{1}{2}(1 + \rho) = 0.016\) for this example, and you just found about twice that. Both numbers are correct. What is each one comparing?
Exercise 2.3: A Control Variate
Use \(Y = U\), whose mean \(1/2\) is known. The optimal coefficient \(c^* = \Cov(e^U, U)/\Var(U)\) is the slope of the regression of \(e^U\) on \(U\). Estimate it from the same draws.
Report. Your variance ratio, next to \(1 - \rho^2\) computed from the targets, where \(\rho\) is now the correlation of \(e^U\) and \(U\).
Exercise 2.4: Stratified Sampling
Split \((0,1)\) into \(k\) equal strata and draw \(n/k\) uniforms inside each. With equal strata and equal allocation, the stratified estimator is the plain average of all \(n\) evaluations.
set.seed(536)
strat <- function(n, fn, k) {
m <- n / k
# TODO: n draws, m of them uniform on each stratum ((j-1)/k, j/k)
# HINT: rep(0:(k - 1), each = m) gives the stratum index of every draw
U <- ___________
mean(fn(U))
}
ks <- c(4, 10, 50)
v_strat <- sapply(ks, function(k) rep_var(strat, n, integrand, k))
strat_tab <- data.frame(k = ks, ratio_to_plain = v_strat / v_plain)
print(signif(strat_tab, 3))
slope <- coef(lm(log(ratio_to_plain) ~ log(k), data = strat_tab))[2]
print(round(slope, 2))p_strat <- ggplot(strat_tab, aes(k, ratio_to_plain)) +
geom_abline(intercept = log10(1.10), slope = -2, linetype = "dashed", colour = pal[["grey"]]) +
geom_point(size = 3.5, colour = pal[["blue"]]) +
scale_x_log10(labels = tex_num) + scale_y_log10(labels = tex_num) +
labs(x = "number of strata $k$", y = "variance ratio to plain Monte Carlo",
title = "Stratified sampling against $1.10/k^2$ (dashed)") +
tex_theme
knitr::include_graphics(tikz_png(p_strat, "figures/lab04/p2-strat.png", width = 8, height = 5))Report. The table, the log-log slope, and the plot. State the rate at which the variance falls in \(k\), and which of the three methods wins at \(k = 10\) and at \(k = 50\).
Question. Ten strata per coordinate in \(d = 6\) dimensions is how many strata? What does that say about where stratification is useful?
Exercise 2.5: Predict, Then Run
Replace the integrand with \(h(x) = (x - 1/2)^2\), whose integral over \((0,1)\) is \(1/12\).
Before running anything, write down your prediction for the variance ratio of (a) antithetic pairs and (b) the control variate \(Y = U\), each against plain Monte Carlo. A sketch of \(h\) is enough to get both.
Report. Your two predictions as you wrote them down, then the three ratios. For each of antithetic pairs and the control variate, one sentence on why the ratio is what it is, referring to your numbers. Then one sentence on why stratification is unaffected.
Question. Here \(h(U)\) is a deterministic function of \(U\), and yet \(U\) is useless as a control variate. What does that tell you about what a control variate needs?
Summary and Reflection
Takeaways
- One identity: Week 1, basic Monte Carlo, and importance sampling all average \(u/g\) over draws from \(g\)
- Choice of \(g\): Every valid \(g\) gives an unbiased estimator. The flatter \(u/g\), the smaller the variance
- Failure: Look at the largest weights and the effective sample size before you trust a reported standard error
- Equal budget: Compare variance reduction methods at the same number of function evaluations
- Requirements: Each variance reduction method needs something from the integrand. Know what it is before you use the method
Self-Assessment Questions
- Why is \(\frac{1}{n}\sum u(X_i)/g(X_i)\) unbiased for every valid \(g\)?
- Which density gives zero variance, and why can we not use it?
- How do you recognize an infinite-variance importance sampler from a single run?
- Why does the antithetic reduction factor on the slide differ from the one you computed?
- When do antithetic pairs increase the variance?