Importance sampling and Monte Carlo estimators

Worked solutions

Problem 1

1. $\mathbb{E}[\hat\mu_n] = \tfrac1n\sum_i\mathbb{E}_p[f(x_i)] = \tfrac1n\cdot n\mu = \mu$.Linearity; every $x_i$ is a draw from $p$.
2. $\operatorname{Var}(\hat\mu_n) = \tfrac1{n^2}\sum_i\operatorname{Var}_p(f(x_i)) = \dfrac{\operatorname{Var}_p(f)}n$.Independence makes the variance of the sum the sum of the variances; $\operatorname{Var}(aX) = a^2\operatorname{Var}(X)$ with $a = 1/n$.
3. The standard error $\sqrt{\operatorname{Var}_p(f)/n}$ is halved when $n$ is multiplied by $4$.$\sqrt{1/(4n)} = \tfrac12\sqrt{1/n}$.
4. $\mathbb{E}[\hat\mu_n] = \mu$, $\operatorname{Var}(\hat\mu_n) = \operatorname{Var}_p(f)/n$, standard error $\sqrt{\operatorname{Var}_p(f)/n}$; halving it takes four times the samplesMonte Carlo converges at the rate $1/\sqrt n$ in every dimension, which is its whole appeal against quadrature and its whole weakness against anything with a closed form: each extra digit of accuracy costs a hundred times the samples. The constant $\operatorname{Var}_p(f)$ is the only thing that can be changed, and Problems 2 to 10 are all ways of changing it. For a probability $\mu = \mathbb{E}_p[\mathbf{1}\{A\}]$ the variance is $\mu(1 - \mu)$ and the relative error $\sqrt{(1 - \mu)/(n\mu)}$, which for a rare event is enormous; Problem 5 makes this concrete.

Problem 2

1. $\mathbb{E}_p[f] = \sum_xp(x)f(x) = \sum_{x : q(x) > 0}q(x)\dfrac{p(x)}{q(x)}f(x)$.Multiply and divide by $q(x)$, which is allowed only where $q(x) > 0$; the terms with $q(x) = 0$ are dropped from the sum.
2. The dropped terms have $f(x)p(x) = 0$ by the support condition, so they contributed nothing.The condition says $q = 0$ only where the summand $p(x)f(x)$ is already zero.
3. $\sum_{x : q(x) > 0}q(x)\,w(x)f(x) = \mathbb{E}_q[w(X)f(X)]$.A $q$-weighted sum over the support of $q$ is an expectation under $q$.
4. $\mathbb{E}[\hat\mu_{\mathrm{IS}}] = \tfrac1n\sum_i\mathbb{E}_q[w(x_i)f(x_i)] = \mu$.Linearity with $x_i \sim q$, then step 3.
5. $\mathbb{E}_p[f] = \mathbb{E}_q[wf]$ and $\hat\mu_{\mathrm{IS}}$ is unbiased, provided $q > 0$ wherever $fp \neq 0$; if the condition fails, the estimator is biased by the missing mass $\sum_{x : q(x) = 0}p(x)f(x)$The identity is a change of measure: samples from $q$ are reweighted to look like samples from $p$, and the weight is large where $q$ undersamples $p$ and small where it oversamples. The condition is exactly what the ELBO page's $w = p(x, z)/q(z)$ and the policy-gradient page's behaviour policy need. It fails silently: a proposal that never visits a region where $p$ has mass returns an estimate that is biased by everything in that region, with no large weight to signal it, because the samples that would have carried the large weight are never drawn. Unbiased is the easy part; Problem 3 is the hard part.

Problem 3

