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