dsum JAGS

[Jose Jimenez]

Hi. Thank you very much for your NIMBLE software. It works fast and smoothly.
I haven't found JAGS "dsum" distribution. What is the equivalence of in NIMBLE?

Thank you,
Jose

[Chris Paciorek]

Hi Jose,

You can probably use our new (as of version 0.4) dconstraint functionality to mimic what dsum does. However, note (as also noted in the JAGS manual) that introducing an equality constraint can cause problems in MCMC sampling as a proposal would need to obey the constraints exactly.  So while you can probably introduce the constraint, I suspect you will not be able to fit the model.  You may be able to reparameterize your model to automatically obey the constraint and reduce the number of parameters by one.

--
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/3e0fd0f8-1c3f-4922-ab06-9604bd701ce8%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jose Jimenez]

Thank you Chris for your fast answer.

Regards,

Jose

[nimble.stats]

Hi Jose,

I think it would also be possible to set up a sampler that would respect the constraint.  Adding custom samplers is covered in section 7.8.1 of the User Manual (version 0.4), but we could provide a more specific example if that would help.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/9a5ab9dc-23c1-4506-a1b7-0e298c848a1c%40googlegroups.com.

[Jose Jimenez]

Hi Perry. That would be great.

Thank you in advance,
Jose

[nimble.stats]

Here is an example.  I want to echo what Chris said: I’m not sure this is the best path since representing constraints in different ways can lead to different results.  But for illustration here is an example:

## Example of adding a sampler that obeys a linear constraint

library(nimble)

## Make a model with two nodes, a and b.

## Suppose we want to sample with constraint a + b = c, where c will be data

## and we'll take care that we provide initial values that follow the constraint

sillyModel <- nimbleModel( nimbleCode({a ~ dnorm(0, 1); b ~ dnorm(0,1); c ~ dnorm(0, sd = 1000)}), inits = list(a = 0.3, b = 0.4), data = list(c = 0.7) )

## We don't really need a distribution for c.  It is just written as a way to get c as a node into the model.

## Maybe we should make a better way to do that.

## To do constrained sampling, we'll work temporarily in the space

## w1 = (a + b) ## i.e. w1 = c

## w2 = (a - b)

## Then we'll sample w2 as a random walk Metropolis-Hastings and convert back to a and b

RW_constrained <- nimbleFunction(

    contains = sampler_BASE,

    ## assume target is a vector of two nodes and control$constraint gives the constraint node

    setup = function(model, mvSaved, target, control) {

        scale = control$scale

        constraintNode <- control$constraintNode

        calcNodes <- model$getDependencies(target)

    },

    run = function() {

        target_vals <- values(model, target)

        current_w2 <- target_vals[2] - target_vals[1] ## could use linear algebra for more general case

        current_w1 <- sum(target_vals)

        if(abs(model[[constraintNode]] - current_w1) > 1e-6) print('Looks like the constraint is being violated\n'); 

        proposal_w2 <- rnorm(1, current_w2, scale)

        target_vals[1] <- 0.5*(current_w1 - proposal_w2)

        target_vals[2] <- 0.5*(current_w1 + proposal_w2)

        values(model, target) <<- target_vals

        accept <- decide(calculateDiff(model, calcNodes))

        if(accept) {

            copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)

        } else {

            copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)

        }

    },

    methods = list(

        reset = function() {}

    )

)

## Make an empty configuration by creating one and then removing the samplers

confMCMC <- configureMCMC(sillyModel)

confMCMC$getSamplers()

confMCMC$removeSamplers('a')

confMCMC$removeSamplers('b')

## alternatively we could do configureMCMC(sillyModel, nodes = NULL) to get an empty MCMC configuration directly (i.e. one with no samplers)

## now add one instance of our constrained sampler

confMCMC$addSampler(target = c('a','b'), type = RW_constrained, control = list(scale = 0.5, constraintNode = 'c'))

sillyMCMC <- buildMCMC(confMCMC)

cConstrainedExample <- compileNimble(sillyModel, sillyMCMC)

cConstrainedExample$sillyMCMC$run(10000)

samp <- as.matrix(cConstrainedExample$sillyMCMC$mvSamples)

plot(samp[,'a'])

plot(samp[,'b'])

plot(samp[,'a'], samp[, 'b'])

abline(0.7, -1, col = 'red') ## show the constraint

## end of example

Now one might want to use adaptive random walk instead of choosing a scale parameter as above.  To do that one could copy from samplers_RW in MCMC_samplers.R in the nimble/R directory of the source code.  One would then replace the first four lines of the run function with the first nine lines of the run function above.  (The error check could be removed).  A more general case of constraints could be done using linear algebra.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/2fd539f5-1d26-4cc4-9cbc-384c40d62f84%40googlegroups.com.

[Jose Jimenez]

Dear Perry:

Thank you very much for your example. It is also very informative about how Nimble works.

Thanks again,
Jose