An Initialisation Puzzle

[Dominic Steinitz]

I have the world’s simplest LibBI script:

> 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 I am generating synthetic data thus:

> T <- 10
>
> MV <- bi_model("MeanAndVar.bi")
> synthetic_dataset <- bi_generate_dataset(end_time=T,
> model=MV,
> seed=44,
> init=list(mu=0.5, sigma=0.1))
>
> rdata <- bi_read(synthetic_dataset)

Which gives

> > rdata
> $eta
> time value
> 1 0 0.000000e+00
> 2 1 -6.565509e-07
> 3 2 -3.938302e-07
> 4 3 -4.371722e-07
> 5 4 1.363798e-06
> 6 5 -7.483461e-07
> 7 6 -4.805008e-07
> 8 7 1.106146e-06
> 9 8 -1.670766e-06
> 10 9 6.707771e-07
> 11 10 5.461491e-07
>
> $sigma
> [1] 0.1
>
> $mu
> time value
> 1 0 0.5000000
> 2 1 0.4999993
> 3 2 0.4999989
> 4 3 0.4999985
> 5 4 0.4999999
> 6 5 0.4999991
> 7 6 0.4999986
> 8 7 0.4999998
> 9 8 0.4999981
> 10 9 0.4999988
> 11 10 0.4999993
>
> $x
> time value
> 1 0 0.9231253
> 2 1 0.5377911
> 3 2 0.4956203
> 4 3 0.4757358
> 5 4 0.5661350
> 6 5 0.5317747
> 7 6 0.5319090
> 8 7 0.4755841
> 9 8 0.3562161
> 10 9 0.3595353
> 11 10 0.5349840

Notice the first value is 0.9231253 which is over 4 standard deviations away from the mean. I’ve tried a few different seeds and this phenomenon (of the first value being an improbably large way from the mean) always happens. In fact with a seed of 43 I get -1.3211570 which 8 standard deviations away!

Does anyone have any idea why this happens?

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

[Sebastian Funk]

I think this might be due to the order in which model blocks and initialisation files are sampled. In your case the 'initial' and 'observation' blocks get processed before mu is updated from the init file (and written to the output file as the value of mu at time 0). It looks like your x at time 0 is distributed according to the distributions given in the model file, whereas at later times it's distributed with a mean around 0.5 and standard deviation of 0.1 as expected.

[Lawrence Murray]

Indeed, Seb is correct. If you run the following:

libbi rewrite --model-file MV.bi --with-transform-obs-to-state

you will see the internal view of your model after transformation (when
you use --target joint, as when generating data, the
--with-transform-obs-to-state flag is automatically applied). Note how
the body of the observation block is shifted into your initial and
transition blocks, and your init file will be read *after* the initial
block is executed. So the first observation is using the value of mu
simulated from the prior, not that from your init file.

This behaviour would be tricky to change. A workaround would be to make
the initial value of mu a parameter itself. Call this mu0, say. Then in
your parameter block:

mu0 ~ normal(0.0, 1.0)

and in your initial block

mu <- mu0

Then put mu0 in your init instead of mu.

Cheers,

Lawrence