1. $\operatorname{Var}(\hat\mu_{\mathrm{IS}}) = \dfrac1n\operatorname{Var}_q(wf) = \dfrac1n\big(\mathbb{E}_q[(wf)^2] - (\mathbb{E}_q[wf])^2\big)$.Problem 1's argument with $wf$ in place of $f$ and $q$ in place of $p$; the variance as second moment minus squared mean.
2. $\mathbb{E}_q[wf] = \mu$, so $\operatorname{Var}(\hat\mu_{\mathrm{IS}}) = \dfrac1n\big(\mathbb{E}_q[w^2f^2] - \mu^2\big)$.Problem 2.
3. $\mathbb{E}_q[w^2f^2] = \sum_xq(x)\dfrac{p(x)^2}{q(x)^2}f(x)^2 = \sum_xp(x)\dfrac{p(x)}{q(x)}f(x)^2 = \mathbb{E}_p[wf^2]$.Cancel one $q$ against one of the two in $w^2$; what remains is a $p$-weighted sum of $wf^2$.
4. $\operatorname{Var}(\hat\mu_{\mathrm{IS}}) = \dfrac1n\big(\mathbb{E}_q[w^2f^2] - \mu^2\big) = \dfrac1n\big(\mathbb{E}_p[wf^2] - \mu^2\big)$; it is infinite when $\mathbb{E}_p[wf^2] = \mathbb{E}_p[f^2p/q]$ diverges, which happens when $q$ has lighter tails than $f^2p$The variance is controlled by $\mathbb{E}_p[w f^2]$: the weight $p/q$ averaged under $p$ and amplified by $f^2$. Where $q$ is much smaller than $p$ the weight is huge and the rare samples that land there dominate the average; if $p/q$ grows fast enough in the tails the second moment is infinite and the estimator, though unbiased, has no central limit theorem, so the sample average lurches whenever a tail sample arrives and its empirical standard error is meaningless. Problem 7 computes $\mathbb{E}_p[w]$ for Gaussians and finds exactly this: a proposal narrower than the target by a factor of $\sqrt2$ in standard deviation already gives infinite variance (Mistake 3). The rule of thumb is that $q$ must have tails at least as heavy as $p$'s, and Problem 4 says what the best $q$ is.

Problem 4

1. $\mathbb{E}_q[w^2f^2] = \mathbb{E}_q\big[(w\lvert f\rvert)^2\big] \ge \big(\mathbb{E}_q[w\lvert f\rvert]\big)^2$.$f^2 = \lvert f\rvert^2$; then Jensen's inequality for the convex function $t \mapsto t^2$, the Jensen page's Problem 5(a): a second moment is at least the squared mean.
2. $\mathbb{E}_q[w\lvert f\rvert] = \mathbb{E}_p\lvert f\rvert$.Problem 2's identity applied to $\lvert f\rvert$.
3. For $q^* = \lvert f\rvert p/\mathbb{E}_p\lvert f\rvert$: $w^*\lvert f\rvert = \dfrac{p}{q^*}\lvert f\rvert = \mathbb{E}_p\lvert f\rvert$, a constant, so equality holds in step 1.$p/q^* = \mathbb{E}_p\lvert f\rvert/\lvert f\rvert$ wherever $f \neq 0$, and the $\lvert f\rvert$ cancels; a constant has zero Jensen gap. $q^*$ is a valid proposal: nonnegative, sums to one, and positive wherever $fp \neq 0$.
4. For $f \ge 0$: $\mathbb{E}_p\lvert f\rvert = \mu$, $q^* = fp/\mu$, and every sample gives $w^*(x)f(x) = \mu$ exactly, so $\operatorname{Var}(\hat\mu_{\mathrm{IS}}) = 0$.Steps 2 and 3 with $\lvert f\rvert = f$; the estimator is the constant $\mu$ whatever is drawn.
5. $\mathbb{E}_q[w^2f^2] \ge (\mathbb{E}_p\lvert f\rvert)^2$ with equality at $q^* \propto \lvert f\rvert p$; for $f \ge 0$, $q^* = fp/\mu$ has zero variance; it cannot be used because normalising it requires $\mu$, the unknownThe best proposal is not the target but the target reshaped by the integrand: sample where $\lvert f\rvert p$ is large, which for a tail probability means sampling in the tail (Problem 5) and for the ELBO means sampling near the posterior. The zero-variance proposal is a statement about direction rather than a recipe, since writing it down solves the problem, but every good proposal is an approximation to it, and the bound in step 1 gives the floor that any approximation is measured against. For $f$ taking both signs the minimum variance is $(\mathbb{E}_p\lvert f\rvert)^2 - \mu^2 > 0$, so the floor is not zero. Mistake 4 on the Jensen page is the reason $q = p$ is usually far from optimal, which Problem 5 quantifies.

Problem 5

