# Bayesian updating using $n$ noisy observations of Brownian motion

I am very new to Bayesian inference and can't figure out what may be an elementary problem. Also, please forgive me if I am screwing up the notation -- this is my first foray into Bayesian statistics.

Set-up: At time $i=0$ I start with a random variable $X_0\sim \mathcal{N}(\mu_0,\sigma^2_0)$. Over time I observe the Brownian motion process such that for $j>i$, $X_j-X_i\sim \mathcal{N}(0,(j-i)D)$. Starting at time 0, I collect a sequence of $n$ observations $\{Y_i\}_{i=0}^{n-1}$ of $\{X_i\}_{i=0}^{n-1}$, which are subject to Gaussian noise such that $p(y_i|x_i)=\frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(y_i-x_i)^2}{2N_0}}$. Noise is independent from observation to observation. I am interested in the mean and variance (of the estimator) of $X_{n-1}$ given $\{Y_i\}_{i=0}^{n-1}$ in terms of $n$, $\mu_0$, $\sigma^2_0$, $D$, and $N_0$.

What I've done: Using Bayes rule and the properties of Gaussian distribution, I've found that for the first observation ($n=0$):

$$p(x_0|y_0)=\frac{1}{\sqrt{2\pi S}}e^{-\frac{(x_0-M)^2}{2S}}$$

where $M=\frac{y_0\sigma^2_0+\mu_0N_0}{\sigma^2_0+N_0}$ and $S=\frac{\sigma^2_0N_0}{\sigma^2_0+N_0}$. I am not completely sure if the above is correct, but I am pretty confident that it is, as it looks like the mean squared error of the estimate of $X_0$ has decreased with the observation. However, I am having trouble extending this to $n>0$ and would appreciate any guidance.

EDIT: I've misspecified the Brownian motion process. I think it's correct now.

-
Note that your distribution for $X_i$ is conditional on only knowing $X_{0}$. if you know say $X_{i-s}$ then the distribution is $X_{i}|X_{i-s}\sim N(X_{i-s},(i-s)\sigma^2)$. although i'm not sure if they mean parameter is "drift" or the starting point. If its drift then you need to add $(i-s)\mu_{0}$ to the mean. Interesting question though. –  probabilityislogic Feb 15 '12 at 0:29
Remember brownian motion is about independent increments, and is usually best described that way. –  probabilityislogic Feb 15 '12 at 0:31
@probabilityislogic I made a mistake in specifying the Brownian motion process. I think it's fixed now. I'm still lost though. I think I need to use the updated mean and variance as a prior for $n=1$ and proceed recursively, but I'm not sure how to do that... –  M.B.M. Feb 15 '12 at 3:59

So $X_{n-1}$ is the target of inference. Then the correct Bayesian procedure is to calculate this posterior:

$$p(X_{n-1}|Y_0\dots Y_{n-1}\mu_0\sigma_0N_0DI)$$

(Note I have added the symbol $I$ to explicitly indicate the assumptions and model structure being used). Now because we also know that the other $X_i$ are important, but not present in the above posterior, then they must have been integrated out (via the law of total probability):

$$\int \dots \int p(X_{n-1}X_{n-2}\dots X_0|Y_0\dots Y_{n-1}\mu_0\sigma_0N_0DI)dX_{0}\dots dX_{n-2}$$

Now I will suppress the $\mu_0\sigma_0N_0DI$ from the probabilities for brevity, but they are still there hiding. Now the above is not directly calculable using the information given, but we can use the rules of probability theory to manipulate this into something we do know how to calculate. We can use bayes theorem to get

$$p(X_{n-1}|Y_0\dots Y_{n-1})$$ $$=\frac{\int \dots \int p(X_{n-1}\dots X_0)p(Y_{n-1}\dots Y_0|X_{n-1}\dots X_0)dX_{0}\dots dX_{n-2}}{\int \dots \int p(X_{n-1}\dots X_0)p(Y_{n-1}\dots Y_0|X_{n-1}\dots X_0) dX_{0}\dots dX_{n-1}}$$

