How can we compute the ratio between the distributions if we don't know one of the distributions?

Question

Here is my understanding of importance sampling. If we have two distributions $p(x)$ and $q(x)$, where we have a way of sampling from $p(x)$ but not from $q(x)$, but we want to compute the expectation wrt $q(x)$, then we use importance sampling.

The formula goes as follows:

$$ E_q[x] = E_p\Big[x\frac{q(x)}{p(x)}\Big] $$

The only limitation is that we need a way to compute the ratio. Now, here is what I don't understand. Without knowing the density function $q(x)$, how can we compute the ratio $\frac{q(x)}{p(x)}$?

Because if we know $q(x)$, then we can compute the expectation directly.

I am sure I am missing something here, but I am not sure what. Can someone help me understand this?

Answer 1

The rationale behind importance sampling is that $q(x)$ is difficult to sample from but easy to evaluate. Or at least you can easily evaluate some $\tilde{q}$ such that: $$ \tilde{q}(z) = Zq(z) $$ where $Z$ (scalar) might be unknown. The geometrical example would be here e.g. sampling uniformly from an area under the curve $q(x)$ (in general it's not easy).

Because if we know $q(x)$, then we can compute the expectation directly.

That's the task we're trying to solve to begin with. And calculating expectation might be hard if we can't sample efficiently from $q$.

Say you want to compute an expectation of $x$, $E[x]$. For this you need to calculate the following integral: $$ E[x] = \int{xq(x)dx} $$ where $q$ is a probability distribution of $x$ for which you have an expression - so you can evaluate $q(x)$ (up to the constant of proportionality). This integral might be hard to evaluate analytically so we need to use other methods such as Monte Carlo. Let's say it is hard to generate samples from $q$ (as per example above, e.g. generating samples from the area under the curve $q(x)$ uniformly).

What you can do is to calculate an expectation under a simple distribution $p$ (proposal distribution) which is a distribution of your choice that needs to allow you to easily sample from it (say Gaussian). Then you can rewrite your integral as: $$ E_q[x] = \int{xq(x)dx} = \int{xq(x) \color{blue}{\frac{p(x)}{p(x)}} dx} = \int{x \frac{q(x)}{\color{blue}{p(x)}} \color{blue}{p(x)} dx} = E_p \Big[{x\frac{q(x)}{p(x)}}\Big] $$ (added index $p$ and $q$ to expectation to denote the sampling distribution). Now you can approximate the last expectation by Monte Carlo: $$ E_p \Big[{x\frac{q(x)}{p(x)}}\Big] = \frac{1}{S} \sum_{s}{x^{(s)} \frac{q(x^{(s)})}{p(x^{(s)})} }, \ x^{(s)} \sim q(x) $$

Answer 2

It is common in Bayesian statistics to only know the posterior up to a constant of proportionality. This means that we can't directly sample from the posterior. However, using importance sample we are able to.

Consider our posterior density $\pi$ is only known up to some constant, i.e. $\pi(x) = K \tilde{\pi}(x)$, where $K$ is some constant and we only have $\tilde{\pi}$. Then by importance sampling we can evaluate the expectation of $X$ (or any function thereof) as follows by using a proposal density $q$:

\begin{align} \mathbb{E}_\pi[X] & = \int_\mathbb{R} x \frac{\pi(x)}{q(x)}q(x)dx \; ; \\ & = \frac{\int_\mathbb{R} x \frac{\pi(x)}{q(x)}q(x)dx}{\int_\mathbb{R}\frac{\pi(x)q(x)}{q(x)}dx} \; ;\\ & = \frac{\int_\mathbb{R} x \frac{\tilde{\pi}(x)}{q(x)}q(x)dx}{\int_\mathbb{R}\frac{\tilde{\pi}(x)q(x)}{q(x)}dx} \; ; \\ & = \frac{\mathbb{E}_q[xw(x)]}{\mathbb{E}_q[w(x)]} \; ; \end{align} where $w(x) = \frac{\tilde{\pi}(x)}{q(x)}$. Note that on line two we have not done anything crazy - as $\pi$ is a density we know that it integrates to one and then we multiply the integral by $1 = \frac{q(x)}{q(x)}$. The thing to notice is that the if we were to write $\pi(x) = K \tilde{\pi}(x)$ then the constants $K$ in the integrals would cancel, and so we have our result.

To summarise - we can sample from a distribution that is difficult/impossible to sample from (e.g. because we only know the density up to a constant of proportionality) by using importance sampling, as this allows us to calculate the importance ratio and use samples that are generated from a distribution of our choosing that is easier to sample from.

Note that importance sampling isn't just used in Bayesian statistics - for instance it could be used in Reinforcement Learning as an off policy way of sampling from the environment whilst still evaluating the value of the policy you're interested in.

edit: as requested I have added a concrete example

As an example to make things concrete - suppose we have $Y_i | \theta \sim \text{Poisson}(\theta)$ and we are interested in $\theta \in (0, \infty)$. The likelihood function for the Poisson distribution is $$ f(\textbf{y} | \theta) = \prod\limits_{i=1}^n \frac{\theta^{y_i}\exp(-\theta)}{y_i\!}\;.$$

We can then assign a gamma prior to $\theta$, that is we say that $\theta \sim \text{Gamma}(a,b)$ with density $$\pi(\theta) \propto \theta^{a-1} \exp(-b\theta)\;.$$

By applying Bayes rule our posterior is then \begin{align} \pi(\theta|\textbf{y}) & \propto f(\textbf{y} | \theta) \pi(\theta) \\ & = \prod\limits_{i=1}^n \frac{\theta^{y_i}\exp(-\theta)}{y_i\!} \times \theta^{a-1} \exp(-b\theta) \\ & = \theta^{\sum\limits_{i=1}^n y_i + a - 1} \exp(-[n+b]\theta)\;. \end{align} Now we know that this is the kernel of a Gamma($\sum\limits_{i=1}^n y_i + a$, $n+b$) distribution, but assume that we didn't know this and didn't want to calculate the normalising integral. This would mean that we are not able to calculate the mean of our posterior density, or even sample from it. This is where we can use importance sampling, for instance we could choose an Exponential(1) proposal distribution.

We would sample say 5000 times from the exponential distribution and then calculate the two expectations using MC integration and obtain an estimate for the mean of the posterior. NB that in this example $X$ from earlier would be $\theta$ in this example.

Below is some Python code to further demonstrate this.

import numpy as np

np.random.seed(1)

# sample our data
y = np.random.poisson(lam=0.5,size = 100)

# sample from proposal
samples_from_proposal = np.random.exponential(scale=1,size=5000)

# set parameters for the prior
a = 5; b = 3

def w(x, y, a, b):
    # calculates the ratio between our posterior kernel and proposal density
    pi = x ** (np.sum(y) + a - 1) * np.exp(-(len(y) + b) * x)
    q = np.exp(-x)
    return pi/q

# calculate the top expectation
top = np.mean(samples_from_proposal * w(samples_from_proposal,y,a,b))

# calculate the bottom expectation
bottom = np.mean(w(samples_from_proposal,y,a,b))

print(top/bottom)

# calculate the true mean since we knew the posterior was actually a gamma density
true_mean = (np.sum(y) + a)/(len(y) + b)
print(true_mean)

Running this you should see that the Expectation from importance sampling is 0.5434 whereas the true mean is 0.5436 (both of which are close to the true value of $\theta$ that I used to simulate the data from) so importance sampling approximates the expectation well.