1. (a) $f = \mathbf{1}\{X > 3\}$ has $f^2 = f$, so $\operatorname{Var}_p(f) = \mu - \mu^2 = \mu(1 - \mu) \approx 1.348\times10^{-3}$, and the relative standard error at $n = 1$ is $\sqrt{(1 - \mu)/\mu} \approx 27.2$.An indicator is Bernoulli with success probability $\mu \approx 1.350\times10^{-3}$; divide the standard error by $\mu$.
2. (b) $w(x) = \dfrac{\varphi(x)}{\varphi(x - 3)} = \exp\Big(-\dfrac{x^2}2 + \dfrac{(x - 3)^2}2\Big) = \exp\Big(\dfrac{-6x + 9}2\Big) = e^{9/2 - 3x}$.The two densities share the $1/\sqrt{2\pi}$; expand $(x - 3)^2 = x^2 - 6x + 9$ and the $x^2$ cancels.
3. $\mathbb{E}_q[w^2f^2] = \int_3^\infty\dfrac{\varphi(x)^2}{\varphi(x - 3)}\,dx = \int_3^\infty\dfrac1{\sqrt{2\pi}}\exp\Big(-x^2 + \dfrac{(x - 3)^2}2\Big)dx$.Problem 3's $\mathbb{E}_q[w^2f^2] = \int q\,(p/q)^2f^2 = \int p^2/q$ over the region where $f = 1$; one factor of $\sqrt{2\pi}$ survives.
4. $-x^2 + \dfrac{x^2 - 6x + 9}2 = -\dfrac{x^2}2 - 3x + \dfrac92 = -\dfrac{(x + 3)^2}2 + 9$.Complete the square: $-\tfrac12(x^2 + 6x) = -\tfrac12(x + 3)^2 + \tfrac92$, plus the $\tfrac92$ already there.
5. $\mathbb{E}_q[w^2f^2] = e^9\int_3^\infty\varphi(x + 3)\,dx = e^9\int_6^\infty\varphi(u)\,du = e^9\big(1 - \Phi(6)\big) \approx 8103.1\times 9.866\times10^{-10} \approx 7.994\times10^{-6}$.Substitute $u = x + 3$; the integral of $\varphi$ from $6$ is the upper tail.
6. (c) $\operatorname{Var}_q(wf) = 7.994\times10^{-6} - \mu^2 \approx 7.994\times10^{-6} - 1.822\times10^{-6} = 6.17\times10^{-6}$, and the ratio is $1.348\times10^{-3}/6.17\times10^{-6} \approx 218$; the relative standard error at $n = 1$ is $\sqrt{6.17\times10^{-6}}/\mu \approx 1.84$.Problem 3; $\mu^2 \approx (1.350\times10^{-3})^2$.
7. (a) $\operatorname{Var}_p = \mu(1 - \mu) \approx 1.348\times10^{-3}$, relative error $\approx 27.2$ per sample; (b) $w(x) = e^{9/2 - 3x}$ and $\mathbb{E}_q[w^2f^2] = e^9(1 - \Phi(6)) \approx 7.994\times10^{-6}$; (c) IS variance $\approx 6.17\times10^{-6}$ per sample, $218$ times smaller, relative error $\approx 1.84$ per samplePlain Monte Carlo sees a success once in $741$ draws and needs about $n \approx 74{,}000$ samples for a $10\%$ relative error; the shifted proposal lands half its draws in the tail with weights no larger than $e^{9/2 - 9} = e^{-4.5}$, and reaches the same accuracy with $n \approx 340$. The gain comes from sampling where $fp$ is large, Problem 4's prescription, and the shift to the threshold rather than beyond it is the standard choice because it keeps the weights bounded on the region that matters. The check evaluates both variances by quadrature and confirms the factor of $218$.

Problem 6

