pmcmc() move outside the support of the likelihood

[erinclancey]

I am using the pmcmc() function in pomp to estimate parameters from an SEIR model. So far it seems to be working great and I can estimate parameters on simulated data. One thing that does happen is that sometimes the algorithm with propose a move to a new parameter value that is outside the support of the likelihood function. For example, it will try to propose a negative Beta and pfilter will return NaN for the value of the log likelihood and the procedure will stop. All I have figured out how to do so far is to dial down the standard deviations in the proposal distribution that controls the random walk so the chain cannot wonder too far. This happens whether I use a prior or not. Do you have any suggestions on how to deal with this issue? Or is this a red flag indicating I am doing something wrong?

[kingaa]

Yes, in MCMC, it's important to propose only parameters for which the likelihood (and the prior probability density) functions are defined as non-negative real numbers. As you've discovered, the simple and general proposal distributions supplied with pomp presuppose that every parameter can be varied across the whole real line (i.e., from -∞ to ∞) without violating some assumption of the model.

In your case (which is very far from uncommon), this is not true. In particular, as you've found, you have a parameter, $\beta$, such that the model only makes sense—and the likelihood is therefore only defined—when $\beta>0$. You haven't told us what proposal distribution you're using in pomp's pmcmc, but you have determined that it is occasionally proposing values of $\beta\le{0}$.

There are at least three solutions to this problem, all of which are classical.

1. The barrier method: redefine your likelihood (and prior if necessary) so that the likelihood is zero (and log likelihood = -∞) whenever parameters are outside of their allowable ranges. In effect, you are expanding the parameter space to the full Euclidean space $\mathbb{R}^n$ by defining the likelihood and prior over this whole range.
2. Modify the proposal: in pomp, you can specify your own proposal distribution. As the pmcmc manual page describes, you can specify the proposal function of your choice, it doesn't have to be one of the prebuilt proposal distributions. It does, however, have to be a symmetric proposal. That is, the probability of going from $\theta_1\to\theta_2$ must be the same as that of going from $\theta_2\to\theta_1$, for all $\theta_1,\theta_2$. Future pomp development may alleviate this requirement (pull requests welcome!!).
3. The transformation method: transform your parameter space so that the parameters to be estimated are allowed to range over the full Euclidean space. If you take this approach, you must also transform your prior distribution appropriately. Note also that the pomp parameter transformation facility won't help you here. That is, pmcmc doesn't make use of parameter transformations, even if you supply them. This is because, if it did, users would routinely assume that their priors would be transformed appropriately, which is not something that can be done for an arbitrary user-supplied prior.

[erinclancey]

Thank you for the response! I am using a mvn.diag.rw() as the proposal. I think a Gaussian proposal is fine for my problem. I would like to use the barrier method. I understand what I need to do, I am just a little confused how to redefine the likelihood function within the pomp syntax. Do I do this within dmeasure(), pfilter(), pmcmc(), or the pomp() object?

[erinclancey]

My guess is I need to modify dmeasure().

[kingaa]

dmeasure reports the likelihood to pfilter and pmcmc. So you can insert something like

if (<bad parameters>) {
  lik = (give_log) ? R_NegInf : 0.0;
} else {
  <likelihood or log-likelihood for good parameters>
}

Perhaps you can share your revised dmeasure with us when you're done?

[erinclancey]

I actually just wrote the following and it worked!

dmeas <- Csnippet("
  if (Beta_p<0 || Beta_w<0 || Beta_m<0 || rho_w<0 || rho_p<0 || k<0) {
    lik=log(0);
  } else {
    lik=dnbinom_mu(reports_p,k,rho_p*H_p,give_log)+dnbinom_mu(reports_w,k,rho_w*H_w,give_log);
  }
")

I don't know how to code in C. I will also try your code. Thank you!

[erinclancey]

This code also works:

dmeas <- Csnippet("
if (Beta_p<0 || Beta_w<0 || Beta_m<0 || rho_w<0 || rho_p<0 || k<0) {
  lik = (give_log) ? R_NegInf : 0.0;
} else {
  lik=dnbinom_mu(reports_p,k,rho_p*H_p,1)+dnbinom_mu(reports_w,k,rho_w*H_w,1);
  lik = (give_log) ? lik : exp(lik);
}
")

[erinclancey]

Everything is working great now! Thank you so much!

[kingaa]

The second one is correct, but the first one is not. Remember that dmeasure has to supply the likelihood when give_log = 0 (i.e., FALSE) and the log likelihood when give_log = 1 (TRUE). As written, your first example will give the wrong result when give_log = 0.