A Kalman Question

[Dominic Steinitz]

I have samples from a normal distribution. If I run a Kalman filter then the estimate of the mean should be the cumulative mean of the data. We can try this out:

> model MV {
> param sigma;
> state mu;
> obs x;
> noise eta;
>
> sub parameter {
> sigma ~ uniform(0.0, 1.0);
> }
>
> sub proposal_parameter {
> sigma ~ truncated_gaussian(sigma, 0.1, 0.0, 1.0);
> }
>
> sub initial {
> mu ~ normal(0.0, 1.0);
> }
>
> sub transition {
> eta ~ normal(0.0, 1.0e-6);
> mu <- mu + eta;
> }
>
> sub observation {
> x ~ normal(mu, sigma);
> }
> }

And then run this

> rm(list = ls(all.names=TRUE))
> unlink(".RData")
>
> library('rbi')
> library('ggplot2')
>
> T <- 200
>
> MV <- bi_model("MeanAndVar.bi")
> synthetic_dataset <- bi_generate_dataset(end_time=T,
> model=MV,
> seed=43,
> init=list(mu=0.5, sigma=0.1, x=0.5))
>
> rdata <- bi_read(synthetic_dataset)
>
> MVInferKalman <- bi_model("MeanAndVar.bi")
> bi_kalman_object <- libbi(client="filter", model=MVInferKalman, obs=synthetic_dataset)
>
> bi_kalman_object$run(add_options = list(
> "end-time" = T,
> noutputs = T),
> init=list(mu=0.5, sigma=0.1, x=0.5),
> filter="kalman",
> stdoutput_file_name = tempfile(pattern="kalmanoutput", fileext=".txt"))
>
> mu_kalman <- bi_read(bi_kalman_object, "mu")


which gives

> > mean(rdata$x$value[1:201]) - mu_kalman$value[201]
> [1] 2.526777e-05
> > mean(rdata$x$value[1:200]) - mu_kalman$value[200]
> [1] 2.543114e-05
> > mean(rdata$x$value[1:191]) - mu_kalman$value[191]
> [1] 2.655763e-05


So agreement to 10^5 but should I expect better? I guess there is a prior which is causing this discrepancy but how is the prior calculated?

If I just have 11 observations, the discrepancy is worse as one might expect as there is less data to modify the prior.

> > mean(rdata$x$value[1:11]) - mu_kalman$value[11]
> [1] 0.0004699153
> > mean(rdata$x$value[1:10]) - mu_kalman$value[10]
> [1] 0.0005267053
> > mean(rdata$x$value[1:9]) - mu_kalman$value[9]
> [1] 0.0006022458


Dominic Steinitz
dom...@steinitz.org
http://idontgetoutmuch.wordpress.com

[Lawrence Murray]

If you get rid of your transition model and make mu a parameter, not a
state variable, then you should expect that the estimate of the mean
will converge to the posterior mean as the number of samples increases.
Because your prior is uniform, the posterior mean equals the maximum
likelihood estimate (MLE) of the mean in this case (assuming that the
MLE is in the interval [0,1]), and with any other prior you would expect
the posterior mean to converge to the MLE as the number of data points
increases.

Your model here, however, is slightly different. You have a state-space
model that introduces a new mean for each observation, where that new
mean is a slight perturbation of the previous mean. Because that
perturbation is relatively small, the estimate of the mean at the final
time may still approach the MLE of the original model (the same mean
parameter for all observations), but the perturbation will add variance.
The posterior will not contract as quickly.

Lawrence