Filtering.lhs

% Fun with Filters
% Dominic Steinitz
% 3rd July 2014

---
bibliography: Kalman.bib
---

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

Noisy Observation
=================

Suppose we wish to estimate the mean of a sample drawn from a normal
distribution. In the Bayesian approach, we know the prior distribution
for the mean (it could be a non-informative prior) and then we update
this with our observations to create the posterior, the latter giving
us improved information about the distribution of the mean. In symbols

$$
p(\theta \,\vert\, x) \propto p(x \,\vert\, \theta)p(\theta)
$$

Typically, the samples are chosen to be independent, and all of the
data is used to perform the update but, given independence, there is
no particular reason to do that, updates can performed one at a time
and the result is the same; nor is the order of update
important. Being a bit imprecise, we have

$$
p(z \,\vert\, x, y) = p(z, x, y)p(x, y) = p(z, x, y)p(x)p(y) =
p((z \,\vert\, x) \,\vert\, y) =
p((z \,\vert\, y) \,\vert\, x)
$$

The standard notation in Bayesian statistics is to denote the
parameters of interest as $\theta \in \mathbb{R}^p$ and the
observations as $x \in \mathbb{R}^n$. For reasons that will become
apparent, let us change notation and label the parameters as $x$ and
the observations as $y$.


Let us take a very simple example of a prior $X \sim {\cal{N}}(0,
\sigma^2)$ where $\sigma^2$ is known and then sample from a normal
distribution with mean $x$ and variance for the $i$-th sample $c_i^2$
where $c_i$ is known (normally we would not know the variance but
adding this generality would only clutter the exposition
unnecessarily).

$$
p(y_i \,\vert\, x) = \frac{1}{\sqrt{2\pi c_i^2}}\exp\bigg(\frac{(y_i - x)^2}{2c_i^2}\bigg)
$$

The likelihood is then

$$
p(\boldsymbol{y} \,\vert\, x) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi c_i^2}}\exp\bigg(\frac{(y_i - x)^2}{2c_i^2}\bigg)
$$

As we have already noted, instead of using this with the prior to
calculate the posterior, we can update the prior with each observation
separately. Suppose that we have obtained the posterior given $i - 1$
samples (we do not know this is normally distributed yet but we soon
will):

$$
p(x \,\vert\, y_1,\ldots,y_{i-1}) = {\cal{N}}(\hat{x}_{i-1}, \hat{\sigma}^2_{i-1})
$$

Then we have

$$
\begin{aligned}
p(x \,\vert\, y_1,\ldots,y_{i}) &\propto p(y_i \,\vert\, x)p(x \,\vert\, y_1,\ldots,y_{i-1}) \\
&\propto \exp-\bigg(\frac{(y_i - x)^2}{2c_i^2}\bigg) \exp-\bigg(\frac{(x - \hat{x}_{i-1})^2}{2\hat{\sigma}_{i-1}^2}\bigg) \\
&\propto \exp-\Bigg(\frac{x^2}{c_i^2} - \frac{2xy_i}{c_i^2} + \frac{x^2}{\hat{\sigma}_{i-1}^2} - \frac{2x\hat{x}_{i-1}}{\hat{\sigma}_{i-1}^2}\Bigg) \\
&\propto \exp-\Bigg( x^2\Bigg(\frac{1}{c_i^2} + \frac{1}{\hat{\sigma}_{i-1}^2}\Bigg) - 2x\Bigg(\frac{y_i}{c_i^2} + \frac{\hat{x}_{i-1}}{\hat{\sigma}_{i-1}^2}\Bigg)\Bigg)
\end{aligned}
$$

Writing

$$
\frac{1}{\hat{\sigma}_{i}^2} \triangleq \frac{1}{c_i^2} + \frac{1}{\hat{\sigma}_{i-1}^2}
$$

and then completing the square we also obtain

$$
\frac{\hat{x}_{i}}{\hat{\sigma}_{i}^2} \triangleq \frac{y_i}{c_i^2} + \frac{\hat{x}_{i-1}}{\hat{\sigma}_{i-1}^2}
$$

More Formally
-------------