1. $\mathbb{E}_q[\tilde w] = \sum_xq(x)\dfrac{Zp(x)}{q(x)} = Z\sum_xp(x) = Z$.The $q$ cancels and $p$ sums to one.
2. $\mathbb{E}_q[\tilde wf] = Z\,\mathbb{E}_q[wf] = Z\mu$.$\tilde w = Zw$ and Problem 2.
3. $\hat\mu_{\mathrm{SN}} = \dfrac{\tfrac1n\sum_i\tilde w_if(x_i)}{\tfrac1n\sum_i\tilde w_i}$, and as $n \to \infty$ the numerator tends to $Z\mu$ and the denominator to $Z$, so the ratio tends to $\mu$.Divide top and bottom by $n$; each is a plain Monte Carlo average under $q$ and converges to its mean (the law of large numbers); the ratio of the limits is the limit of the ratio because $Z > 0$.
4. For $n = 1$: $\hat\mu_{\mathrm{SN}} = \dfrac{\tilde w(x_1)f(x_1)}{\tilde w(x_1)} = f(x_1)$ with $x_1 \sim q$, so $\mathbb{E}[\hat\mu_{\mathrm{SN}}] = \mathbb{E}_q[f] \neq \mu$ in general.The single weight cancels; the estimate is the integrand at a draw from the proposal, not from the target.
5. $\mathbb{E}_q[\tilde w] = Z$, $\mathbb{E}_q[\tilde wf] = Z\mu$, so $\hat\mu_{\mathrm{SN}} \to \mu$; but $\mathbb{E}[\hat\mu_{\mathrm{SN}}] \neq \mu$ for finite $n$, and for $n = 1$ it equals $\mathbb{E}_q[f]$Self-normalisation is what makes importance sampling usable for posteriors, energy-based models and any $\tilde p$ whose normaliser is unknown: the unknown $Z$ cancels in the ratio. The price is bias of order $1/n$, because the ratio of two unbiased estimates is not unbiased (the Jensen page's Mistake 3: the denominator's fluctuations do not average out through a reciprocal), and the check enumerates every pair of draws for $n = 2$ to measure it. The bias vanishes faster than the standard error, which is $O(1/\sqrt n)$, so it is harmless for large $n$ and the variance of Problem 3, with $f$ replaced by $f - \mu$, is the right measure of quality; Mistake 2 overstates the guarantee.

Problem 7

1. In one dimension, $\mathbb{E}_q[w^2] = \int\dfrac{p(x)^2}{q(x)}\,dx = \int\dfrac{(2\pi)^{-1}e^{-x^2}}{(2\pi\sigma^2)^{-1/2}e^{-x^2/(2\sigma^2)}}\,dx = \dfrac{\sigma}{\sqrt{2\pi}}\int\exp\Big(-x^2\Big(1 - \dfrac1{2\sigma^2}\Big)\Big)dx$.Problem 3's $\mathbb{E}_q[w^2] = \mathbb{E}_p[w] = \int p^2/q$; collect the constants and the exponents.
2. $\int e^{-\alpha x^2}\,dx = \sqrt{\pi/\alpha}$ for $\alpha > 0$, with $\alpha = 1 - \dfrac1{2\sigma^2} = \dfrac{2\sigma^2 - 1}{2\sigma^2}$.The Gaussian integral; it converges only when $\alpha > 0$, that is $\sigma^2 > \tfrac12$.
3. $\mathbb{E}_q[w^2] = \dfrac{\sigma}{\sqrt{2\pi}}\sqrt{\dfrac{2\pi\sigma^2}{2\sigma^2 - 1}} = \dfrac{\sigma^2}{\sqrt{2\sigma^2 - 1}}$.Multiply out; $\sqrt{\pi/\alpha} = \sqrt{2\pi\sigma^2/(2\sigma^2 - 1)}$.
4. In $d$ dimensions $p$ and $q$ factor over coordinates, so $w(x) = \prod_{j=1}^dw_j(x_j)$ and $\mathbb{E}_q[w^2] = \prod_j\mathbb{E}_q[w_j^2] = \Big(\dfrac{\sigma^2}{\sqrt{2\sigma^2 - 1}}\Big)^d$.Independent coordinates: the expectation of a product of independent factors is the product of the expectations (the variance page).
5. $\sigma^2 = 2$: $\dfrac{2}{\sqrt3} \approx 1.1547$, so $\mathbb{E}_q[w^2] \approx 4.2$ for $d = 10$ and $\approx 1.3\times10^3$ for $d = 50$.$1.1547^{10} \approx 4.21$ and $1.1547^{50} \approx 1329$.
6. $\mathbb{E}_q[w^2] = \Big(\dfrac{\sigma^2}{\sqrt{2\sigma^2 - 1}}\Big)^d$ for $\sigma^2 > \tfrac12$; for $\sigma^2 = 2$ it is $\approx 4.2$ at $d = 10$ and $\approx 1.3\times10^3$ at $d = 50$; for $\sigma^2 \le \tfrac12$ it is infiniteThe base $\sigma^2/\sqrt{2\sigma^2 - 1}$ equals $1$ only at $\sigma^2 = 1$ and exceeds $1$ on both sides, so any mismatch between proposal and target is raised to the power $d$: the weights' variance, $\mathbb{E}_q[w^2] - 1$, grows exponentially with dimension, and with it the variance of every importance-sampling estimate (Problem 3). Problem 8 turns this into an effective sample size of roughly $n/1329$ at $d = 50$, so a million samples are worth about $750$. Below $\sigma^2 = \tfrac12$ the proposal's tails are too light and the integral diverges: unbiased, infinite variance, Mistake 3. This is why importance sampling is a low-dimensional tool, why the ELBO page's single-sample bound is loose in high dimensions, and why Markov chain methods replace it when $d$ is large.

