Great book!
I am trying to run BernMetropolisTemplate.R, as is, to start playing with it, and I get this error: Error in histInfo$density : $ operator is invalid for atomic vectors In looking at plotPost.R, it seems like the returned value is different than what the calling code is looking for. It's a onerow matrix, with no column named "density", even. I can get by using another value, it's only using it for positioning the text, but it seemed worth mentioning. Thanks! 
Thank you very much for alerting me to this problem! I recently changed the output of the plotPost() function, knowing that it might cause a problem somewhere, but in all the dozens of programs I couldn't find any conflicts. Sorry you had to find it for me! I have now updated BernMetropolisTemplate.R (and the zip file containing all the programs) at the program web site: http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/Programs/?M=D Thank you again for notifying me, and again apologies for the glitch. 
Thank you!
On Fri, Aug 3, 2012 at 7:35 AM, John K. Kruschke [via Doing Bayesian Data Analysis] <[hidden email]> wrote:

In reply to this post by John K. Kruschke
Dr. Kruschke,
Great book! I'm learning it by myself. I have no prior training in statistics (not contaminated by NHST yet :). I still got the error of "Error in histInfo$density : $ operator is invalid for atomic vectors" for Exercise 7.5 code. I'm using JAGS instead of BUGS. See my code below: require(rjags) # Kruschke, J. K. (2011). Doing Bayesian Data Analysis: # A Tutorial with R and BUGS. Academic Press / Elsevier. # # THE MODEL. # Specify the model in JAGS language, but save it as a string in R: modelString = " # JAGS model specification begins ... model { # Likelihood: for ( i in 1:nFlips ) { y[i] ~ dbern( theta ) } # Prior distribution: ##### CHANGE01 start theta ~ dunif( priorA , priorB ) priorA < 0 priorB < 0.4 ##### CHANGE01 end } # ... JAGS model specification ends. " # close quote to end modelString # Write the modelString to a file, using R commands: ##### CHANG02 comment out the line below for inlining JAGS # writeLines(modelString,con="model.txt") # # THE DATA. # Specify the data in R, using a list format compatible with JAGS: dataList = list( nFlips = 14 , y = c( 1,1,1,1,1,1,1,1,1,1,1,0,0,0 ) ) # # INTIALIZE THE CHAIN. # Can be done automatically in jags.model() by commenting out inits argument. # Otherwise could be established as: # initsList = list( theta = sum(dataList$y)/length(dataList$y) ) # # RUN THE CHAINS. parameters = c( "theta" ) # The parameter(s) to be monitored. adaptSteps = 500 # Number of steps to "tune" the samplers. burnInSteps = 1000 # Number of steps to "burnin" 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: ##### CHANGE03 use inlining JAGS  textConnection jagsModel = jags.model( textConnection(modelString) , data=dataList , # inits=initsList , n.chains=nChains , n.adapt=adaptSteps ) # Burnin: 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. # Convert codaobject 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 ) thetaSample = mcmcChain # Make a graph using R commands: # dev.new(10,6) layout( matrix( c(1,2) , nrow=1 ) ) plot( thetaSample[1:500] , 1:length(thetaSample[1:500]) , type="o" , xlim=c(0,1) , xlab=bquote(theta) , ylab="Position in Chain" , cex.lab=1.25 , main="JAGS Results" ) source("plotPost.R") histInfo = plotPost( thetaSample , xlim=c(0,1) , col="skyblue" ) acceptedTraj = thetaSample meanTraj = mean( thetaSample ) sdTraj = apply( thetaSample, 2, sd ) myData = dataList$y # Define the Bernoulli likelihood function, p(Dtheta). # The argument theta could be a vector, not just a scalar. likelihood = function( theta , data ) { z = sum( data == 1 ) N = length( data ) pDataGivenTheta = theta^z * (1theta)^(Nz) # The theta values passed into this function are generated at random, # and therefore might be inadvertently greater than 1 or less than 0. # The likelihood for theta > 1 or for theta < 0 is zero: pDataGivenTheta[ theta > 1  theta < 0 ] = 0 return( pDataGivenTheta ) } prior = function( theta ) { prior = dunif( theta , min=0 , max=0.4 ) # uniform density # The theta values passed into this function are generated at random, # and therefore might be inadvertently greater than 1 or less than 0. # The prior for theta > 1 or for theta < 0 is zero: prior[ theta > 1  theta < 0 ] = 0 return( prior ) } # Compute a,b parameters for beta distribution that has the same mean # and stdev as the sample from the posterior. This is a useful choice # when the likelihood function is Bernoulli. a = meanTraj * ( (meanTraj*(1meanTraj)/sdTraj^2)  1 ) b = (1meanTraj) * ( (meanTraj*(1meanTraj)/sdTraj^2)  1 ) # For every theta value in the posterior sample, compute # dbeta(theta,a,b) / likelihood(theta)*prior(theta) # This computation assumes that likelihood and prior are proper densities, # i.e., not just relative probabilities. This computation also assumes that # the likelihood and prior functions were defined to accept a vector argument, # not just a singlecomponent scalar argument. wtdEvid = dbeta( acceptedTraj , a , b ) / ( likelihood( acceptedTraj , myData ) * prior( acceptedTraj ) ) pData = 1 / mean( wtdEvid ) # Display p(D) in the graph if ( meanTraj > .5 ) { xpos = 0.0 ; xadj = 0.0 } else { xpos = 1.0 ; xadj = 1.0 } densMax = max( histInfo$density ) text( xpos , 0.9*densMax , bquote( p(D)==.( signif(pData,3) ) ) , adj=c(xadj,0) , cex=1.5 ) thanks, Hong 
Administrator

