Exercise 9.2A

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Exercise 9.2A

Greg
I'm sure I'm missing something simple, but I'm struggling with this problem.  When I make changes to FilconJags.R, I get the following message:

Error in jags.model("model.txt", data = dataList, inits = initsList, n.chains = nChains,  :
  RUNTIME ERROR:
Compilation error on line 10.
Subset out of range: kappa[2]


I went and looked at the answers, and they match, but of course, the answers are given in BRUGS, rather than JAGS.  I'm pretty sure the problem is with the "mean(kappa)" command, but I can't find the JAGS equivalent.  Any help would be appreciated.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Exercise 9.2A

John K. Kruschke
Administrator

I suspect it's the initialization, not the model specification. Try this: Don't "manually" initialize the parameters, instead let JAGS do it automatically. In the jags.model command, just comment out (or delete) the argument inits=...  Try the program. If it works, you know it was a problem with initialization. If you don't mind letting JAGS automatically initialize, you're done!



On Fri, Apr 12, 2013 at 3:58 PM, Greg [via Doing Bayesian Data Analysis] <[hidden email]> wrote:
I'm sure I'm missing something simple, but I'm struggling with this problem.  When I make changes to FilconJags.R, I get the following message:

Error in jags.model("model.txt", data = dataList, inits = initsList, n.chains = nChains,  :
  RUNTIME ERROR:
Compilation error on line 10.
Subset out of range: kappa[2]


I went and looked at the answers, and they match, but of course, the answers are given in BRUGS, rather than JAGS.  I'm pretty sure the problem is with the "mean(kappa)" command, but I can't find the JAGS equivalent.  Any help would be provided.


If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.x6.nabble.com/Exercise-9-2A-tp5000636.html
To start a new topic under Doing Bayesian Data Analysis, email [hidden email]
To unsubscribe from Doing Bayesian Data Analysis, click here.
NAML

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Exercise 9.2A

Greg Hayes
Here's what I'm getting now...

> source('~/Bayesian Analysis/ProgramsDoingBayesianDataAnalysis/FilconJagsExercise9_2a.R')
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Deleting model

Error in jags.model("model.txt", data = dataList, n.chains = nChains,  :
  RUNTIME ERROR:
Compilation error on line 10.
Subset out of range: kappa[2]

The following is the code changes I've made.



 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 ~ dgamma( Skappa , Rkappa )
......
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.

sqzData = .01+.98*dataList$z/dataList$N
mu = aggregate( sqzData , list(dataList$cond) , "mean" )[,"x"]
sd = aggregate( sqzData , list(dataList$cond) , "sd" )[,"x"]
kappa = mu*(1-mu)/sd^2 - 1
initsList = list( theta = sqzData , mu = mu , kappa = mean(kappa ))
*******

Changed the jagsModel command to:  
jagsModel = jags.model( "model.txt" , data=dataList ,  #inits=initsList ,
                        n.chains=nChains , n.adapt=adaptSteps )

Thanks!
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Exercise 9.2A

John K. Kruschke
Administrator

kappa is a single (scalar, non-vector) value, so you can't have any uses of kappa[condIdx]



On Fri, Apr 12, 2013 at 10:01 PM, Greg Hayes [via Doing Bayesian Data Analysis] <[hidden email]> wrote:
Here's what I'm getting now...

> source('~/Bayesian Analysis/ProgramsDoingBayesianDataAnalysis/FilconJagsExercise9_2a.R')
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Deleting model

Error in jags.model("model.txt", data = dataList, n.chains = nChains,  :
  RUNTIME ERROR:
Compilation error on line 10.
Subset out of range: kappa[2]

The following is the code changes I've made.



 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 ~ dgamma( Skappa , Rkappa )
......
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.

sqzData = .01+.98*dataList$z/dataList$N
mu = aggregate( sqzData , list(dataList$cond) , "mean" )[,"x"]
sd = aggregate( sqzData , list(dataList$cond) , "sd" )[,"x"]
kappa = mu*(1-mu)/sd^2 - 1
initsList = list( theta = sqzData , mu = mu , kappa = mean(kappa ))
*******

Changed the jagsModel command to:  
jagsModel = jags.model( "model.txt" , data=dataList ,  #inits=initsList ,
                        n.chains=nChains , n.adapt=adaptSteps )

Thanks!


If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.x6.nabble.com/Exercise-9-2A-tp5000636p5000638.html
To start a new topic under Doing Bayesian Data Analysis, email [hidden email]
To unsubscribe from Doing Bayesian Data Analysis, click here.
NAML

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Exercise 9.2A

Rothrock
There is also another place in FilconJAGS.R where there is a condIdx, down where you are extracting the parameter values from the mcmc chain returned by JAGS.

Change this line:

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

to this line:

kappa = rbind( kappa , mcmcChain[, "kappa"] )
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Exercise 9.2A

John K. Kruschke
Administrator
# Here's a complete program.
# The original lines with kappa are commented out, with the replacement nearby in the code.

graphics.off()
rm(list=ls(all=TRUE))
fileNameRoot="FilconJags-Ex9.2A-" # for constructing output filenames
if ( .Platform$OS.type != "windows" ) {
  windows <- function( ... ) X11( ... )
}
require(rjags)         # Kruschke, J. K. (2011). Doing Bayesian Data Analysis:
                       # A Tutorial with R and BUGS. Academic Press / Elsevier.
#------------------------------------------------------------------------------
# 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]] )T(0.001,0.999)
   }
   for ( condIdx in 1:nCond ) {
#      a[condIdx] <- mu[condIdx] * kappa[condIdx]
#      b[condIdx] <- (1-mu[condIdx]) * kappa[condIdx]
      a[condIdx] <- mu[condIdx] * kappa
      b[condIdx] <- (1-mu[condIdx]) * kappa
      # Hyperprior on mu and kappa:
      mu[condIdx] ~ dbeta( Amu , Bmu )
#      kappa[condIdx] ~ dgamma( Skappa , Rkappa )
   }
   kappa ~ dgamma( Skappa , Rkappa )
   # Constants for hyperprior:
   Amu <- 1
   Bmu <- 1
   Skappa <- pow(meanGamma,2)/pow(sdGamma,2)
   Rkappa <- meanGamma/pow(sdGamma,2)
   meanGamma <- 10
   sdGamma <- 10
}
" # close quote for modelstring
writeLines(modelstring,con="model.txt") # Write model to a file.

#------------------------------------------------------------------------------
# THE DATA.

# For each subject, specify the condition s/he was in,
# the number of trials s/he experienced, and the number correct.
# (These lines intentionally exceed the margins so that they don't take up
# excessive space on the printed page.)
cond = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4)
N = c(64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64)
z = c(45,63,58,64,58,63,51,60,59,47,63,61,60,51,59,45,61,59,60,58,63,56,63,64,64,60,64,62,49,64,64,58,64,52,64,64,64,62,64,61,59,59,55,62,51,58,55,54,59,57,58,60,54,42,59,57,59,53,53,42,59,57,29,36,51,64,60,54,54,38,61,60,61,60,62,55,38,43,58,60,44,44,32,56,43,36,38,48,32,40,40,34,45,42,41,32,48,36,29,37,53,55,50,47,46,44,50,56,58,42,58,54,57,54,51,49,52,51,49,51,46,46,42,49,46,56,42,53,55,51,55,49,53,55,40,46,56,47,54,54,42,34,35,41,48,46,39,55,30,49,27,51,41,36,45,41,53,32,43,33)
nSubj = length(cond)
nCond = length(unique(cond))

# Specify the data in a form that is compatible with BRugs model, as a list:
dataList = list(
 nCond = nCond ,
 nSubj = nSubj ,
 cond = cond ,
 N = N ,
 z = z
)

#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.

sqzData = .01+.98*dataList$z/dataList$N
mu = aggregate( sqzData , list(dataList$cond) , "mean" )[,"x"]
sd = aggregate( sqzData , list(dataList$cond) , "sd" )[,"x"]
# kappa = mu*(1-mu)/sd^2 - 1
kappa = mean( mu*(1-mu)/sd^2 - 1 )
initsList = list( theta = sqzData , mu = mu , kappa = kappa )

#------------------------------------------------------------------------------
# RUN THE CHAINS.

parameters = c("mu","kappa","theta","a","b")  # 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=5000           # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = 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=nPerChain , thin=thinSteps )
# resulting codaSamples object has these indices:
#   codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]

#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.

checkConvergence = F
if ( checkConvergence ) {
  show( summary( codaSamples ) )
  plot( codaSamples , ask=T ) 
  autocorr.plot( codaSamples , ask=T )
}

# Convert coda-object codaSamples to matrix object for easier handling.
# But note that this concatenates the different chains into one long chain.
# Result is mcmcChain[ stepIdx , paramIdx ]
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="") ] )
}
kappa = mcmcChain[,"kappa"]
save( mu , kappa , mcmcChain , file=paste(fileNameRoot,"MuKappa.Rdata",sep="") )
chainLength = NCOL(mu)