Now let's be a bit more formal about conditional probability and use
the notation of $\sigma$-algebras to define ${\cal{F}}_i =
\sigma\{Y_1,\ldots, Y_i\}$ and $M_i \triangleq \mathbb{E}(X \,\vert\,
{\cal{F}}_i)$ where $Y_i = X + \epsilon_i$, $X$ is as before and
$\epsilon_i \sim {\cal{N}}(0, c_k^2)$. We have previously calculated
that $M_i = \hat{x}_i$ and that ${\cal{E}}((X - M_i)^2 \,\vert\, Y_1,
\ldots Y_i) = \hat{\sigma}_{i}^2$ and the tower law for conditional
probabilities then allows us to conclude ${\cal{E}}((X - M_i)^2) =
\hat{\sigma}_{i}^2$. By [Jensen's
inequality](http://en.wikipedia.org/wiki/Jensen%27s_inequality), we have

$$
{\cal{E}}(M_i^2) = {\cal{E}}({\cal{E}}(X \,\vert\, {\cal{F}}_i)^2)) \leq
{\cal{E}}({\cal{E}}(X^2 \,\vert\, {\cal{F}}_i))) =
{\cal{E}}(X^2) = \sigma^2
$$

Hence $M$ is bounded in $L^2$ and therefore converges in $L^2$ and
almost surely to $M_\infty \triangleq {\cal{E}}(X \,\vert\,
{\cal{F}}_\infty)$.

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         #-}

> {-# LANGUAGE FlexibleContexts             #-}

> module Filtering where

> import Control.Monad
> import Data.Random.Source.PureMT
> import Data.Random
> import Data.Random.Distribution.Binomial
> import Data.Random.Distribution.Multinomial
> 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.

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

As before we can re-write this

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

We can now sample $X^{(i)} \sim \pi({x} \, \vert \,
y_1, \ldots y_T)$ repeatedly to obtain

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

where the weights $w_i$ are defined as before by

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

We follow [Alex
Cook](http://blog.nus.edu.sg/alexcook/teaching/sph6004/) and use the
example from [@citeulike:5986027].

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

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

> sampleSize :: Int
> sampleSize = 10000

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

> weightsRaw :: [Double]
> weightsRaw = map (\p -> pdf (Binomial nv p) xv) pv

> weightsSum :: Double
> weightsSum = sum weightsRaw

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

> meanPv :: Double
> meanPv = sum $ zipWith (*) pv weights

But

    [ghci]
    length $ filter (>= 1e-6) weights

so we may not be getting a very good estimate.

> sampleNormal :: RVarT (W.Writer [Double]) ()
> sampleNormal = do
>   let xvd   = fromIntegral xv
>       nvd   = fromIntegral nv
>       mu    = xvd / nvd
>       sigma = (sqrt (xvd / nvd) * (1 - xvd / nvd)) / nvd
>   x <- rvarT $ Normal mu sigma
>   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))

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

> weightsRawC :: [Double]
> weightsRawC = map (\p -> pdf (Binomial nv p) xv) pvC

> weightsSumC :: Double
> weightsSumC = sum weightsRawC

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

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

    [ghci]
    length $ filter (>= 1e-6) weightsC

Sequential Importance Sampling
==============================

Now let us generalize our Bayesian model and allow
${X}_i$ to depend on ${X}_{i-1}$ and
${Y}_i$ depend on ${X}_i$, that is

$$
\begin{aligned}
{X}_i &\sim \condprob{p}{{x}_i}{{x}_{i-1}} \\
{Y}_i &\sim \condprob{p}{{y}_i}{{x}_{i}} \\
\end{aligned}
$$

Recall in our original model we had

$$
\begin{aligned}
{X}_i &\sim p({{x}}) \\
{Y}_i &\sim \condprob{p}{{y}_i}{{x}} \\
\end{aligned}
$$

And, as before, by Bayes and conditional independence we also had

$$
\begin{aligned}
\condprob{p}{{x}}{{y}_1, \ldots, {y}_T} & \propto
p({x})\condprob{p}{{y}_1, \ldots, {y}_T}{{x}} \\
&  =p({x})\prod_{i=1}^T\condprob{p}{{y}_i}{{x}} \\
& = p({x})\condprob{p}{{y}_T}{{x}}\prod_{i=1}^{T-1}\condprob{p}{{y}_i}{{x}} \\
& \propto \condprob{p}{{y}_T}{{x}}\condprob{p}{{x}}{{y}_1, \ldots, {y}_{T-1}}
\end{aligned}
$$

