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