Why do we need importance sampling?

Question

Let's say we want to calculate the following expectation:

$$ \mathbb{E}_{z\sim p_z(z)}[f(z)] $$

One issue, is that the samples from $p_z(z)$ could be not very informative:

We see here that $f(z)$ only takes on non-zero values near the edges of $p_z(z)$.

We we can solve this by using importance sampling, which will take more samples from regions with non-zero contribution.

What I don't understand is why we would want to use importance sampling. I.e., the samples of $z$ that cause $f(z)$ to be zero are going to drag down the expectation towards zero, which seems good since indeed, for most values of $z$, $f(z)$ will be zero.

Note, I'm assuming that the expectation is a single value and roughly, the expression means: "The expected value of $f(z)$ when we randomly sample from $p_z(z)$ is $E$", where $E$ is a single number.

Answer 1

The short version is that there are two different reasons why we might use importance sampling. First, importance sampling is useful if we don't know how to sample from $p_z(z)$. Second, it is also useful in the situation you draw, where $f(z)$ has very large peaks in the tail of $p_z(z)$; it gives us a more accurate estimate of the expectation. These are two different motivations why we might use importance sampling. Both are valid. I'll explain below in depth.

---

Before I jump into the details of your question, let me start with some background, namely, the standard setting that motivates importance sampling. The standard setting is that we wish to estimate

$$\mathbb{E}_{z \sim p_z(z)}[z],$$

but we do not have a way to sample from $p_z(z)$. Since we don't have any way to sample from $p_z(z)$, there is no way to use the standard estimator for the mean (drawing a few samples from the distribution and then computing the sample average), so we need some other way to estimate this expectation.

Importance sampling gives us a solution to this problem. Importance sampling allows us to pick a distribution $q(z)$ that we do know how to sample from, and then estimate the expectation above by choosing samples distributed according to $z$. In particular, it uses the fact that

$$\mathbb{E}_{z \sim p_z(z)}[z] = \mathbb{E}_{z \sim q(z)}\left[{z p_z(x) \over q(z)}\right].$$

Notice that we do have a way to estimate the right-hand side using standard methods, by drawing samples from $q(z)$, as long as we have a way to compute $p_z(z)$ for any given $z$. There are many situations where we do know how to compute $p_z(z)$ for any given $z$ of interest, but we don't know how to efficiently draw a random sample from the distribution $z \sim p_z(z)$.

---

Now, your question asks about a slightly different situation. The above exposition talks about computing $\mathbb{E}[z]$, but your question asks about $\mathbb{E}[f(z)]$, which is slightly different. So now let's address your specific question. Why do we use importance sampling for this task? The answer is: it depends. There are two different reasons why we might use importance sampling, for this task.

First situation: If you don't have a way to sample from $p_z(z)$, but you do have a way to compute $p_z(z)$, then all of the above discussion still applies, and we can use importance sampling in this generalized setting. In particular, we use the equality

$$\mathbb{E}_{z \sim p_z(z)}[f(z)] = \mathbb{E}_{z \sim q(z)}\left[{f(z) p_z(x) \over q(z)}\right],$$

and we use standard methods to estimate the right-hand side by drawing several samples from $q(z)$. So, in this setting, the reason we use importance sampling is because there is no other choice: we can't use simpler methods, so we're forced to use this more sophisticated method.

Second situation: If we do have a way to sample from $p_z(z)$, then we don't have to use importance sampling, but it might still be to our benefit to do so. In particular, your picture gives exactly the sort of situation that can occur, but you've drawn the wrong conclusion, because your intuition has led you astray.

To help you improve your intuition, I suggest drawing the same picture, but now imagine that the height of the peak of $f(x)$ is $1000\times$ higher than the peak of $p_z(z)$. Then it is hopefully clear that sampling from $p_z(z)$ is sub-optimal: most samples will be from the center of $p_z(z)$, where $f(z)$ is close to zero. However, the expectation will be dominated by the area in the tail where $f(z)$ is large. Those samples will occur rarely if you sample from $p_z(z)$, but their value will be so huge that they will dominate. The result is that naive sampling is inefficient: you need a lot of samples to make sure you obtain a few that are in the tail. The way that this shows up is that the variance of the standard estimator is very high, and the estimator will be unbiased in theory but in practice if you don't have a very large sample size, the estimator will often give values that are very wrong.