In our new and generalized model we have a whole sequence of
parameters for which we wish to find the posterior distribution.

$$
\condprob{p}{{x}_0,\ldots,{x}_T}{{y}_1,\ldots,{y}_T}
$$

Applying Bayes and the restrictions of our model that the parameters are Markovian

$$
\begin{aligned}
\condprob{p}{{x}_0,\ldots,{x}_T}{{y}_1,\ldots,{y}_T}
& \propto
\condprob{p}{{y}_T}{{x}_0,\ldots,{x}_T, {y}_1,\ldots,{y}_{T-1}}
\condprob{p}{{x}_0,\ldots,{x}_T}{{y}_1,\ldots,{y}_{T-1}} \\
& =
\condprob{p}{{y}_T}{{x}_T}
\condprob{p}{{x}_T}{{x}_0,\ldots,{x}_{T-1}, {y}_1,\ldots,{y}_{T-1}}\condprob{p}{{x}_0,\ldots,{x}_{T-1}}{{y}_1,\ldots,{y}_{T-1}} \\
& =
\condprob{p}{{y}_T}{{x}_T}
\condprob{p}{{x}_T}{{x}_{T-1}}\condprob{p}{{x}_0,\ldots,{x}_{T-1}}{{y}_1,\ldots,{y}_{T-1}} \\
\end{aligned}
$$

Now let us take an importance distribution
$\condprob{\pi}{x_0,\ldots,x_T}{y_1,\ldots,y_T}$ and as before compute
the un-normalized weights by substituting in the posterior we have
just calculated.

$$
w_T^{(i)} = \frac
{\condprob{p}{{y}_T}{{X_T^{(i)}}}
 \condprob{p}{{X_T^{(i)}}}{{X_{T-1}^{(i)}}}\condprob{p}{{X_0^{(i)}},\ldots,{X_{T-1}^{(i)}}}{{y}_1,\ldots,{y}_{T-1}}
}
{\condprob{\pi}{X^{(i)}_0,\ldots,X^{(i)}_T}{y_1,\ldots,y_T}}
$$

The problem now is how to select an appropriate importance distribution.
Suppose we can do this recursively.

$$
{\condprob{\pi}{x_0,\ldots,x_T}{y_1,\ldots,y_T}} =
{\condprob{\pi}{x_T}{x_0,\ldots,x_{T-1}, y_1,\ldots,y_T}}
{\condprob{\pi}{x_0,\ldots,x_{T-1}}{y_1,\ldots,y_{T-1}}}
$$

then

$$
\begin{aligned}
w_T^{(i)} & = \frac
{\condprob{p}{{y}_T}{{X_T^{(i)}}}
 \condprob{p}{{X_T^{(i)}}}{{X_{T-1}^{(i)}}}
}
{
{\condprob{\pi}{X_T^{(i)}}{X_0^{(i)},\ldots,X_{T-1}^{(i)}, y_1,\ldots,y_T}}
}
\frac
\condprob{p}{{X_0^{(i)}},\ldots,{X_{T-1}^{(i)}}}{{y}_1,\ldots,{y}_{T-1}}
{\condprob{\pi}{X_0^{(i)},\ldots,X_{T-1}^{(i)}}{y_1,\ldots,y_{T-1}}} \\
& =
\frac
{\condprob{p}{{y}_T}{{X_T^{(i)}}}
 \condprob{p}{{X_T^{(i)}}}{{X_{T-1}^{(i)}}}
}
{
{\condprob{\pi}{X_T^{(i)}}{X_0^{(i)},\ldots,X_{T-1}^{(i)}, y_1,\ldots,y_T}}
}w^{(i)}_{T-1}
\end{aligned}
$$

We can start the whole sampling process by sampling $X_0^{(i)} \sim
p(x_0)$ from the prior and setting $w_0^{(i)} = 1/N$ where $N$ is the
total number of samples.

