**I converted all the BUGS programs to JAGS versions back on 01 January 2012. Please see the previous post for information about JAGS.**

A reader reports this problem:

I [the reader] have received this error several times:When I try running the program from Exercise 12.3 with BUGS, I also get the error, although I did not get any error a few years ago when I created the program and reported its results in the solutions manual. I have tried to solve the problem by arbitrarily censoring/truncating the distributions in BUGS (as I did in some other programs in the book), but without success in this application.

"update error for node <theta[66]> algorithm conjugate beta updater error argument two has too small a value "

... In particular, I [the reader] noticed for exercise 12.3, if I remove the 0's from the data, I do not get the error.

I quickly re-wrote the BUGS+BRugs program for JAGS+rjags, going step-by-step through the guidelines from this previous post. The resulting program,

*using the identical model statement*, runs fine in JAGS and produces results identical to those in the solutions manual. Here is the complete program and results:

graphics.off()

rm(list=ls(all=TRUE))

fileNameRoot="Exercise12.3.A.JAGS" # for constructing output filenames

require(rjags)

if ( .Platform$OS.type != "windows" ) {

windows <- function( ... ) X11( ... )

}

#------------------------------------------------------------------------------

# THE MODEL.

modelstring = "

model {

for ( subjIdx in 1:nSubj ) {

# Likelihood:

z[subjIdx] ~ dbin( theta[subjIdx] , N[subjIdx] )

# Prior on theta: Notice nested indexing.

theta[subjIdx] ~ dbeta( a[cond[subjIdx]] , b[cond[subjIdx]] )

}

for ( condIdx in 1:nCond ) {

a[condIdx] <- mu[condIdx] * kappa[condIdx]

b[condIdx] <- (1-mu[condIdx]) * kappa[condIdx]

# Hyperprior on mu and kappa:

mu[condIdx] ~ dbeta( Amu , Bmu )

kappa[condIdx] ~ dgamma( Skappa , Rkappa )

}

# Constants for hyperprior:

Amu <- 0.35 * 2

Bmu <- 0.65 * 2

Skappa <- pow(meanGamma,2)/pow(sdGamma,2)

Rkappa <- meanGamma/pow(sdGamma,2)

meanGamma <- 10

sdGamma <- 10

}

" # close quote for modelstring

# Write model to a file:

writeLines(modelstring,con="model.txt")

#------------------------------------------------------------------------------

# THE DATA.

nCorrOfSubj = c( 8, 4, 6, 3, 1, 4, 4, 6, 4, 2, 2, 1, 1, 4, 3 ,3 ,2 ,6, 3, 4,

2, 1, 1, 3, 2, 7, 2, 1, 3, 1, 0, 2, 4, 2, 3, 3, 0, 1, 2, 2 )

CondOfSubj = c( rep(1,20) , rep(2,20) )

nTrlOfSubj = rep(10,40)

nSubj = length(CondOfSubj)

nCond = length(unique(CondOfSubj))

# Specify the data in a form that is compatible with BRugs model, as a list:

dataList = list(

nCond = nCond ,

nSubj = nSubj ,

cond = CondOfSubj ,

N = nTrlOfSubj ,

z = nCorrOfSubj

)

#------------------------------------------------------------------------------

# INTIALIZE THE CHAINS.

#------------------------------------------------------------------------------

# RUN THE CHAINS.

parameters = c("mu","kappa","theta") # The parameter(s) to be monitored.

adaptSteps = 500 # Number of steps to "tune" the samplers.

burnInSteps = 1000 # Number of steps to "burn-in" the samplers.

nChains = 3 # Number of chains to run.

numSavedSteps=50000 # Total number of steps in chains to save.

thinSteps=1 # Number of steps to "thin" (1=keep every step).

nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.

# Create, initialize, and adapt the model:

jagsModel = jags.model( "model.txt" , data=dataList , # inits=initsList ,

n.chains=nChains , n.adapt=adaptSteps )

# Burn-in:

cat( "Burning in the MCMC chain...\n" )

update( jagsModel , n.iter=burnInSteps )

# The saved MCMC chain:

cat( "Sampling final MCMC chain...\n" )

codaSamples = coda.samples( jagsModel , variable.names=parameters ,

n.iter=nIter , thin=thinSteps )

# resulting codaSamples object has these indices:

# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]

#------------------------------------------------------------------------------

# EXAMINE THE RESULTS.

source("plotPost.R")

mcmcChain = as.matrix( codaSamples )

# Extract parameter values and save them.

mu = NULL

kappa = NULL

for ( condIdx in 1:nCond ) {

mu = rbind( mu , mcmcChain[, paste("mu[",condIdx,"]",sep="") ] )

kappa = rbind( kappa , mcmcChain[, paste("kappa[",condIdx,"]",sep="") ] )

}

save( mu , kappa , file=paste(fileNameRoot,"MuKappa.Rdata",sep="") )

chainLength = NCOL(mu)

# Histograms of mu differences:

windows(width=7,height=7)

plotPost( mu[1,]-mu[2,] , xlab=expression(mu[1]-mu[2]) , main="" ,

breaks=20 , compVal=0 , ROPE=c(-.1,.1) )

savePlot( file=paste(fileNameRoot,"MuDiffs.jpg",sep="") , type="jpg" )

# Scatterplot of mu,kappa in each condition:

windows(width=7,height=7)

muLim = c(min(mu),max(mu))

kappaLim = c(min(kappa),max(kappa))

mainLab="Posterior"

thindex = round( seq( 1 , chainLength , len=300 ) )

plot( mu[1,thindex] , kappa[1,thindex] , main=mainLab ,

xlab=expression(mu[c]) , ylab=expression(kappa[c]) , cex.lab=1.75 ,

xlim=muLim , ylim=kappaLim , log="y" , col="red" , pch="1" )

points( mu[2,thindex] , kappa[2,thindex] , col="blue" , pch="2" )

savePlot( file=paste(fileNameRoot,"Scatter.jpg",sep="") , type="jpg" )

What about Stan?

ReplyDeleteI'd love it if people posted some translations of DBDA programs into Stan!

ReplyDelete