Initializing each EnKF ensemble member

[Quinn Thomas]

Hi all,

I am exploring the use of Nimble in iterative ecological forecasting.  The EnKF is particularly exciting in the nimbleSMC package. In the examples provided, the initial conditions for the latent states are defined as distributions within the nimble code and the number of ensemble members is set at the run step (Cmy_enKF$run(m = 1000)).  Is it possible to pass the initial value for each ensemble member (i.e., the last step of a previous run of the EnKF is the initial conditions for a new run of the EnKF).  This would allow an iterative approach to running the EnKF.

Here is the code I am running with the initial conditions defined within the nimble code

library("nimbleSMC")

## For illustration only. 

exampleCode <- nimbleCode({

  x0 ~ dnorm(0, var = 1)

  x[1] ~ dnorm(a * x0, var = 1) 

  y[1] ~ dnorm(x[1], var = .5) 

  for(t in 2:10){

    x[t] ~ dnorm(a * x[t-1], var = 1)

    y[t] ~ dnorm(x[t], var = .01)

    }

  })

y <- rnorm(10)

model <- nimbleModel(code = exampleCode, data = list(y = y) , inits = list(x0 = 0, x = rnorm(10)))

model$a <- 0.8

my_enKF <- buildEnsembleKF(model, 'x',

                           control = list(saveAll = TRUE, thresh = 1))

## Now compile and run, e.g.,

Cmodel <- compileNimble(model)

Cmy_enKF <- compileNimble(my_enKF, project = model)

Cmy_enKF$run(m = 1000)

enKF_X <- as.matrix(Cmy_enKF$mvSamples, 'x')

plot(enKF_X[1, ], type = "l")

for(i in 2:100){

  points(enKF_X[i, ], type = "l")

}

points(y, col = "red")

Thanks for any insight from the group!

[Perry de Valpine]

Hi Quinn,

Thanks for making this request.  It's great to see interest in the sequential Monte Carlo methods.

I made a rough draft for how to do this on a branch on nimble-dev/nimbleSMC on GitHub.  It took only a few small tweaks.  The methods are all written as nimbleFunctions, so you can inspect the changes without too much trouble in the source code.  Below is code to install from the branch and an expansion of the example to see (I think) that it works.  Can you give it a try and see if it does what you were looking for?

Are there any easy suggestions on how people will want to use a feature like this beyond what I drafted (one can set the first samples and then tell the EnKF to use them)?