Importance sampling offers a better method. Suppose we can find a distribution $q(z)$ that is approximately proportional to $p_z(z) f(z)$, and suppose we can sample from $q(z)$. Then it turns out that using importance sampling with $q(z)$ yields a better estimator than standard sampling from $p_z(z)$. In particular, the variance of the importance sampling estimator will be lower (its estimates will tend to be more accurate). Formally, the optimal (most efficient) estimator sets $q(z)$ to be exactly proportional to $p_z(z) f(z)$, and any distribution $q(z)$ that is close to that will be close to optimal.

Hopefully, it is intuitive why this might be so. If the peak of $f(z)$ is $1000\times$ higher than you drew it here, it's best to bias your sampling towards the area where $f(z)$ is large, and then you want to correct for that bias by dividing by some appropriate constant. That's exactly what importance sampling does, if you use a distribution $q(z)$ that is proportional to $p_z(z) f(z)$.

At this point, you might be wondering, why not simply define $q^*(z) = p_z(z) f(z)$, and then sample from $q^*(z)$? Well, two reasons. First, this is not a distribution; we need to renormalize it so it sums to one. So, the correct definition is $q^*(z) = p_z(z) f(z) / \int_t p_z(t) f(t) \; dt$. This brings us to the second reason: we might not have any efficient way to sample from the resulting $q^*(z)$. For starters, we might not have any efficient way to compute $\int_t p_z(t) f(t) \; dt$, and in addition, just because you can compute $q^*(z)$ for any $z$ of your choice is not enough to let you sample from such a distribution of your choice. Importance sampling provides a solution to all of these challenges. It doesn't require us to be able to compute the optimal distribution $q^*(z)$ or to be able to sample from it. It lets us use any convenient distribution that is close to $q^*(z)$.

Answer 2

Importance sampling returns an unbiased estimate of the expectation of interest, $E$ say, for all valid importance densities $q$. The choice of $q$ is driven by (i) the possibility to simulate from $q$, (ii) the ability to compute the importance ratio $f(\cdot)p_z(\cdot)/q(\cdot)$, and (iii) a reduced variance${}^\star$ of the Monte Carlo estimate, $$\widehat E = \frac{1}{n}\sum_{i=1}^n f(z_i) \frac{p_z(z_i)}{q(z_i)}\qquad z_i\stackrel{\text{iid}}{\sim} q\tag{1}$$ which means $q(\cdot)$ must be as similar as possible to $|f(\cdot)|p_z(\cdot)$, which formally is the optimal choice, albeit of no use since it requires the value of $E$.

Avoiding the kernel of $f$, i.e., the region of the space $\mathfrak X$ where$$f(z)=0$$is a good thing for estimating $E$, because this reduces the variance (e.g., in not producing a random number of zeroes in the Monte Carlo estimate), while not inducing bias. In other words, since $$\int f(z)p_z(z)\,\text dz=\int f(z)\mathbb I_{f(z)\ne 0}p_z(z)\,\text dz+\underbrace{\int f(z)\mathbb I_{f(z)=0}p_z(z)\,\text dz}_{\text{equal to }0}$$ using a random variable to approximate the second term is creating unnecessary noise.

For instance, take the ideal case of an arbitrary density $q(\cdot)$ that is positive everywhere and consider $q^+(\cdot)$, its truncation to the support of $f(\cdot)$, i.e., $$q^+(x) = \dfrac{q(x)\mathbb I_{f(x)\ne 0}}{\int q(z)\mathbb I_{f(z)\ne 0}\text dz}\stackrel{\text{def}}{=}\dfrac{q(x)\mathbb I_{f(x)\ne 0}}{\mathfrak Q(f(X)\ne 0)}$$ as an alternative (valid) importance density. (I assume for comparison sake that $\mathfrak Q(f(X)\ne 0)>0$ can be computed exactly.)