Further, if we choose ${\condprob{\pi}{X_T}{X_0,\ldots,X_{T-1},
y_1,\ldots,y_T}} = {\condprob{\pi}{X_T}{X_{T-1}, y_1,\ldots,y_T}}$ then we
do not need to keep the whole history of parameters
$X_0^{(i)},\ldots,X_{T-1}^{(i)}$ only the previous parameters.

This simplifies the recursive weight equations

$$
\begin{aligned}
w_T^{(i)}
& =
\frac
{\condprob{p}{{y}_T}{{X_T^{(i)}}}
 \condprob{p}{{X_T^{(i)}}}{{X_{T-1}^{(i)}}}
}
{
{\condprob{\pi}{X_T^{(i)}}{X_{T-1}^{(i)}, y_1,\ldots,y_T}}
}w^{(i)}_{T-1}
\end{aligned}
$$

If we further take $\condprob{\pi}{X_T^{(i)}}{X_{T-1}^{(i)},
y_1,\ldots,y_T}$ to be $\condprob{p}{{X_T^{(i)}}}{{X_{T-1}^{(i)}}}$
weighted by $w_{T-1}^{(i)}$ then

$$
\begin{aligned}
w_T^{(i)}
& =
\condprob{p}{{y}_T}{{\tilde{X}_T^{(i)}}}
\end{aligned}
$$

where

$$
\tilde{X}_T \sim \sum_{i=1}^N w_{T-1}^{(i)} \condprob{p}{X_T}{X_{T-1}^{(i)}}
$$

This second sampling according to the weights is known as **resampling**.

> muPrior, sigmaPrior, muLikelihood :: Double
> muPrior = 0.0
> sigmaPrior = 1.0
> muLikelihood = 0.0

> cs :: [Double]
> cs = repeat 1.0

> bigN :: Int
> bigN = 100

> normalPdf :: Double -> Double -> Double -> Double
> normalPdf mu sigma x =
>   (recip (sqrt (2 * pi * sigma2))) * (exp ((-(x - mu)^2) / (2 * sigma2)))
>   where
>     sigma2 = sigma^2

> sir :: [(Double, Double)] -> Double -> Double ->
>        RVarT (W.Writer [Double]) [(Double, Double)]
> sir weightsMusPrev y sigma = do
>   let n           = length weightsMusPrev
>       weightsPrev = map fst weightsMusPrev
>
>   let musPrev  = map snd weightsMusPrev
>   let m = sum $ zipWith (*) weightsPrev musPrev
>
>   nParticless <- rvarT $ Multinomial weightsPrev n
>
>   let musTilde = concatMap (\i -> replicate (nParticless!!i) (musPrev!!i)) [0..n - 1]
>
>   musNew <- return musTilde
>
>   let weightsNew = map (\i -> normalPdf (musNew!!i) sigma y) [0..n - 1]
>       sumWeights = sum weightsNew
>       weightsNorm = map (/ sumWeights) weightsNew
>
>   lift $ W.tell [m]
>   return (zip weightsNorm musNew)

> initSir :: Int -> RVarT (W.Writer [Double]) [(Double, Double)]
> initSir n = do
>   mus <- replicateM n (rvarT $ Normal muPrior sigmaPrior)
>   let weights = replicate n (recip (fromIntegral n))
>   let m = sum $ zipWith (*) weights mus
>   lift $ W.tell [m]
>   return (zip weights mus)

> createObs :: Int -> RVar (Double, [Double])
> createObs n = do
>   x <- rvarT (Normal muPrior sigmaPrior)
>   ys <- mapM (\c -> rvarT (Normal x c)) (take n cs)
>   return (x, ys)

> obss :: (Double, [Double])
> obss = evalState (sample (createObs 100)) (pureMT 2)

> testPF :: RVarT (W.Writer [Double]) [(Double, Double)]
> testPF = do
>   wms <- initSir bigN
>   let sirs = zipWith (\obs c -> (\wm -> sir wm obs c)) (snd obss) cs
>   foldr (>=>) return sirs wms

> runPF :: ([(Double, Double)], [Double])
> runPF = W.runWriter (evalStateT (sample testPF) (pureMT 2))

```{.dia width='900'}
dia = image "diagrams/SingleRvNoisy.png" 1.0 1.0
```

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