How can I prevent parameters getting negative in pmcmc?

[tillweigt]

I'm fitting a standard sv model.

I get negative parameters for the sd of the transition equation.

I use uniform priors with positive support for this parameter.

How can I prevent the parameter from getting negtiave?

[kingaa]

I suggest either writing a proposal that eschews negative s.d.s or reparameterizing (but you'll have to adjust the prior to match).

[tillweigt]

Thank you for your fast answer.

I thought that dprior would do that as a negative proposal would lead to negative infinity.

But it does not work that way.

Proposals in pmcmc are always normal. Thus, I do not see how to do it this way.

Reparameterizing would be possible, but a bit tedious.

Ok, thank you. I thought there would be an easy setting with dprior or something like that.

[kingaa]

Perhaps I responded too fast! I'm afraid I'm not familiar with "sv" models. Perhaps that's just jargon I don't recognize.... Can you say a bit more about the model for the uninitiated? In particular, what is the prior distribution?

If your dprior gives zero probability to negative values, then those proposals should be always rejected. Are you saying that you are observing them being accepted?

Yes, the proposal distributions that are provided in pomp are all normal, but you have the option of providing a proposal of your own, as an R function. See the documentation for pmcmc. In your case, however, this may not be enough, since the underlying Metropolis Hastings algorithm assumes that the proposal is symmetric. Support for non-symmetric proposals is something on the TODO list that I haven't got round to adding.

[tillweigt]

Hi again,

thank you for your fast reply.

Right now it works. Sorry. I really did not write you directly.

Here is the minimal example.

By the way, I do have another question.

Is it possible to start with a vector or a matrix of state particles. I only foun the possibility to define a distribution.

Now, I wrote my own rinit function with metaprogramming. This works. I just wondered, if there is a cleaner approach.

Kind regards, Till

library(pomp)

N = 100 Np = 50 Nmcmc = 1000

modelDefinition = pomp( data = data.frame( y = rep(NA, N), time = 1:N ), times = "time", t0 = 0, rmeasure = Csnippet("y = rnorm(0, sqrt(exp(x1)));"), dmeasure = Csnippet("lik = dnorm(y, 0, sqrt(exp(x1)), give_log);"), rprocess = discrete_time( Csnippet("x1 = rnorm(x1, p1); "), delta.t = 1 ), obsnames = "y", statenames = c("x1"), params = c(p1 = 0.1, x1.0 = 1.0), paramnames = c("p1"), dprior = Csnippet("lik = dunif(p1, 0, 10, give_log);") )

modelDefinition@data[1, ] = simulate(modelDefinition)@data[1, ]

modelDefinition@params["p1"] = 1.0

estimate = pmcmc( data = modelDefinition, Nmcmc = Nmcmc, Np = Np,

proposal = mvn.diag.rw(c(p1 = 0.5)),

proposal = mvn.rw.adaptive( rw.sd = c(p1 = 0.5), scale.start = 1, shape.start = 100 ), params = c(p1 = 1.0, x1.0 = 1.0) )

[kingaa]

Thanks for the MWE: it makes everything much clearer.

You ask about initializing with a given matrix of state particles. May I ask out of curiosity, why do you wish to do this?

In effect, you are making the initial values of the state particles part of your prior, so perhaps a more principled approach would be to conceptualize the prior as consisting of a mixture of Dirac delta functions at the points in your matrix, then write an rprior that samples from them. One could also use rinit for this purpose: the rinit function would sample from the columns of the matrix. If you wanted to dispense with the random element for some reason, it could be done, but I would be concerned about the validity of the resulting inference: I would check to make sure that the resulting model made sense as a model.

P.S. If you use the Issues web interface (as opposed to the email reply), our correspondence can be formatted with markdown, which, for example, will allow code to be formatted in a more readable way.

[tillweigt]

We are working in a live environment.

The workflow is pmcmc, then pfilter for some datapoints, then again pmcmc, then pfilter for some datapoints, and so on.

Every time a new datapoint comes in, we apply one step of pfilter using the last state particles. That's why we do not need something stochastic as the old particles were just resampled.

I do not understand your matrix solution. The problem we had with pomp is that params (in the pomp object) has to be a named vector. So I can not pass a matrix to params. Otherwise it would be easy.

Normally, we should use SMC^2 or something like that due to the live environment, but there is no code in R.

If you have time for another question. Is this the right way to define the dprior for, say, two parameters named p1 and p2:

Csnippet("lik = dunif(p1, 0, 1) * dunif(p2, 0, 1) ;")

I really appreciated your help and fast answers. Thank you very much.

[kingaa]

I do not understand your matrix solution. The problem we had with pomp is that params (in the pomp object) has to be a named vector. So I can not pass a matrix to params. Otherwise it would be easy.

No, but R functions can access objects in the environment in which they were created. So you could create a matrix to hold your states, then access this within your rinit function (which would not then depend on parameters, I suppose).

Is this the right way to define the dprior for, say, two parameters named p1 and p2:

Csnippet("lik = dunif(p1, 0, 1) * dunif(p2, 0, 1) ;")

Recall from the documentation (?prior_spec) that the dprior has to return the value of the density, or its log, according to whether give_log is set to 0 or 1, respectively. The code above handles only the first case. You will have to handle the other case, e.g.,

Csnippet("lik = (give_log) ? dunif(p1,0,1,1) + dunif(p2,0,1,1) : dunif(p1, 0, 1, 0) * dunif(p2, 0, 1, 0) ;")

[tillweigt]

Okay, that's possible. We solved it now with some metaprogramming without going to the parent environment.

Thank you for the dprior definition.

[kingaa]

I am closing this now. Feel free to re-open if more discussion is warranted.

[tillweigt]

Thank you very much for your help

Am Fr., 11. Sept. 2020 um 17:54 Uhr schrieb Aaron A King < notifications@github.com>:

I am closing this now. Feel free to re-open if more discussion is warranted.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/kingaa/pomp/issues/124#issuecomment-691176141, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGDFEPIMYG3FKFI5BXEFUQLSFJBZ3ANCNFSM4RHKNSQQ .