The importance Monte Carlo estimator associated with $q^+(\cdot)$ then writes as $$\frac{1}{n}\sum_{i=1}^n f(z_i) \frac{p_z(z_i)}{q^+(z_i)}\qquad z_i\stackrel{\text{iid}}{\sim} q^+$$ which can be rewritten as $$\frac{\mathfrak Q(f(X)\ne 0)}{n}\sum_{i=1}^n f(z_i) \frac{p_z(z_i)}{q(z_i)}\qquad z_i\stackrel{\text{iid}}{\sim} q^+\tag{2}$$ to compare with (1), which also writes as $$\frac{M}{n}\frac{1}{M}\sum_{i=1}^M f(z_i) \frac{p_z(z_i)}{q(z_i)}\qquad z_i\stackrel{\text{iid}}{\sim} q^+$$ where$$M\sim\mathcal Bin(n,\mathfrak Q(f(X)\ne 0))$$is the random number of simulations $z_i$ from $q$ (out of $n$) such that $f(z_i)\ne 0$, which are then simulations from $q^+$. One can spot two sources of additional variability in (1), against (2):

1. The random ratio $\frac{M}{n}$ has expectation $\mathfrak Q(f(X)>0)$ and variance$$\mathfrak Q(f(X)\ne 0)\mathfrak Q(f(X)=0)/n$$
2. (1) is based on less simulated random variables than (2), meaning that the variance of the empirical average is $n/M$ times larger for (1), compared with (2).

---

${}^\star$It is important to remind readers that some importance functions $q$ lead to infinite variance solutions for pairs $(f,p_z)$.

Answer 3

So long as importance sampling is done correctly, it always give an unbiased estimator for the true integral of interest --- this occurs irrespective of the generative distribution you use in the method, so long as it meets the basic requirements. So the main remaining issue of concern is the variability of the estimator. As a general rule, the closer the kernel of the generative distribution is to the integrand in the integral (i.e., to the function $x \mapsto f(x) p_X(x)$) the lower the variability of the estimator. That is largely the point of importance sampling.

Answer 4

Let me give you a very simple example. Assume that $p$ is a standard Gaussian and that $f : x \mapsto 10^{8} \cdot \mathbf{1}_{x \geq 5}$. Then the expectation you're trying to compute is $10^8 \cdot \mathbb{P}(X \geq 5) \approx 28.7$. The intuition that the expectation will be $0$ because $f$ is 0 in the important regions of $p$ is wrong. If $f$ has extremely large values in regions of small mass under $p$ then you may have an expectation that is far from $0$.

With this in mind, of course if you use standard Monte Carlo to approximate this expectation, your estimates will be 0 almost all the time even if you use $100000$ samples to compute each MC estimate. This is quite far from $28.7$.

Meanwhile, if you do importance sampling with proposal $\mathcal{N}(5, 1)$ you get much much better estimates, as I illustrate through the following snippet. I estimate the root mean square error of such estimates by computing n_rep estimates.

from torch.distributions import Normal
from scipy.stats import norm 

n_samples = 5000
n_rep = 1000
f = lambda x : 10.**8 * (x >= 5)
p = Normal(0., 1.)
q = Normal(5., 1.)

p_samples = p.sample((n_rep, n_samples))
q_samples = p_samples + 5

#Monte Carlo estimates, sampling from p
mc = f(p_samples).mean(-1)

#Importance Sampling estimates, sampling from q
log_weights = p.log_prob(q_samples) - q.log_prob(q_samples) #importance weights
IS = (f(q_samples) * log_weights.exp()).mean(-1)

print(f'Expectation value: {10**8 * (1 - norm.cdf(5))}')
print(f'RMSE Monte Carlo estimates: {mc.std()}')
print(f'RMSE IS estimates: {IS.std()}')

The code returns

Expectation value: 28.66515719235352
RMSE Monte Carlo estimates: 632.4555053710938
RMSE IS estimates: 0.9297240972518921

In short, importance sampling is very useful when you don't know how to sample from the distribution $p$ but it is also as much useful when you want to estimate "rare events" as is the case here.