Statistical Sampling

The average height of all people in Cambridge can be obtained by measuring the heights of all people in Cambridge and then averaging them:

$$E_{p\in C} [h(p)] =\frac{1}{\vert C\vert} \sum_{p\in C} h(p)$$

Alternatively, this can be approximated by sampling the population:

$$E_{p\in C} [h(p)] \approx \frac{1}{S} \sum_{s=1}^S h(p^{(s)})$$

where $\{ p^{(s)}\vert s=1,\cdots,S) \} \in C$ is a subset of all people in Cambridge $C$.

Similarly, the following expectation can be approximated by

$$E_{p(x)} f(x)=\int f(x) p(x) dx \approx \frac{1}{S} \sum_{s=1}^S f(x^{(s)})$$

where $x^{(s)}$ are drawn from the distribution $p(x)$. This is called simple Monte Carlo method. But how to sample a given distribution $p(x)$?

• Direct sampling

  We can sample from $p(x)$ directly. First find the accumlative density of $p(x)$:
  

  $$F(x)=\int_{-\infty}^x p(x')dx'$$
  
  
  whose value will range from 0 ( $x\rightarrow -\infty$) to 1 ( $x\rightarrow \infty$). We then sample $u \sim Uniform[0,1]$ (corresponding to $-\infty <x< \infty$) and map $u$ back to $x$: $x(u)=F^{-1}(u)$. This is not applicable if the true $p(x)=p^*(x)/Z_x$ is not normalized, i.e., we only have the unnormalized density $p^*(x)$ without the normalizing factor $Z_x=\int p^*(x)dx$.

• Rejection (acceptance) sampling

  Assume we want to sample a complicated distribution $p^*(x)$ which is not normalized, i.e., the normalizing factor $Z_p=\int p^*(x)dx$ (not necessarily equal to 1) is unknown.The rejection sampling method makes use of a simpler proposal distribution $q^*(x)$, also not normalized with unknown $Z_q=\int q^*(x)dx$. Assume a scalar $c$ can be obtained so that $cq^*(x)$ is the upper bound of $p^*(x)$ for all range of $x$:
  

  $$cq^*(x) > p^*(x)\;\;\;\;\;\mbox{for all $x$}$$
  
  
  We then draw two random numbers $u\sim Uniform[0,cq^*(x)]$ and $x\sim q(x)$. This $x$ is accepted as a sample from $p^*(x)$ if $u<p^*(x)$ or rejected if otherwise. As a special case, $q^*(x)$ could be a uniform distribution with constant values over $x$.

  The accepted $x's$ represent the desired distribution $p^*(x)$, because the rate of acceptance is proportional to the height of $p^*(x)$. Rejection sampling works best if the proposal distribution $q(x)$ is similar to the actual distribution $p(x)$, as otherwise many $x$ values will be rejected.

• Importance sampling

  Importance sampling is not a method for generating samples from $p(x)$, instead it is used to estimate the expectation of a function $f(x)$. Importance sampling can avoid the rejection (therefore wasted effort) in previous method. To find the following expectation of $f(x)$ with respect to a complicated distribution $p(x)=p^*(x)/Z_p$ with unknow scalar $Z_p$:
  

  $$E_{p(x)}[f(x)]=\int f(x) p(x) dx$$
  
  
  we can sample a simpler proposal distribution $q^*(x)=q(x)/Z_q$ to approximate the expectation:
  
  $$\int f(x) p(x) dx =\int f(x) \frac{p(x)}{q(x)} q(x) dx \appr... ...frac{1}{S}\sum_{s=1}^S f(x^{(s)})\frac{p(x^{(s)})}{q(x^{(s)})}$$
  
  
  where the samples $x^{(s)}$ are drawn from $q(x)$.

  If the normalizing scalars $Z_p=\int p^*(x)dx$ and $Z_q=\int q^*(x)dx$ are unknown, the above can still be obtained by
  

  $$\int f(x) p(x) dx \approx \frac{1}{S}\sum_{s=1}^S f(x^{(s)})\frac{p(x^{(s)})}{q(x^{(s)})} =\frac{Z_q}{Z_p}\frac{1}{S}\sum_{s=1}^S f(x^{(s)})\frac{p^*(x^{(s)})}{q^*(x^{(s)})}$$

  $$\approx \sum_{s=1}^S f(x^{(s)})\frac{w^{(s)}}{\sum_{s'} w^{(s')}}$$

  
  
  where $w^{(s')}$ is the weight for sample $s'$ defined as
  
  $$w^{(s')} \stackrel{\triangle}{=}\frac{p^*(x^{(s')})}{q^*(x^{(s')})}$$
  
  
  The last step is due to the following:
  
  $$\frac{Z_p}{Z_q}=\frac{1}{Z_q} \int p^*(x)dx =\frac{1}{Z_q} \i... ...p^*(x)}{q^*(x)}q(x)dx\approx \frac{1}{S}\sum_{s'=1}^S w^{(s')}$$
  
  
  Importance sampling approximates the expectation of $f(x)$ by a weighted average of the function evaluated at samples $x^{(s)}$ taken from $q(x)$, instead of $p(x)$. The function of the weights $w^{(s)}$ is to compensate for sampling the wrong distribution. When $p(x^{(s)})>q(x^{(s)})$, $w^{(s)}>1$ will increase the importance of the sample $x^{(s)}$ where $p(x^{(s)})$ is under-sampled, but when $p(x^{(s)})<q(x^{(s)})$, $w^{(s)}<1$ will derease the importance of the sample $x^{(s)}$ where $p(x^{(s)})$ is over-sampled.

Importance sampling and rejection sampling work well only if the proposal density $q(x)$ is similar to $p(x)$, which may be difficult in practice.