How to get the particles?

[Dominic Steinitz]

I am running almost the simplest model possible: estimating the mean of sample drawn from a normal. I treat the mean as the state with the standard deviation as a parameter.

> 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);
> }
> }

I’d like to draw the particle paths for each Markov chain step but if I do e.g.

> bi_object$run(add_options = list(
> "end-time" = T,
> noutputs = T,
> nsamples = 8,
> nparticles = 32,
> seed=42,
> nthreads = 1),
> verbose = TRUE,
> sampler="smc2",
> stdoutput_file_name = tempfile(pattern="pmmhoutput", fileext=".txt"))
>
> bi_file_summary(bi_object$result$output_file_name)


I get

> Summary of file /private/var/folders/3p/m593dprn5snbjz6sc45c77f80000gn/T/RtmpqWLBkv/MV_outputeaff7423f4da.nc
> [1] "2016-12-20 09:46:38 GMT"
> File /private/var/folders/3p/m593dprn5snbjz6sc45c77f80000gn/T/RtmpqWLBkv/MV_outputeaff7423f4da.nc (NC_FORMAT_NETCDF4):
>
> 8 variables (excluding dimension variables):
> double time[nr] (Contiguous storage)
> double eta[np,nr] (Contiguous storage)
> double mu[np,nr] (Contiguous storage)
> double sigma[np] (Contiguous storage)
> double loglikelihood[np] (Contiguous storage)
> double logprior[np] (Contiguous storage)
> double logweight[np] (Contiguous storage)
> double logevidence[nr] (Contiguous storage)
>
> 2 dimensions:
> nr Size:51
> np Size:8
>
> 3 global attributes:
> libbi_schema: SMC
> libbi_schema_version: 1
> libbi_version: 1.2.0


The LibBI manual tantalisingly talks about `np` being the number of particles when a particle filter is chosen. I assume SMC (`libbi_schema: SMC`) is somehow not a particle filter? I have also tried `pmmh` and get `libbi_schema: MCMC` again with the dimensions of the time steps and the chains.

Anyhow the TL;DR is: how do I get hold of the particle paths? In case it isn’t obvious I am using the `R` `rbi` package.

Many thanks,

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

[dom...@steinitz.org]

Path storage is apparently implemented in LibBI:

• P. E. Jacob, L. M. Murray, and S. Rubenthaler. Path storage in the particle filter, 2013. Statistics and Computing, to appear. [doi] [arXiv]

I just can’t see how to get hold them

[Sebastian Funk]

Hi Dominic,

The particle filter yields one trajectory per sample: If you do

> trajectories <- bi_read(bi_object, "mu")

You get the 8 trajectories for mu (along the 'np' dimension).

Seb

[dom...@steinitz.org]

Dumb question - thanks for your indulgence

On 20 Dec 2016, at 10:14, Sebastian Funk <sebasti...@lshtm.ac.uk> wrote:

trajectories <- bi_read(bi_object, "mu")

[dom...@steinitz.org]

Ah no - I still am missing something. I should have 32 particles so 32 paths per Markov chain step? I think what this gives me is the mean path per Markov chain step (I think median would be better but I am guessing it is the mean).

[Sebastian Funk]

See section 3.1.2 in Murray(2013):
https://arxiv.org/abs/1306.3277

The particle filter returns a single sample from the posterior distribution p(X|\theta, y) via ancestral tracing. I don't think there is a way to retrieve the individual particles/weights at each data point form a libbi run.

dom...@steinitz.org wrote:
> Ah no - I still am missing something. I should have 32 particles so 32 paths per Markov chain step? I think what this gives me is the mean path per Markov chain step (I think median would be better but I am guessing it is the mean).
>
>> On 20 Dec 2016, at 10:16, dom...@steinitz.org wrote:
>>
>> Dumb question - thanks for your indulgence
>>

>>> On 20 Dec 2016, at 10:14, Sebastian Funk <sebasti...@lshtm.ac.uk <mailto:sebasti...@lshtm.ac.uk>> wrote:
>>>
>>>> trajectories <- bi_read(bi_object, "mu")
>>
>> Dominic Steinitz

>> dom...@steinitz.org <mailto:dom...@steinitz.org>
>> http://idontgetoutmuch.wordpress.com <http://idontgetoutmuch.wordpress.com/>

[Lawrence Murray]

Indeed, there is no way to retrieve all of the particles used for each
sample, only the single path that is traced out. The storage required to
output everything could be huge for some models.

What you can do is run the "filter" command for one set of parameters,
and the output file for that will contain all of the particles and
weights. Given this is a simple model, perhaps you could even implement
a simple PMMH in R that calls "libbi filter" on each iteration with the
proposed set of parameters?

[Dominic Steinitz]

I had just reached that conclusion and was trying to find out what I had to say in `rbi` to do this *but* I will get different paths unless I can find some way of saving the RNG state as the random numbers used will be different (not too much of a problem for my toy example).

> --
> You received this message because you are subscribed to the Google Groups "LibBi Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to libbi-users...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Dominic Steinitz]

This works. All I need to do now is figure out how to set sigma to the values I got in the PMMH run I did previously.

> bi_filter_object <- libbi(client="filter", model=MV, obs = synthetic_dataset)
>
> bi_filter_object$run(add_options = list(

> "end-time" = T,
> noutputs = T,

> nparticles = 32,
> seed=42,
> nthreads = 1),
> verbose = TRUE,

> stdoutput_file_name = tempfile(pattern="pmmhoutput", fileext=".txt"))
>

> bi_file_summary(bi_filter_object$result$output_file_name)
> sigma_filter <- bi_read(bi_filter_object, "sigma”)

[Sebastian Funk]

If you have a variable `bi_object` used to run libbi,
> samples <- 100
> bi_object <- libbi(client="sample", target=posterior, nsamples=samples, {...})
you can feed this directly into the new run:

> bi_filter_object$run(init=bi_object, {...})
will use the samples of the PMMH run.

bi_filter_object$run(init=bi_object, init_np=sample-1, {...})
will use the last sample of the PMMH run.

Or if you want to fix sigma to some other value, you could use
bi_filter_object$run(init=list(sigma=0.5), {...})

The same would work for `input` instead of `init` in all these cases.