Answers.Body: <p>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. </p>
  
  <p>Consider our posterior density <span class="math-container">$\pi$</span> is only known up to some constant, i.e. <span class="math-container">$\pi(x) = K \tilde{\pi}(x)$</span>, where <span class="math-container">$K$</span> is some constant and we only have <span class="math-container">$\tilde{\pi}$</span>. Then by importance sampling we can evaluate the expectation of <span class="math-container">$X$</span> (or any function thereof) as follows by using a proposal density <span class="math-container">$q$</span>: </p>
  
  <p><span class="math-container">\begin{align}
   \mathbb{E}_\pi[X] &amp; = \int_\mathbb{R} x \frac{\pi(x)}{q(x)}q(x)dx \; ; \\
  &amp; = \frac{\int_\mathbb{R} x \frac{\pi(x)}{q(x)}q(x)dx}{\int_\mathbb{R}\frac{\pi(x)q(x)}{q(x)}dx} \; ;\\
  &amp; = \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} \; ; \\
  &amp; = \frac{\mathbb{E}_q[xw(x)]}{\mathbb{E}_q[w(x)]} \; ;
  \end{align}</span>
  where <span class="math-container">$w(x) = \frac{\tilde{\pi}(x)}{q(x)}$</span>. Note that on line two we have not done anything crazy - as <span class="math-container">$\pi$</span> is a density we know that it integrates to one and then we multiply the integral by <span class="math-container">$1 = \frac{q(x)}{q(x)}$</span>. The thing to notice is that the if we were to write <span class="math-container">$\pi(x) = K \tilde{\pi}(x)$</span> then the constants <span class="math-container">$K$</span> in the integrals would cancel, and so we have our result. </p>
  
  <p>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.</p>
  
  <p>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. </p>
  
  <p>edit: as requested I have added a concrete example</p>
  
  <p>As an example to make things concrete - suppose we have <span class="math-container">$Y_i | \theta \sim \text{Poisson}(\theta)$</span> and we are interested in <span class="math-container">$\theta \in (0, \infty)$</span>. The likelihood function for the Poisson distribution is 
      <span class="math-container">$$ f(\textbf{y} | \theta) = \prod\limits_{i=1}^n \frac{\theta^{y_i}\exp(-\theta)}{y_i\!}\;.$$</span></p>
  
  <p>We can then assign a gamma prior to <span class="math-container">$\theta$</span>, that is we say that <span class="math-container">$\theta \sim \text{Gamma}(a,b)$</span> with density 
      <span class="math-container">$$\pi(\theta) \propto \theta^{a-1} \exp(-b\theta)\;.$$</span></p>
  
  <p>By applying Bayes rule our posterior is then 
      <span class="math-container">\begin{align}
      \pi(\theta|\textbf{y}) &amp; \propto f(\textbf{y} | \theta) \pi(\theta) \\
      &amp; = \prod\limits_{i=1}^n \frac{\theta^{y_i}\exp(-\theta)}{y_i\!} \times \theta^{a-1} \exp(-b\theta) \\
      &amp; = \theta^{\sum\limits_{i=1}^n y_i + a - 1} \exp(-[n+b]\theta)\;.
  \end{align}</span>
  Now we know that this is the kernel of a Gamma(<span class="math-container">$\sum\limits_{i=1}^n y_i + a$</span>, <span class="math-container">$n+b$</span>) 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. </p>
  
  <p>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 <span class="math-container">$X$</span> from earlier would be <span class="math-container">$\theta$</span> in this example.</p>
  
  <p>Below is some Python code to further demonstrate this. </p>
  
  <pre><code>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)
  </code></pre>
  
  <p>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 <span class="math-container">$\theta$</span> that I used to simulate the data from) so importance sampling approximates the expectation well. </p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-20T22:45:08.730Z
