Importance.lhs

% Importance Sampling
% Dominic Steinitz
% 14th August 2014

---
bibliography: Kalman.bib
---

\newcommand{\condprob} [3] {#1\left(#2 \,\vert\, #3\right)}

Importance Sampling
===================

Suppose we have an random variable $X$ with pdf $1/2\exp{-\lvert
x\rvert}$ and we wish to find its second moment numerically. However,
the [random-fu](https://hackage.haskell.org/package/random-fu) package
does not support sampling from such as distribution. We notice that

$$
\int_{-\infty}^\infty x^2 \frac{1}{2} \exp{-\lvert x\rvert} \mathrm{d}x =
\int_{-\infty}^\infty x^2 \frac{\frac{1}{2} \exp{-\lvert x\rvert}}
                               {\frac{1}{\sqrt{8\pi}}{\exp{-x^2/8}}}
                      \frac{1}{\sqrt{8\pi}}{\exp{-x^2/8}}
\,\mathrm{d}x
$$

So we can sample from ${\cal{N}}(0, 4)$ and evaluate

$$
x^2 \frac{\frac{1}{2} \exp{-\lvert x\rvert}}
         {\frac{1}{\sqrt{8\pi}}{\exp{-x^2/8}}}
$$

> {-# OPTIONS_GHC -Wall                     #-}
> {-# OPTIONS_GHC -fno-warn-name-shadowing  #-}
> {-# OPTIONS_GHC -fno-warn-type-defaults   #-}
> {-# OPTIONS_GHC -fno-warn-unused-do-bind  #-}
> {-# OPTIONS_GHC -fno-warn-missing-methods #-}
> {-# OPTIONS_GHC -fno-warn-orphans         #-}

> module Importance where

> import Control.Monad
> import Data.Random.Source.PureMT
> import Data.Random
> import Data.Random.Distribution.Binomial
> import Data.Random.Distribution.Beta
> import Control.Monad.State
> import qualified Control.Monad.Writer as W


> sampleImportance :: RVarT (W.Writer [Double]) ()
> sampleImportance = do
>   x <- rvarT $ Normal 0.0 2.0
>   let x2 = x^2
>       u = x2 * 0.5 * exp (-(abs x))
>       v = (exp ((-x2)/8)) * (recip (sqrt (8*pi)))
>       w = u / v
>   lift $ W.tell [w]
>   return ()

> runImportance :: Int -> [Double]
> runImportance n =
>   snd $
>   W.runWriter $
>   evalStateT (sample (replicateM n sampleImportance))
>              (pureMT 2)

We can run this 10,000 times to get an estimate.

    [ghci]
    import Formatting
    format (fixed 2) (sum (runImportance 10000) / 10000)

Since we know that the $n$-th moment of the exponential distribution
is $n! / \lambda^n$ where $\lambda$ is the rate (1 in this example),
the exact answer is 2 which is not too far from our estimate using
importance sampling.

The value of

$$
w(x) = \frac{1}{N}\frac{\frac{1}{2} \exp{-\lvert x\rvert}}
                       {\frac{1}{\sqrt{8\pi}}{\exp{-x^2/8}}}
     = \frac{p(x)}{\pi(x)}
$$

is called the weight, $p$ is the pdf from which we wish to sample and
$\pi$ is the pdf of the importance distribution.

Importance Sampling Approximation of the Posterior
==================================================

Suppose that the posterior distribution of a model in which we are
interested has a complicated functional form and that we therefore
wish to approximate it in some way. First assume that we wish to
calculate the expectation of some arbitrary function $f$ of the
parameters.

$$
{\mathbb{E}}(f({x}) \,\vert\, y_1, \ldots y_T) =
\int_\Omega f({x}) p({x} \, \vert \, y_1, \ldots y_T) \,\mathrm{d}{x}
$$

Using Bayes

$$
\int_\Omega f({x}) \condprob{p}{x}{y_1, \ldots y_T} \,\mathrm{d}{x} =
\frac{1}{Z}\int_\Omega f({x}) \condprob{p}{y_1, \ldots y_T}{x}p(x) \,\mathrm{d}{x}
$$

where $Z$ is some normalizing constant.

As before we can re-write this using a proposal distribution $\pi(x)$

$$
\frac{1}{Z}\int_\Omega f({x}) \condprob{p}{y_1, \ldots y_T}{x}p(x) \,\mathrm{d}{x} =
\frac{1}{Z}\int_\Omega \frac{f({x}) \condprob{p}{y_1, \ldots y_T}{x}p(x)}{\pi(x)}\pi(x) \,\mathrm{d}{x}
$$

We can now sample $X^{(i)} \sim \pi({x})$ repeatedly to obtain

$$
{\mathbb{E}}(f({x}) \,\vert\, y_1, \ldots y_T) \approx \frac{1}{ZN}\sum_1^N
f({X^{(i)}}) \frac{p(y_1, \ldots y_T \, \vert \, {X^{(i)}})p({X^{(i)}})}
                            {\pi({X^{(i)}})} =
\sum_1^N w_if({X^{(i)}})
$$

where the weights $w_i$ are defined as before by

$$
w_i = \frac{1}{ZN} \frac{p(y_1, \ldots y_T \, \vert \, {X^{(i)}})p({X^{(i)}})}
                        {\pi({X^{(i)}})}
$$

We follow [Alex
Cook](http://blog.nus.edu.sg/alexcook/teaching/sph6004/) and use the
example from [@citeulike:5986027]. We take the prior as $\sim
{\cal{Be}}(1,1)$ and use ${\cal{U}}(0.0,1.0)$ as the proposal
distribution. In this case the proposal and the prior are identical
just expressed differently and therefore cancel.

Note that we use the log of the pdf in our calculations otherwise we
suffer from (silent) underflow, e.g.,

    [ghci]
    pdf (Binomial nv (0.4 :: Double)) xv

On the other hand if we use the log pdf form

    [ghci]
    logPdf (Binomial nv (0.4 :: Double)) xv

> xv, nv :: Int
> xv = 51
> nv = 8197

> sampleUniform :: RVarT (W.Writer [Double]) ()
> sampleUniform = do
>   x <- rvarT StdUniform
>   lift $ W.tell [x]
>   return ()

> runSampler :: RVarT (W.Writer [Double]) () ->
>               Int -> Int -> [Double]
> runSampler sampler seed n =
>   snd $
>   W.runWriter $
>   evalStateT (sample (replicateM n sampler))
>              (pureMT (fromIntegral seed))

> sampleSize :: Int
> sampleSize = 1000

> pv :: [Double]
> pv = runSampler sampleUniform 2 sampleSize

> logWeightsRaw :: [Double]
> logWeightsRaw = map (\p -> logPdf (Beta 1.0 1.0) p +
>                            logPdf (Binomial nv p) xv -
>                            logPdf StdUniform p) pv

> logWeightsMax :: Double
> logWeightsMax = maximum logWeightsRaw
>
> weightsRaw :: [Double]
> weightsRaw = map (\w -> exp (w - logWeightsMax)) logWeightsRaw

> weightsSum :: Double
> weightsSum = sum weightsRaw

> weights :: [Double]
> weights = map (/ weightsSum) weightsRaw

> meanPv :: Double
> meanPv = sum $ zipWith (*) pv weights
>
> meanPv2 :: Double
> meanPv2 = sum $ zipWith (\p w -> p * p * w) pv weights
>
> varPv :: Double
> varPv = meanPv2 - meanPv * meanPv

We get the answer

    [ghci]
    meanPv

But if we look at the size of the weights and the effective sample size

    [ghci]
    length $ filter (>= 1e-6) weights
    (sum weights)^2 / (sum $ map (^2) weights)

so we may not be getting a very good estimate. Let's try

> sampleNormal :: RVarT (W.Writer [Double]) ()
> sampleNormal = do
>   x <- rvarT $ Normal meanPv (sqrt varPv)
>   lift $ W.tell [x]
>   return ()

> pvC :: [Double]
> pvC = runSampler sampleNormal 3 sampleSize

> logWeightsRawC :: [Double]
> logWeightsRawC = map (\p -> logPdf (Beta 1.0 1.0) p +
>                             logPdf (Binomial nv p) xv -
>                             logPdf (Normal meanPv (sqrt varPv)) p) pvC

> logWeightsMaxC :: Double
> logWeightsMaxC = maximum logWeightsRawC
>
> weightsRawC :: [Double]
> weightsRawC = map (\w -> exp (w - logWeightsMaxC)) logWeightsRawC

> weightsSumC :: Double
> weightsSumC = sum weightsRawC

> weightsC :: [Double]
> weightsC = map (/ weightsSumC) weightsRawC

> meanPvC :: Double
> meanPvC = sum $ zipWith (*) pvC weightsC

> meanPvC2 :: Double
> meanPvC2 = sum $ zipWith (\p w -> p * p * w) pvC weightsC
>
> varPvC :: Double
> varPvC = meanPvC2 - meanPvC * meanPvC

Now the weights and the effective size are more re-assuring

    [ghci]
    length $ filter (>= 1e-6) weightsC
    (sum weightsC)^2 / (sum $ map (^2) weightsC)

And we can take more confidence in the estimate

    [ghci]
    meanPvC

Bibliography
============