Importance sampling
Attribution
Much of this material is from Chapter 3 of Robert, Casella, and Casella (1999) and Chapter 9 of Owen (2013)
Recap
Last lecture we talked about the Accept/Reject algorithm for generating draws from a target distribution.
Importance sampling
As is our focus in the first half of the course, our target estimand is the: \[ \mu = \ExpA{h(X)}{f}. \]
Importance sampling is based on the following identity for \(X \sim f\) on the left-hand side and \(X \sim g\) on the RHS:
\[ \ExpA{h(X)}{f} = \ExpA{h(X) \frac{f(X)}{g(X)}}{g} \] The ratios \(f(x) / g(x)\) are called the importance weights.
We have to be somewhat careful about the support of the distributions \(f\) and \(g\) in order for the equality to hold.
Specifically, let \(\mathcal{D} = \{x \mid f(x) > 0\}\) be the support of \(f\) and let \(\mathcal{Q} = \{x \mid g(x) > 0\}\) be the support of \(g\). Furthermore let \(\mathcal{Q}\) be defined such that \(h(x) f(x) \neq 0 \implies g(x) > 0\).
Let’s write the RHS expectation as an integral: \[ \begin{aligned} \int_{\mathcal{Q}} \frac{f(x)}{g(x)}h(x)g(x) dx & = \int_{\mathcal{Q}} h(x)f(x) dx \\ & = \int_{\mathcal{D}} h(x) f(x) dx \\ & + \int_{\mathcal{D}^c \cap \mathcal{Q}} h(x) f(x) dx \\ & - \int_{\mathcal{D} \cap \mathcal{Q}^c} h(x) f(x) dx \\ & = \int_{\mathcal{D}} h(x) f(x) dx \end{aligned} \] The second term is zero because \(x\) is outside the support of \(f\) and the third term is zero because \(g(x)\) needs support only on the set \(h(x)f(x) \neq 0\). Provided we have chosen \(g\) with proper support, we can define an estimator \(\hat{\mu}_{\text{IS}}\): \[ \hat{\mu}_{\text{IS}} = \frac{1}{S}\sum_{i=1}^S \frac{f(X_i)}{g(X_i)}h(X_i),\, X_i \sim g \]
By the SLLN, \[ \hat{\mu}_{\text{IS}} \overset{\text{a.s.}}{\to} \ExpA{\frac{f(X)}{g(X)}h(X)}{g} = \ExpA{h(X)}{f} \] This relationship and estimator will hold for many choices of \(g\). A good choice for \(g\) is one that has low variance.
An investigation of \(\Var{\hat{\mu}_{\text{IS}}}\) shows that the only term that depends on \(g\) is the first term. \[ \Var{\hat{\mu}_{\text{IS}}} = \ExpA{\frac{f^2(X)}{g^2(X)}h^2(X)}{g} - \ExpA{\frac{f(X)}{g(X)}h(X)}{g}^2 \tag{1}\] When we write out the integral we get: \[ \int_{\mathcal{Q}} \frac{f^2(x)}{g^2(x)}h^2(x)g(x) dx = \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx \]
Now we write the variance: \[ \begin{aligned} \Var{\hat{\mu}_\text{IS}} &= \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx - \mu^2 \\ & = \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx - \mu^2 \\ & = \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx + \mu^2 - 2 \mu^2 \\ & = \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx + \int_{\mathcal{Q}} \mu^2 g(x) dx - 2\int_{\mathcal{Q}} h(x) f(x) \mu dx \\ & = \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx + \int_{\mathcal{Q}} \frac{\mu^2 g^2(x)}{g(x)} dx - 2\int_{\mathcal{Q}} \frac{h(x) f(x) g(x)\mu}{g(x)} dx \\ & = \int_{\mathcal{Q}} \frac{(f(x)h(x) - \mu g(x))^2}{g(x)}dx \end{aligned} \] This expression shows that when \(h(x) \geq 0\), and suppose \(\mu > 0\) we get that \(g^*(x) = f(x) h(x) / \mu\), which yields a variance of zero. One can show that the optimal \(g(x)\) is proportional to \(f(x)h(x)\) when \(h(x)\) takes positive and negative values.
In order for the variance of \(\hat{\mu}_{\text{IS}}\) to be finite, the term \(\int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx\) must not diverge.
Example 1 (Normal importance sampling) Let’s look at using a normal importance density with variance \(\sigma^2\) for a standard normal target and \(h(x) = x\). Let’s look at the first term in the variance. Note that in this case \(\mathcal{D} \equiv \mathcal{Q}\):
\[ \begin{aligned} \int_{\mathcal{Q}} \frac{f^2(x)}{g(x)}h^2(x)dx & = \frac{\sigma}{\sqrt{2\pi}} \int_{\R} x^2 \exp \lp -\frac{x^2}{2}(2 - \sigma^{-2}) \rp \\ & = \frac{\sigma}{\sqrt{2\pi}} \int_{\R} -\frac{2}{2-\sigma^{-2}} e^{-\frac{x^2}{2}(2-\sigma^{-2}) } \\ & = -\frac{2\sigma}{2-\sigma^{-2}}\frac{d}{d (2-\sigma^{-2})} \lp (2-\sigma^{-2})^{-1/2}\rp \\ & = \sigma\lp 2-\sigma^{-2}\rp^{-3/2} \\ \end{aligned} \] if \(\sigma^2 > \frac{1}{2}\).
Self-normalized importance sampling
When we don’t know the normalizing constant for \(f(x)\), we can’t use importance sampling because we can’t compute the importance weights. It turns out that we don’t even need to know the normalizing constant for \(g(x)\) if we use self-normalized importance sampling.
Let \(g(x) = c g_u(x)\) and let \(f(x) = d f_u(x)\) so we have \(g_u(x) = g(x) / c\) and \(f_u(x) = f(x) / d\). Then the following estimator can be used to approximate the estimand \(\mu\): \[ \tilde{\mu}_{\text{IS}} = \frac{\sum_{i=1}^S h(X_i) f_u(X_i) / g_u(X_i)}{\sum_{i=1}^S f_u(X_i) / g_u(X_i)}, X_i \sim \text{g} \]
Simple substitution shows that the estimator may be written as: \[ \begin{aligned} \tilde{\mu}_{\text{IS}} & = \frac{c / d \sum_{i=1}^S h(X_i) f(X_i) / g(X_i)}{c / d \sum_{i=1}^S f(X_i) / g(X_i)} \\ & = \frac{\sum_{i=1}^S h(X_i) f(X_i) / g(X_i)}{\sum_{i=1}^S f(X_i) / g(X_i)} \end{aligned} \] Mutliplying the numerator and the denominator by \(1 / S\) shows that the numerator is the vanilla importance sampling estimator, which we know converges a.s. to \(\mu\), and the bottom estimator converges a.s. to \[ \ExpA{\frac{f(X)}{g(X)}}{g} = 1, \]
provided that \(\mathcal{D} \equiv \mathcal{Q}\). This is a downside of the self-normalized importance sampler, namely that we can’t take advantage of the domain of \(h\) when picking \(g\).
This estimator might not seem as useful, but it is widely used in Bayesian statistics.
Example 2 (Leave-one-out cross-validation) A common question we modelers ask is how well a model can be expected to generalize to new data. While this is hard to answer as we typically have only one dataset at the time we need to make a decision about which model to use, we can approximate the answer by iteratively leaving one observation out and seeing how well our model predicts that left out observation. This is called leave-one-out cross-validation, or LOO-CV for short.
In Bayesian statistics, we often rely on approximate computation of the posterior distribution \(\pi(\theta \mid y)\). For large datsets, that would mean that we need to approximate \(\pi(\theta \mid y_{-i}), i = 1, \dots, n\). Computing this can be prohibitively expensive, but we can approximate it with importance sampling.
The importance density, \(g\), will be the posterior with all the data, \(\pi(\theta \mid y)\). We will typically know that only up to a constant, so \[ g_u(\theta) = \pi(\theta \mid y) \propto \pi(\theta) \prod_{i=1}^n p(y_i \mid \theta) \] where \(p(y \mid \theta)\) is the likelihood function. The leave-one out density, \(f_u(\theta)\) is \[ f_u(\theta) = \pi(\theta \mid y_{-i}) \propto \pi(\theta) \prod_{j \neq i}^n p(y_j \mid \theta) \] The ratio of these two \(f_u(\theta) / g_u(\theta) = (p(y_i \mid \theta))^{-1}\).
The target is the expected predictive density or \(p(y_i \mid y_{-i})\), which can be expressed as an expectation with respect to the leave-one-out posterior \[ p(y_i \mid y_{-i}) = \int_{\Theta} p(y_i \mid \theta) \pi(\theta \mid y_{-i}) d\theta \] or \[ \ExpA{p(y_i \mid \theta)}{\pi(\theta \mid y_{-i})} \] To compute this quantity we can use importance sampling, with \(h(\theta) = p(y_i \mid \theta)\).
With draws \(\theta_j \sim \pi(\theta \mid y)\) the estimator then looks like: \[ \begin{aligned} \tilde{\mu}_{\text{IS}} & = \frac{\sum_{j=1}^S p(y_i \mid \theta_j) (p(y_i \mid \theta_j))^{-1}}{\sum_{j=1}^S (p(y_i \mid \theta_j))^{-1}} \\ & = \frac{1}{\frac{1}{S} \sum_{j=1}^S \frac{1}{p(y_i \mid \theta_j)}} \end{aligned} \]
A problem with this is that the variance of the these terms may be unbounded, so the estimator is not useful without modifications.
Variance of the self-normalized IS
In order to compute the variance, we may use the delta method, assuming that the conditions needed to apply the delta method are met.
The estimator takes the form of a ratio of means: \[ \tilde{\mu}_{\text{IS}} = \frac{\frac{1}{S} \sum_{i=1}^S h(X_i) f_u(X_i) / g_u(X_i)}{\frac{1}{S} \sum_{i=1}^S f_u(X_i) / g_u(X_i)}, \] so if these means can be argued to be asymptotically joint multivariate normal, then we can apply the delta method.
Aside on the delta method for ratios
As a reminder, the multivariate delta method says that for a random vector \(Z_n\) such that \(\sqrt{n} (\frac{\sum_{i=1}^n Z_i}{n} - \mu) \overset{p}{\to} N(0, \Sigma)\) as \(n\to \infty\), where \(\Sigma = \Var{Z_1}\) then for a univariate function \(g(Z_n)\), we may derive the limiting distribution as: \[ \sqrt{n} (q(Z_n) - q(\mu)) \overset{p}{\to} N(0, (\nabla q(\mu))^T \Sigma\nabla q(\mu)) \] In our case we’ll write \(q(Z_n)\) as \(\bar{y} / \bar{x}\) and \(q(\mu) = \mu_y / \mu_x\). \[ \nabla q(\mu) = \left(-\frac{\mu_y}{\mu_x^2}, \frac{1}{\mu_x}\right)^T \]
Then we can write the variance term as: \[ \frac{\mu_y^2}{\mu_x^4} \sigma^2_{x} + \frac{\sigma^2_{y}}{\mu_x^2} - 2 \frac{\sigma_{xy} \mu_y}{\mu_x^3} \]
or, pulling out a \(\mu_x^2\) from the denominator and subbing in \(q(\mu)\) for \(\mu_y / \mu_x\) we get
\[ \frac{1}{\mu_x^2} \lp q^2(\mu) \sigma^2_{x} + \sigma^2_{y} - 2 \sigma_{xy}q(\mu)\rp \] If we sub \(m_x^2 - \mu_x^2\),\(m_y^2 - \mu_y^2\), \(m_{xy} - \mu_x \mu_y\) for the cumulants above, we’ll get:
\[ \frac{1}{\mu_x^2} \lp q^2(\mu)(m_x^2 - \mu_x^2) + m_y^2 - \mu_y^2 - 2 (m_{xy} - \mu_x \mu_y)q(\mu)\rp = \frac{q^2(\mu)m_x^2 + m_y^2 - 2 m_{xy}q(\mu)}{\mu_x^2} \] where we can write the numerator as a quadratic form: \[ \Var{\frac{\bar{y}}{\bar{x}}} \approx \frac{\Exp{(Y - \frac{\mu_y}{\mu_x} X)^2}}{\mu_x^2} \]
Delta method for self-normalized importance sampling
The delta method leads to a tidy expression for the variance, which then yields a computable variance estimator:
\[ \Var{\sqrt{S} (\frac{\bar{\frac{h f}{g}}}{\bar{\frac{f}{g}}} - \mu)} \approx \frac{\ExpA{(h(X) f(X) / g(X) - \mu f(X)/g(X))^2}{g}}{\ExpA{\frac{f(X)}{g(X)}}{g}^2} \]
or \[ \Var{\sqrt{S} \lp \frac{\bar{\frac{h f}{g}}}{\bar{\frac{f}{g}}} - \mu\rp} \approx \ExpA{\frac{f^2(X)}{g^2(X)} (h(X) - \mu)^2}{g} \]
Given that we don’t know \(f(X_i) / g(X_i)\) we have to use the SLLN again to get a computable variance estimate, which we apply to the first estimator.
We need a slightly less ungainly expression, so we’ll let \(w_u(X_i) = f_u(X_i) / g_u(X_i)\). We’ll also multiply our variance estimate by \(1/S\), which will give us the variance of the IS estimator, whereas we have derived the variance of the scaled, centered estimator. This is the proper quantity to apply a central limit theorem to, but for a given sample size \(S\), we’ll want to use \(\ExpA{\frac{f^2(X)}{g^2(X)} (h(X) - \mu)^2}{g} / S\) in our confidence set construction:
\[ \frac{\ExpA{w_u^2(X)(h(X) - \mu)^2}{g}}{S\ExpA{w_u(X)}{g}^2} \approx \frac{\frac{1}{S^2}\sum_{i=1}^S w_u^2(X_i)(h(X_i) - \tilde{\mu}_{\text{IS}})^2}{\lp \frac{1}{S} \sum_{i=1}^S w_u(X_i) \rp^2} \] Letting \[ w_n(X_i) = \frac{w_u(X_i)}{\sum_{i=1}^S w_u(X_i)} \]
we can write the final estimator for the variance as: \[ \widehat{\Var{\tilde{\mu}_{\text{IS}}}} = \sum_{i=1}^S w^2_n(X_i) (h(X_i) - \tilde{\mu}_{\text{IS}})^2 \tag{2}\]
Example 3 (LOO-CV-IS estimator) Our unnormalized importance weights are \[ w_u(\theta_j) = 1/p(y_i \mid \theta_j) \] This gives us \[ w_n(\theta_j) = \frac{1}{p(y_i \mid \theta_j)\sum_{k=1}^S 1/p(y_i \mid \theta_k)} = \frac{1}{\sum_{k\neq j}^S 1/p(y_i \mid \theta_k)} \] In this case \(h(X_i) = p(y_i \mid \theta_j)\), so we’re left with \[ \widehat{\Var{\tilde{\mu}_{\text{IS}}}} = \sum_{k=1}^S \lp\sum_{j\neq k}^S \frac{1}{p(y_i \mid \theta_j)}\rp^{-2}(p(y_i \mid \theta_k) - \tilde{\mu}_{\text{IS}})^2 \]