Help debugging convergence problems with MCEM

[Russell Almond]

I've set up a model to use with MCEM.  I can get the bootstrap filter working with this model and it compiles, however, when running the EM algorithm it diverges:

```

## Compilation
brownianModel <-
  nimbleModel(brownianModelCode1,
              constants=list(N=N,J=J,NT=NT,
                             deltaT=deltaT,
                             tasks=tasks),
              dimensions=list(p=c(N,NT,3L)),
              inits=bminits,check=FALSE,
              data=list(y=y),
              buildDerivs = TRUE)

cBrownianModel <- compileNimble(brownianModel,showCompilerOutput=FALSE)
brownianMCEM <- buildMCEM(model=brownianModel,latentNodes="theta")
cBrownianMCEM <- compileNimble(brownianMCEM,project=brownianModel)

## Run MCEM (this step fails)
cBrownianMCEM$findMLE(initM=1000)
```

This last gives an error:

```

## Error in cBrownianMCEM$findMLE(initM = 1000) :

##  Non-finite log-likelihood occurred; the MCEM optimization cannot continue.  
## Please check the state of the compiled model (accessible as 'name_of_model$CobjectInterface') 
## and determine which parameter values are causing the invalid log-likelihood by calling 
## 'calculate' with subsets of the model parameters 
## (e.g., 'name_of_model$CobjectInterface$calculate('y[3]')' 
## to see if node 'y[3]' is the cause of the problem). Note that if your model is maximizing 
## over parameters whose bounds are not constant (i.e., depend on other parameters), 
## this is one possible cause of such problems; in that case you might try running the MCEM 
## without bounds, by setting 'forceNoConstraints = TRUE'.

```

This error message doesn't reflect the current version of nimble, as `cBrownianMCEM` does not have a
`$CobjectInterface` field, although `brownianMCEM$.CobjectInterface` is an object. But this
doesn't have a `$calculate()` method.

The `brownianModel` and `cBrownianModel` objects do have a `calculate()` method, but I'm
not sure that their parameter values are referencing.

Q1. How do I access the values at the point in which MCEM blew up to diagnose the problem?

Q2. Is there a way to turn on tracing for the MCEM algorithm?

There seems to be an undocumented `optimControl` parameter which can be set in `control`,
but tracing at that level would trace within the M-step and I want to trace the EM algorithm.
(As that error message clearly needs to be fixed, I posted Issue #1535 on github. I've got the complete
model and some test data there if you are curious).

Thanks,
--Russell

-- 
Russell Almond
https://ralmond.net/

[Perry de Valpine]

Hi Russell,

Thanks for the question.

Some clarifications: buildMCEM will by default build an MCMC with default samplers for the MCMC steps of the MCEM iterations. You can provide your own MCMC configuration function if you want to change that. But it would not by default build a particle filter, so I wasn't clear what you were referring to (or simply if you were referring to an entirely separate case but using the same model).

The <name_of_model> in the error message would be referring to the model object, i.e. brownianModel, not the algorithm objects BrownianMCEM (uncompiled) or cBrownianMCEM (compiled). After calling compileNimble, then brownianModel$CobjectInterface should exist. But that will simply be cBrownianModel, which in your workflow you already have your hands on, so you do not need to use brownianModel$CobjectInterface to find it. Since the compiled algorithm will use the compiled model, the suggestion is that you can inspect the compiled model (cBrownianModel) to look for where there is a non-finite log probability (called log-likelihood in the message). E.g. you can look at node values (cBrownianModel$x or whatever) or calculations (cBrownianModel$calculate("x") ).

I don't see a built-in way to trace the MCEM iteration by iteration, but here are a few ideas:

(i) Make a small version of the problem and run MCEM uncompiled. Then you can set debug(brownianMCEM$findMLE) and step through it.

(ii) Make a copy of the source code and modify it. Copy buildMCEM from the source code. Rename it in your own local code file (also change the "name=" argument). When you run it the setup code may break because it may use non-exported functions. Insert nimble::: in front of those to access them from the nimble namespace. You might also need to copy, rename, and insert nimble::: as needed in other nimbleFunctions used by buildMCEM. Define a nimbleRcall and put browser() inside the called R function (or debug() it). Insert calls to that nimbleRcall inside the findMLE method's source code. Then you will be thrown from the compiled algorithm into a browser() where you've inserted those calls, and you can inspect the compiledModel to see what is going on with it. You could also insert print statements etc.

(iii) Look up the documentation in the User Manual or help(nimOptim) about providing your own method to optim. (In help(nimOptim), this is described as "experimental" but should work fine). You can provide an R function that simply makes your own call to a regular optimizer (e.g. optim) but also lets you generate output, insert browser() statements, or do other debugging).

HTH

Perry

 

--
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 visit https://groups.google.com/d/msgid/nimble-users/5a47513f-c23b-4ba1-b0a3-68b7c00bf815%40gmail.com.

[Russell Almond]

Thanks for your reply, this is helpful.

I think I was a bit confused as the MCEM algorithm is listed under the SMC chapter in the user manual. 

Looks like what I wanted to do was not coded anywhere.  As my model is pretty much a straight HMM, I should be able to run one of the SMC filters to get an approximate expected distribution for my latent variable, and then use that to build expected values of the sufficient statistics, and the use optim to do M-steps to get updated parameters for my emission models and growth models. 

I think I'll probably start with code for my example rather than try to work through all of the various cases that MCEM supports.

Do check the error message that MCEM is giving when then model diverges as it's advise does not correspond to the current implementation.

  --Russell

[Perry de Valpine]

Thanks Russell.

We'll take a look at how User Manual Chapter 8 is organized. I think the intention is that it covers "other algorithms" and the title names five algorithms without trying to imply one encompasses the others, but I get that it was not clear to you. That's helpful.

In case it's helpful, there are marginal HMM distributions in the nimbleEcology package.

We'll also take a look at clarifying which object is being referred to in the error message you saw.

Good luck with your project.

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/38c5e2b2-7e0b-49bc-a797-b864f3169485n%40googlegroups.com.