Ah, that error is caused by a change in the plotPost() function. It used to return aspects of the histogram, including the probability densities at each bin. The new version of plotPost() does not return that information, and therefore histInfo$density is not defined. Instead of densMax = max( histInfo$density ) text( xpos , 0.9*densMax , bquote( p(D)==.( signif(pData,3) ) ) , adj=c(xadj,0) , cex=1.5 ) perhaps try this: text( xpos , 0 , bquote( p(D)==.( signif(pData,3) ) ) , adj=c(xadj,2) , cex=1.5 ) On Wed, Feb 27, 2013 at 10:44 PM, Hong Xu [via Doing Bayesian Data Analysis] <[hidden email]> wrote: Dr. Kruschke, 
It works. Thanks a lot!
Hong 
In reply to this post by John K. Kruschke
Just an FYI I was working through Ch. 13, example problem 13.3 and have also come across this issue.
If I run FilconJagsPower.R, I get the error: Error in histInfo$density : $ operator is invalid for atomic vectors I assume its in reference to line 106: text( mean(HDIlim) , .25*max(histInfo$density) , adj=c(.5,0) , cex=1.25 , bquote( "HDI width = " * .(signif(HDIlim[2]HDIlim[1],3)) ) ) If I understand this right, you are telling R to put a label on the histogram at a quarter height of the Highest density. I am sure there must be a more elegant way to change the code, but I just changed it to: text( mean(HDIlim) , 4 , adj=c(.5,0) , cex=1.25 , bquote( "HDI width = " * .(signif(HDIlim[2]HDIlim[1],3)) ) ) P.S. I love your book 
Administrator

Sorry for the histInfo error again. One of these days I'll have to do a systematic searchandconvert for every occurrence of histInfo. Meanwhile, thanks for reporting the error and figuring out a solution! On Fri, Jun 28, 2013 at 6:15 PM, Kevin McC [via Doing Bayesian Data Analysis] <[hidden email]> wrote: Just an FYI I was working through Ch. 13, example problem 13.3 and have also come across this issue. 
Free forum by Nabble  Edit this page 