Jacobian adjustment

[Jeremy G Colman]

I am writing a model in which the three parameters of the likelihood function are obtained by a 1-1 continuous deterministic transformation of three independent parameters for which I am able to set specific priors. In Stan I have to include a Jacobian adjustment derived from the determinant of the partial derivatives of the transformation.

I am pretty sure that such an adjustment is needed in my nimble model (because without such an adjustment nimble gives very similar results to the Stan model without the Jacobian adjustment). 

But I don't know how to code the adjustment in nimble. Can anyone help?

[Perry de Valpine]

Hi Jeremy,

Thanks for the question.

I might need some clarification on what you want to achieve.  NIMBLE will use the parameters however you have declared them.  If the three parameters follow a distribution, and then function(s) of the three parameters are used elsewhere in the model, NIMBLE will treat them as following the distribution as declared.  But if what you want is that the distribution should really be defined on the transformed parameters, one way to do that would be to write a user-defined distribution (for the vector of three parameters) and use that in the model.  That user-defined distribution could include a Jacobian adjustment or whatever other math defines the prior that you want.  Writing user-defined distributions is covered in the User Manual.  If I've misunderstood your need here, can you state it more explicitly?

It would also be possible to write a custom sampler that in effect changes the distribution from what is given in the model.  You could copy and modify NIMBLE's samplers to do so, but I think it would be easier to write a custom distribution, so that any sampler will work correctly.

-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 post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/c7bd6048-9883-4488-880b-3db2de18f266%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[David Pleydell]

On Wednesday, 11 July 2018 09:47:37 UTC+2, pdevalpine wrote:

But if what you want is that the distribution should really be defined on the transformed parameters, one way to do that would be to write a user-defined distribution (for the vector of three parameters) and use that in the model.  That user-defined distribution could include a Jacobian adjustment or whatever other math defines the prior that you want.  

It would also be possible to write a custom sampler that in effect changes the distribution from what is given in the model.  You could copy and modify NIMBLE's samplers to do so, but I think it would be easier to write a custom distribution, so that any sampler will work correctly.

There are examples of both these strategies in the code for this paper.  

We write a user defined distribution for log-gamma random variables. These are used to generate dirichlet random variables which are of interest to biologists. We are using these log transformations simply so that default samples can avoid proposing and rejecting negative values. The transform also helped us explore correlations between elements of the dirichlet random variable and other parameters in the model. Here the change of variables correction is in the user defined distribution.

The second strategy (adapting samplers) is adopted for implementing reversible jump MCMC for a finite mixture model with up to four sub-models for the dirichlet rv. Our implementation is adapted from the classic Richardson & Green 1997 RJMCMC paper. Here the Jacobian terms are in the sampler.

Best

David

[Perry de Valpine]

Thanks David!

--
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 post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/fe884769-cc5f-4c0c-b7e1-dba77173f877%40googlegroups.com.

[Jeremy G Colman]

Dear Perry 

Let me be a bit more specific. My priors are defined for three independent parameters q1,q2, q3. (At present these are gamma distributions but that needs to vary in further work). I have written a user-defined distribution for my data given three different parameters mu, sigma and xi. The transformation, f, from the q1,q2,q3 to the mu, sigma , xi is 1-1 deterministic and continuous. I can derive analytically what the joint prior density of the mu, sigma, xi is given my priors for the q1, q2, q3. I attach a pdf showing that derivation and that shows that a Jacobian factor has to be applied.  

What I need is a way of coding the joint prior density of the mu, sigma, xi. I f I just write q1 ~ gamma() etc and then mu = f(q1,q2, q3) etc no Jacobian adjustment is applied and I get the wrong answer. In Stan I just add the log of the Jacobian adjustment to the log likelihood at each iteration. But I don't know how to achieve the equivalent in NIMBLE. Doubtless it will be obvious when I see it!

Very many thanks for your help.

Jeremy

[Perry de Valpine]

Thanks Jeremy.  

I interpret what you are saying as wanting a custom distribution for the triplet (mu, sigma, xi) that you have derived from independent distributions for (q1, q2, q3) along with the function relating the two triplets.

Here are two ways you could implement that in NIMBLE, both using user-defined distributions.

Method 1:  Define (mu, sigma, xi) as a vector of three parameters, say theta[1:3].  Use

theta[1:3] ~ dMyCustomDistribution(parameters of the distributions of q1, q2, q3)

Then in dMyCustomDistribution, you can write the full pdf of theta[1:3], which will include pieces from the distributions of (q1, q2, q3) as well as the Jacobian, as you have derived.  In this method, the distributions for (q1, q2, q3) would not appear in the model code.  They would simply appear in dMyCustomDistribution.

Method 2: You could define a dummy data node, whose value need not even matter, in order to build an additional log density term into the graph.  That term would be the Jacobian need you.  This would look like the following

q1 ~ dgamma( ... )

q2 ~ dgamma( ...)

q3 ~ dgamma( ...)

dummy ~ dMyJacobian(q1, q2, q3, parameters of q1, q2, q3)

Now there will be a node "dummy" in the graph, with parent nodes q1, q2, q3 as well as their parameters.  If you write dMyJacobian (again, following our User Manual on how to write your own distribution) to return the log Jacobian term you need, it will be added into the likelihood.  For example, let's say you have a "RW" (adaptive random-walk Metropolis-Hastings) sampler on a parameter "alpha" of q1.  The stochastic dependencies of alpha will include q1 and dummy.  That means that the sum of log probabilities used in the Metropolis-Hastings acceptance ratio will include a term from dgamma for q1 and a term for dummy that provides the Jacobian.  When you write dMyJacobian, its first argument will be a deviate "x", which will be passed as "dummy", but you need not use it in the calculation.  Other samplers or algorithms should work too, because the graph defined by these declarations include the pieces you need.

Let me know if that gives you enough to move forward or if that was too rough an explanation.

-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 post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0c2a6a91-a191-4649-a707-ab6c194d92a9%40googlegroups.com.