Because we have a ratio of integrals with the same integrand, we can eliminate constants which don't depend on $X_0,\dots,X_{n-1}$. This makes the equations slightly smaller and easier to write. Now we can use the "markovian" property of brownian motion to simplify the joint prior.

$$p(X_{n-1}\dots X_0)=p(X_{n-1}|X_{n-2})p(X_{n-2}|X_{n-3})\dots p(X_1|X_0)p(X_0)$$

Now for brownian motion we have $p(X_i|X_{i-1})\sim N(X_{i-1},D)$ (remembering that conditional on $X_{i-1}$ means it is just a constant, not a random variable) and combined with your prior for $X_0$ we have:

$$p(X_{n-1}\dots X_0)\propto\exp\left(-\frac{1}{2D}\sum_{i=1}^{n-1}(X_i-X_{i-1})^2-\frac{(X_0-\mu_0)^2}{2\sigma_0^2}\right)$$

Now for the likelihood I will make an assumption which I think you have implicitly made. This is that the $Y_i$ are conditionally independent, given $X_i$. I think this is what you mean by "noise is independent from observation to observation". This allows us to factor the likelihood as:

$$p(Y_{n-1}\dots Y_0|X_{n-1}\dots X_0)=\prod_{i=0}^{n-1}p(Y_{i}|X_{i})\propto\exp\left(-\frac{1}{2N_0}\sum_{i=0}^{n-1}(Y_i-X_i)^2\right)$$

Hence the "integrand" is given by the product of the prior and the likelihood, stripped of constants which don't depend on $X_i$. This is given by: $$p(X_{n-1}\dots X_0)p(Y_{n-1}\dots Y_0|X_{n-1}\dots X_0)$$ $$\propto\exp\left(-\frac{1}{2N_0}\sum_{i=0}^{n-1}(Y_i-X_i)^2-\frac{1}{2D}\sum_{i=1}^{n-1}(X_i-X_{i-1})^2-\frac{1}{2\sigma_0^2}(X_0-\mu_0)^2\right)$$ $$\propto\exp\left(-\frac{X_0^2-2X_0\mu_0}{2\sigma_0^2}-\sum_{i=0}^{n-1}\frac{X_i^2-2X_iY_i}{2N_0}-\sum_{i=1}^{n-1}\frac{(X_i-X_{i-1})^2}{2D}\right)$$

Now we "simply" integrate this over $X_0,\dots,X_{n-2}$, noting that we only need to keep the terms which depend in $X_0,\dots,X_{n-1}$. This will be a straight-forward but tedious process if we were to directly integrate out from this equation. An easier way to obtain the integral is to evaluate the posterior for $n=2$ and $n=3$ cases which will establish a recursive relationship between the posterior $p(X_t|Y_0,\dots,Y_t)$ and $p(X_{t+1}|Y_0,\dots,Y_{t+1})$

$$\bf{}n=2\text{ case}$$

The integrand we require is given by:

$$p(X_{1}X_0)p(Y_{1}Y_0|X_{1}X_0)$$ $$\propto\exp\left(-\frac{X_0^2-2X_0\mu_0}{2\sigma_0^2}-\frac{X_1^2+X_0^2-2X_1Y_1-2X_0Y_0}{2N_0}-\frac{(X_1-X_0)^2}{2D}\right)$$ $$=\exp\left(-a_0X_0^2+b_0X_0-\frac{X_1^2-2X_1Y_1}{2N_0}-\frac{X_1^2}{2D}\right)$$

Where $a_0=\frac{1}{2\sigma_0^2}+\frac{1}{2N_0}+\frac{1}{2D}=\frac{1}{2S}+\frac{1}{2D}$ and $b_0=\frac{\mu_0}{\sigma_0^2}+\frac{Y_0}{N_0}+\frac{X_1}{D}=\frac{M}{S}+\frac{X_1}{D}$. I have used $M$ and $S$ from your definition in the question. Now because $-a_0X_0^2+b_0X_0=-a_0\left(X_0-\frac{b_0}{2a_0}\right)^2+\frac{b_0^2}{4a_0}$ we can consider $\hat{X}_0=\frac{b_0}{2a_0}$ as an estimate for $X_0$ with variance equal to $v_0=\frac{1}{2a_0}$. Using these relations we have $\frac{b_0^2}{4a_0}=\frac{\hat{X}_0^2}{2v_0}$. Note that it has the same weighted average form as your $M$ for $n=1$ case. Using this we have a simple gaussian integral, which basically replaces $X_0$ by its estimate. To see this note that $-a_0\hat{X}_0^2+b_0\hat{X}_0=\frac{b_0^2}{4a_0}$, and we have:

$$\int p(X_{1}X_0)p(Y_{1}Y_0|X_{1}X_0)dX_0\propto\sqrt{\frac{\pi}{a_0}}\exp\left(\frac{b_0^2}{4a_0}-\frac{X_1^2-2X_1Y_1}{2N_0}-\frac{X_1^2}{2D}\right)$$ $$\propto\exp\left(-a_1X_1^2+b_1X_1\right)$$

where $a_1=\frac{1}{2N_0}+\frac{1}{2(D+S)}$ and $b_1=\frac{M}{D+S}+\frac{Y_1}{N_0}$. As before, this gives rise to the estimate $M_1=\frac{b_1}{2a_1}$ with variance $S_1=\frac{1}{2a_1}$. Additionally the form for the above expression is a gaussian kernel (i.e. a gaussian pdf stripped of normalising constants). So the posterior for $n=2$ is given by:

$$p(X_1|Y_0Y_1\mu_0\sigma_0N_0DI)=\frac{1}{\sqrt{2\pi S_1}}\exp\left(-\frac{1}{2S_1}\left[X_1-M_1\right]^2\right)$$

Where $M_1=\frac{MN_0+Y_1(D+S)}{N_0+D+S}$ is the estimate and $S_1=\frac{N_0(D+S)}{N_0+D+S}$. Notice here that this is the same form as for $n=1$ except with $M$ as the new "prior estimate" and $S+D$ as the new "prior variance". Note that $D$ in the prior variance accounts for the uncertainty in the inference path $X_0\to X_1$ and $S$ accounts for uncertainty in the inference path $\mu_0,Y_0\to X_0$.

$$\bf{}n=3\text{ case}$$

The integration over $X_0$ is the same, however the integration over $X_1$ is different, because now $X_2$ can be used to help estimate $X_1$ - this is due to the additional term $\frac{(X_2-X_1)^2}{2D}$. So we now have (note the similar form to the $n=2$ case):

$$\int p(X_2X_1X_0)p(Y_2Y_1Y_0|X_2X_1X_0)dX_0$$ $$\propto\exp\left(\frac{b_0^2}{4a_0}-\frac{X_1^2-2X_1Y_1}{2N_0}-\frac{X_1^2}{2D} -\frac{X_2^2-2X_2Y_2}{2N_0}-\frac{(X_2-X_1)^2}{2D}\right)$$ $$\propto\exp\left(-a_1^{'}X_1^2+b_1^{'}X_1-\frac{X_2^2-2X_2Y_2}{2N_0}-\frac{X_2^2}{2D}\right)$$

Where $a_1^{'}=\frac{1}{2S_1}+\frac{1}{2D}$ and $b_1^{'}=\frac{M_1}{S_1}+\frac{X_2}{D}$ We are now ready to integrate out $X_1$ and we get:

$$\int p(X_2X_1X_0)p(Y_2Y_1Y_0|X_2X_1X_0)dX_0dX_1$$ $$\propto\exp\left(\frac{[b_1^{'}]^2}{4a_1^{'}}-\frac{X_2^2-2X_2Y_2}{2N_0}-\frac{X_2^2}{2D}\right)\propto\exp\left(-a_2X_2^2+b_2X_2\right)$$

Where $a_2=\frac{1}{2N_0}+\frac{1}{2(D+S_1)}$ and $b_2=\frac{M_1}{S_1+D}+\frac{Y_2}{N_0}$. As before these give a gaussian kernel:

$$p(X_2|Y_0Y_1Y_2\mu_0\sigma_0N_0DI)=\frac{1}{\sqrt{2\pi S_2}}\exp\left(-\frac{1}{2S_2}\left[X_2-M_2\right]^2\right)$$

Where $M_2=\frac{M_1N_0+Y_2(D+S_1)}{N_0+D+S_1}$ is the estimate and $S_2=\frac{N_0(D+S_1)}{N_0+D+S_1}$. Notice here that this is the same form as for $n=2$ except with $M_1$ as the new "prior estimate" and $S_1+D$ as the new "prior variance". Note that $D$ in the prior variance accounts for the uncertainty in the inference path $X_1\to X_2$ and $S_1$ accounts for uncertainty in the inference path $\mu_0,Y_0,Y_1\to X_1$.

$$\bf{}\text{General case}$$

this has an obvious general case via mathematical induction on the previous two cases. The result can be stated recursively. Set $M_0=\frac{Y_0\sigma_0^2+\mu_0N_0}{\sigma^2_0+N_0}$ and $S_0=\frac{\sigma^2_0N_0}{\sigma^2_0+N_0}$. Then, given the posterior mean and variance at point after observing $Y_0,\dots,Y_t$ are $M_t$ and $S_t$ respectively the posterior mean and variance given $Y_0,\dots,Y_{t+1}$ is given by:

$$M_{t+1}=\frac{M_tN_0+Y_{t+1}(D+S_t)}{N_0+D+S_t}$$ $$S_{t+1}=\frac{N_0(D+S_t)}{N_0+D+S_t}$$

All posterior distributions are normal. Hence for the most recent part of the chain the posterior mean is given by:

$$M_{n-1}=\frac{M_{n-2}N_0+Y_{n-1}(D+S_{n-2})}{N_0+D+S_{n-2}}$$ $$S_{n-1}=\frac{N_0(D+S_{n-2})}{N_0+D+S_{n-2}}$$

You would then recursively "build up" this estimate starting from $M_0$ and $S_0$

-
+1 That is one of the longest algorithmic answers I've ever seen. I almost want to frame it. –  Xeoncross Feb 17 '12 at 19:38
Wow, thank you! Yes, the conditional independence assumption for observation is what I had in mind. I've been gradually building my knowledge about this stuff, but I only got to the $n=2$ case. Was going to do $n=3$ and see if I could see the relationship. I'll read your answer very carefully (and compare to my what I've done) this weekend and accept it then. Again, thanks! –  M.B.M. Feb 17 '12 at 21:21
@bullmose - regarding your comment about throwing away constants my general approach is to do it unless you've anticipated what the answer will be. I tend to make less mistakes this way. and thanks for your compliments. –  probabilityislogic Feb 18 '12 at 8:43
@bullmose - Regarding additonal observations increasing uncertainty we can look at the precision $\frac{1}{S_{t+1}}=\frac{1}{N_{0}}+\frac{1}{D+S_{t}}$. For ease of example i will use $\sigma_0^2>>N_{0}>>D$ so we have $\frac{1}{S_{0}}\approx\frac{1}{N_{0}}$ this means that $\frac{1}{S_{1}}\approx\frac{1}{N_{0}}+\frac{1}{D+N_{0}}\approx\frac_{2}{N_{0}}$‌​. this is a doubling in precision. For next observation we get $\frac{1}{S_{2}}\approx\frac{1}{N_{0}}+\frac{1}{D+N_{0}/2}\approx\frac_{3}{N_{0}‌​}$. hence each observation increases precision by $\frac_{1}{N_{0}}$. –  probabilityislogic Feb 18 '12 at 8:46
@bullmose - alternatively we can use $\sigma_0^2>>D>>N_{0}$ so we have $\frac{1}{S_{0}}\approx\frac{1}{N_{0}}$ this means that $\frac{1}{S_{1}}\approx\frac{1}{N_{0}}+\frac{1}{D+N_{0}}\approx\frac_{1}{N_{0}}+‌​\frac{1}{D}\approx\frac{1}{N_{0}}$. this is a negligible increase in precision. For next observation we get $\frac{1}{S_{2}}\approx\frac{1}{N_{0}}+\frac{1}{D+N_{0}}\approx\frac_{1}{N_{0}}$‌​. this is a negligible increase in precision –  probabilityislogic Feb 18 '12 at 8:58