How to do a direct increment of the log posterior in greta

[jae]

Hi there,

I am trying to implement an operation identified previously here in the context of a “bym2” model:

I need help with some spatio-temporal analysis

Hi dr2p, I see several options: Spatio-temporal models are (naturally) supported in greta via the extension package greta.gp. This extension uses Gaussian Processes to estimate covariance structure across space / time / whatever. This would however lead to different model than the one you pasted. You could check this forum thread for some infos on how to get started. You could check this Stan BYM model, that I guess should be possible to translate into greta. The only tricky line there wo…

Specifically, I am trying to adjust/account for a “summary statistic” which in STAN is permitted via a direct manipulation of the joint log density (“target”):

 target += -0.5 * dot_self(phi[node1] - phi[node2])

I understand this is not a supported operation in greta. From looking at similar questions, it seems I would need to create a new distribution as indicated by:

Usage of variable and opt

Yup, as you say, greta has a slightly higher-level interface than stan, so the user API doesn’t let users increment the joint log density (target, in stan) of the model directly. That’s a deliberate design decision because a) it’s easy to mess up the whole log-Jacobian adjustments thing, and b) I want to motivate people to develop extensions packages for greta, so that new functionality like that can be written once, well, and distributed to other users. The downside of this approach is that it…

and possibly use “greta.funprior” as identified here:

Prior Information On Quantities of Interest AND On Parameters

Right, this is actually something I’ve been thinking about recently, interesting that the post just appeared! In the first version of the post, Gelman casually mentioned needing to account for the log Jacobian adjustment too. The post has since been updated to make clear that in the linear case that doesn’t happen to be necessary, but in general it is. The need (in general) to do a log Jacobian adjustment is the reason that assigning a distribution on the result of an operation is specifically …

Sadly, this is beyond my poor abilities. Would any of you be aware of a better solution or, provide more guidance on how I might try to get this done?

[jae]

Thinking it about it some more, would using the zeros trick from Bugs be a potential solution? Something like:

   zeros <- as_data( rep(0L, N) )   
   nlk <- -log( 0.5 * ( phi[node1] - phi[node2] )^2 )
   distribution( zeros ) <- poisson( nlk )