# Histograms of mu differences:
windows(10,2.75)
layout( matrix(1:3,nrow=1) )
source("plotPost.R")
histInfo = plotPost( mu[1,]-mu[2,] , xlab=expression(mu[1]-mu[2]) , main="" ,
          breaks=20 , compVal=0 , col="skyblue")
histInfo = plotPost( mu[3,]-mu[4,] , xlab=expression(mu[3]-mu[4]) , main="" ,
          breaks=20 , compVal=0 , col="skyblue")
histInfo = plotPost( (mu[1,]+mu[2,])/2 - (mu[3,]+mu[4,])/2 ,
          xlab=expression( (mu[1]+mu[2])/2 - (mu[3]+mu[4])/2 ) ,
          main="" , breaks=20 , compVal=0 , col="skyblue")
#savePlot(file=paste(fileNameRoot,"MuDiffs.eps",sep=""),type="eps")

# Scatterplot of mu,kappa in each condition:
windows()
muLim = c(.60,1) ; kappaLim = c( 2.0 , 40 ) ; 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" )
# points( mu[3,thindex] , kappa[3,thindex] , col="green" , pch="3" )
# points( mu[4,thindex] , kappa[4,thindex] , col="black" , pch="4" )
plot( mu[1,thindex] , kappa[thindex] , main=mainLab ,
      xlab=expression(mu[c]) , ylab=expression(kappa) , cex.lab=1.75 ,
      xlim=muLim , ylim=kappaLim , log="y" , col="red" , pch="1" )
points( mu[2,thindex] , kappa[thindex] , col="blue" , pch="2" )
points( mu[3,thindex] , kappa[thindex] , col="green" , pch="3" )
points( mu[4,thindex] , kappa[thindex] , col="black" , pch="4" )
#savePlot(file=paste(fileNameRoot,"Scatter.eps",sep=""),type="eps")


John K. Kruschke, Professor
Doing Bayesian Data Analysis
The book: http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/

The blog: http://doingbayesiandataanalysis.blogspot.com/





On Sun, May 5, 2013 at 1:56 PM, Rothrock [via Doing Bayesian Data Analysis] <[hidden email]> wrote:
There is also another place in FilconJAGS.R where there is a condIdx, down where you are extracting the parameter values from the mcmc chain returned by JAGS.

Change this line:

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

to this line:

kappa = rbind( kappa , mcmcChain[, "kappa"] )



If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.x6.nabble.com/Exercise-9-2A-tp5000636p5000661.html
To start a new topic under Doing Bayesian Data Analysis, email [hidden email]
To unsubscribe from Doing Bayesian Data Analysis, click here.
NAML

Loading...