Answers.LastActivityDate: 2020-05-21T09:33:40.807Z
Answers.LastEditDate: 2020-05-21T09:33:40.807Z
Answers.UserDisplayName: David
Answers.Body: <p>The rationale behind importance sampling is that <span class="math-container">$q(x)$</span> is difficult to sample from but <strong>easy to evaluate</strong>. Or at least you can easily evaluate some <span class="math-container">$\tilde{q}$</span> such that:
  <span class="math-container">$$
  \tilde{q}(z) = Zq(z)
  $$</span>
  where <span class="math-container">$Z$</span> (scalar) might be unknown. The geometrical example would be here e.g. sampling uniformly from an area under the curve <span class="math-container">$q(x)$</span> (in general it's not easy).</p>
  
  <blockquote>
    <p><em>Because if we know <span class="math-container">$q(x)$</span>, then we can compute the expectation directly.</em></p>
  </blockquote>
  
  <p>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 <span class="math-container">$q$</span>. </p>
  
  <p>Say you want to compute an expectation of <span class="math-container">$x$</span>, <span class="math-container">$E[x]$</span>. For this you need to calculate the following integral:
  <span class="math-container">$$
  E[x] = \int{xq(x)dx}
  $$</span>
  where <span class="math-container">$q$</span> is a probability distribution of <span class="math-container">$x$</span> for which you have an expression - so you can evaluate <span class="math-container">$q(x)$</span> (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 <span class="math-container">$q$</span> (as per example above, e.g. generating samples from the area under the curve <span class="math-container">$q(x)$</span> uniformly). </p>
  
  <p>What you can do is to calculate an expectation under a simple distribution <span class="math-container">$p$</span> (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:
  <span class="math-container">$$
  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]
  $$</span>
  (added index <span class="math-container">$p$</span> and <span class="math-container">$q$</span> to expectation to denote the sampling distribution). Now you can approximate the last expectation by Monte Carlo:
  <span class="math-container">$$
  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)
  $$</span></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-05T13:47:16.747Z
Answers.Text: Comments are not for extended discussion; this conversation has been [moved to chat](https://chat.stackexchange.com/rooms/118058/discussion-on-answer-by-tomasz-bartkowiak-how-can-we-compute-the-ratio-between-t).
Answers.UserDisplayName: nbro
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-20T22:18:21.540Z
Answers.LastActivityDate: 2020-05-21T08:26:13.260Z
Answers.LastEditDate: 2020-05-21T08:26:13.260Z
Answers.UserDisplayName: Tomasz Bartkowiak
Body: <p>Here is my understanding of importance sampling. If we have two distributions <span class="math-container">$p(x)$</span> and <span class="math-container">$q(x)$</span>, where we have a way of sampling from <span class="math-container">$p(x)$</span> but not from <span class="math-container">$q(x)$</span>, but we want to compute the expectation wrt <span class="math-container">$q(x)$</span>, then we use importance sampling. </p>
  
  <p>The formula goes as follows: </p>
  
  <p><span class="math-container">$$
  E_q[x] = E_p\Big[x\frac{q(x)}{p(x)}\Big]
  $$</span></p>
  
  <p>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 <span class="math-container">$q(x)$</span>, how can we compute the ratio <span class="math-container">$\frac{q(x)}{p(x)}$</span>? </p>
  
  <p>Because if we know <span class="math-container">$q(x)$</span>, then we can compute the expectation directly. </p>
  
  <p>I am sure I am missing something here, but I am not sure what. Can someone help me understand this?  </p>
  
Tags: <reinforcement-learning><monte-carlo-methods><importance-sampling>
Title: How can we compute the ratio between the distributions if we don't know one of the distributions?
UserDisplayName: pecey
ThreadText: # How can we compute the ratio between the distributions if we don't know one of the distributions?
  
  (Asked by: username_0 on 2020-05-20)
  
  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?
  
  **Tags:** reinforcement-learning, monte-carlo-methods, importance-sampling
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2020-05-20. Score: 3 | Accepted answer
  
  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)
  $$
  
  > **Comments:**
  > **username_3**: Comments are not for extended discussion; this conversation has been [moved to chat](https://chat.stackexchange.com/rooms/118058/discussion-on-answer-by-tomasz-bartkowiak-how-can-we-compute-the-ratio-between-t).  
  
  ---
  
  ### Answer by username_1 on 2020-05-20. Score: 3
  
  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.
  
  
  ---
  
  