I have a question for anyone on this list about typical use cases for Ensemble KF.   Of all the SMC methods we implemented (which are described here: https://arxiv.org/abs/1703.06206), among the developers involved we have least experience with EnKF.  Other SMC methods return both an approximate log likelihood and the filtered state samples.  The papers we worked from on EnKF emphasized state estimation, and that's what we implemented.  In practice do people make use of an approximate log likelihood from EnKF?

-Perry

library(devtools)
install_github(repo = "nimble-dev/nimbleSMC", ref = "support-recursive-use", subdir = "packages/nimbleSMC")

library("nimbleSMC")
## For illustration only.
exampleCode <- nimbleCode({
  x0 ~ dnorm(0, var = 1)
  x[1] ~ dnorm(a * x0, var = 1)

  y[1] ~ dnorm(x[1], var = .2)
  for(t in 2:10){
    x[t] ~ dnorm(a * x[t-1], sd = .2). # I changed these sds to create a more obvious visual check below
    y[t] ~ dnorm(x[t], sd = .2)
    }
  })

set.seed(3)

y <- rnorm(10)

model <- nimbleModel(code = exampleCode, data = list(y = y) , inits = list(x0 = 0, x = rnorm(10)))

model$a <- 0.8
my_enKF <- buildEnsembleKF(model, 'x',
                           control = list(saveAll = TRUE, thresh = 1))
## Now compile and run, e.g.,
Cmodel <- compileNimble(model)
Cmy_enKF <- compileNimble(my_enKF, project = model)
Cmy_enKF$run(m = 1000)
enKF_X <- as.matrix(Cmy_enKF$mvSamples, 'x')

plot(enKF_X[1, ], type = "l", ylim = c(-2, 2))

for(i in 2:100){
  points(enKF_X[i, ], type = "l")
}
points(y, col = "red")

lastSamples <- enKF_X[, 'x[10]', drop=FALSE] # pull out last samples as a matrix
Cmy_enKF$setFirstSample(lastSamples)        # set them as the first samples
Cmy_enKF$run(m = 1000, skipFirst = 1)          # run again, skipping t = 1
enKF_X2 <- as.matrix(Cmy_enKF$mvSamples, 'x')
for(i in 1:100){
  points(enKF_X2[i, ], type = "l", col = 'cyan')
}

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/2cdf31c4-8c7c-40fe-b225-1b8e90445111n%40googlegroups.com.

[Quinn Thomas]

Hi Perry 

This is great and my initial test worked fine.  In response to your question, the vanilla EnKF assumes equal weight of all the ensemble members so an approximate log likelihood does not apply.  

One question about your update:  what happens if the ensemble size of Cmy_enKF$setFirstSample(lastSamples) differs from m in Cmy_enKF$run(m = 1000, skipFirst = 1)?  It seems like it forces m to equal  Cmy_enKF$setFirstSample(lastSamples), which seems to me to be the correct behavior.   In the example, I was confused a bit because m was included in the argument but doesn't seem necessary (it did break the run if excluded).

Thanks for the quick development!

Quinn  

[Perry de Valpine]

Hi Quinn,

Great that it works for you.  Aha, thanks for explaining that about EnKF.

You are right about the values of m.  It would not make sense to use a different value in run than the size that was used in setFirstSample.  Both methods resize the modelValues object internally, so I think the practical effect will be that if m in run is larger than in setFirstSample, there will be some nonsense values used, and if it is smaller, simply not all of the "lastSamples" will be used.  But you said it broke the run, so I'll check on that later.

Unfortunately we have very limited default-value behavior supported for compiled nimbleFunctions.  We could have a default for m that is -1, which if provided means that an internal default should be used, and that internal default could be updated by setFirstSample.  The default internal default could be 100, the current default for m.  Does that makes sense?

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/6bf6fe4e-5e27-4cc5-a082-85546222d8a9n%40googlegroups.com.

[Quinn Thomas]

Hi Perry

My response was incorrect.  The code did not break when removing the m from the function call or changing m do a different value than in the setFirstSample.  It appeared to maintain the same number of samples.

Quinn

[Quinn Thomas]

Hi Perry

A follow-up question will these updates improve the other sequential methods by allowing them to restart with the ensemble states and log likelihoods? 

Quinn

[Perry de Valpine]

Yes I was thinking about that too.  I think it would apply to bootstrap particle filter, auxiliary particle filter, and Liu-West particle filter (which is not very practical anyway).  I don't think it's relevant for iterated filtering or particle MCMC.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/39fe54ac-7e0c-4ef3-b6f2-50447f67fcb9n%40googlegroups.com.

[Quinn Thomas]

Another thing to consider that would be very powerful for the ecological forecasting applications.  In your current implementation, the state variable x is initialized at x[1], the "forecast step" is applied to transition to x[2], and then the analysis step is applied the "forecasted" x[2].  So the first step is a "forecast step".  This assumes that either x[1] is initialized following an analysis step at the end of the previous run (or a forecast step if no data was available at the last step).  In an ecological forecasting application there may not have been data at the last step of the previous run (because it is a forecast of the future) but data is now available at reinitialization (because it is reinitialized one actual data later).  This means that one would want to initialize and then perform an analysis using data at t =1.  Is this feasible to implement?  (in case you are wonder who might use these features - I am leading an NSF-funded Research Coordination Network focused on ecological forecasting using NEON data.  We are actively looking for (and developing if necessary) tools for near-term iterative forecasting, especially Bayesian tools (see ecoforecast.org)

Thanks!

Quinn

 

[Perry de Valpine]

Hi Quinn,

Thanks for thinking more about this.  Let me say in more detail what I believe happens now and see if that helps sketch out improvements.

With skipFirst = FALSE (i.e. the first run), the state variable at x[1] is initialized from its prior (if you want to look at it that way) as if that is its forecast.  If there is data that depends on x[1], the filter calculations are done to create the sample of x[1] given data.  I think this is what you are calling the "analysis step" and what you were saying does not happen for x[1], by which I think you are talking about the drafted skipFirst = TRUE case (see below) but wasn't totally sure.  (Note that there can be arbitrary data that depend on each x[t], and they need not be the same in structure or distribution for each t).  If one does not want filtering at the first step or the last step as a matter of model structure, one can simply have no data nodes that depend on the state variable at one or both of those times.

Here is how I implemented your suggestion for recursive operation.   I assumed the final sample from the first EnKF run gets filtered by any observations.  That is then set as the first sample for the second EnKF run.  Separately, the skipFirst argument for the second run allows both sampling from the prior and filtering on the data to be skipped for x[1].  Hence, my assumption was that filtering happens at the end of one run and is not desired at the beginning of the next run.

I think one could achieve what you're suggesting now either by careful choice of model or by adding some switches to the algorithm.  On the model side, if you simply have a final time step with no data, you'll end up with unfiltered projections at the last time step of states.  At the first time step, one trick would be to have state "dynamics" from x[1] to x[2] that are essentially deterministic, like x[2] ~ dnorm(x[1], sd = 1e-6).  Then one could use the system I set up yesterday and make the observations at x[2] be the "new data" that really correspond to the final (unfiltered) time-step from the first EnKF run.  The alternative of adding some switches could be to provide control over whether filtering is done at the first and/or last time-steps, so that even if there are relevant data nodes, the filtering can be skipped.  This may have advantages of greater clarity (arguably) and greater flexibility to make the most of a single model and algorithm once they have been compiled (because it could be used in multiple ways rather than having the algorithm intention baked into the model structure).  At least, I think those ideas would work, but I haven't tried them!

Let me know if that all makes sense and what you think would be useful.  We could also take this off-list to go granular and then come back and share what we come up with.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5062e081-c9d0-4ba6-b8ee-bdfccb680006n%40googlegroups.com.