Problem 8

1. $\big(\sum_iw_i\big)^2 = \sum_iw_i^2 + \sum_{i \neq j}w_iw_j \ge \sum_iw_i^2$, so $\operatorname{ESS} \ge 1$.Expand the square; the cross terms are positive.
2. $\big(\sum_iw_i\big)^2 = \big(\sum_i 1\cdot w_i\big)^2 \le \big(\sum_i1^2\big)\big(\sum_iw_i^2\big) = n\sum_iw_i^2$, so $\operatorname{ESS} \le n$.Cauchy–Schwarz for the vectors $\mathbf{1}$ and $w$; equivalently the Jensen page's $(\mathbb{E}[W])^2 \le \mathbb{E}[W^2]$ for the uniform distribution on the $w_i$.
3. Equality in step 2 holds when $w$ is proportional to $\mathbf{1}$, all weights equal; and if $w_1 \gg w_i$ for $i > 1$ then both $(\sum_iw_i)^2$ and $\sum_iw_i^2$ are $\approx w_1^2$, so $\operatorname{ESS} \approx 1$.The Cauchy–Schwarz equality case; dominant-term approximation.
4. $\sum_iW_i^2 = \dfrac{\sum_iw_i^2}{(\sum_jw_j)^2} = \dfrac1{\operatorname{ESS}}$.Substitute $W_i = w_i/\sum_jw_j$.
5. $\dfrac{\operatorname{ESS}}n = \dfrac{\big(\tfrac1n\sum_iw_i\big)^2}{\tfrac1n\sum_iw_i^2} \to \dfrac{(\mathbb{E}_q[w])^2}{\mathbb{E}_q[w^2]} = \dfrac1{\mathbb{E}_q[w^2]}$.Divide top and bottom by $n^2$ and $n$; each average converges to its mean by the law of large numbers, and $\mathbb{E}_q[w] = 1$ for normalised $p$ and $q$ (Problem 6 with $Z = 1$).
6. $1 \le \operatorname{ESS} \le n$, equal to $n$ for equal weights and near $1$ for one dominant weight; $\operatorname{ESS} = 1/\sum_iW_i^2$; $\operatorname{ESS}/n \approx 1/\mathbb{E}_q[w^2]$The effective sample size is the number of plain Monte Carlo samples the weighted sample is worth, in the sense that the self-normalised estimator's variance is approximately $\operatorname{Var}_p(f)/\operatorname{ESS}$ rather than $\operatorname{Var}_p(f)/n$. It is the diagnostic for Problem 7's degeneracy: with $\mathbb{E}_q[w^2] \approx 1329$ at $d = 50$, a million draws have an ESS near $750$, and a histogram of the normalised weights shows a handful of samples carrying almost all the mass. Sequential Monte Carlo resamples whenever $\operatorname{ESS}$ falls below a threshold such as $n/2$. The ESS uses only the weights, not $f$, so it can report a healthy sample for an integrand that is badly estimated, which is the limitation Problem 4's $\lvert f\rvert p$ points at.