BernMetropolisTemplate.R

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

BernMetropolisTemplate.R

Carl Dunham
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 one-row 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!
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: BernMetropolisTemplate.R

John K. Kruschke

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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: BernMetropolisTemplate.R

Carl Dunham
Thank you!

On Fri, Aug 3, 2012 at 7:35 AM, John K. Kruschke [via Doing Bayesian Data Analysis] <[hidden email]> wrote:

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.


If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.n6.nabble.com/BernMetropolisTemplate-R-tp5000491p5000492.html
To unsubscribe from BernMetropolisTemplate.R, click here.
NAML

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

Re: BernMetropolisTemplate.R

Hong Xu
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 "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:
##### CHANGE03 use inlining JAGS - textConnection
jagsModel = jags.model( textConnection(modelString) , 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.

# 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 )

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(D|theta).
# 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 * (1-theta)^(N-z)
  # 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*(1-meanTraj)/sdTraj^2) - 1 )
b = (1-meanTraj) * ( (meanTraj*(1-meanTraj)/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 single-component 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: BernMetropolisTemplate.R

John K. Kruschke
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,

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 "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:
##### CHANGE03 use inlining JAGS - textConnection
jagsModel = jags.model( textConnection(modelString) , 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.

# 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 )

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(D|theta).
# 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 * (1-theta)^(N-z)
  # 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*(1-meanTraj)/sdTraj^2) - 1 )
b = (1-meanTraj) * ( (meanTraj*(1-meanTraj)/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 single-component 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


If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.n6.nabble.com/BernMetropolisTemplate-R-tp5000491p5000602.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: BernMetropolisTemplate.R

Hong Xu
It works. Thanks a lot!

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

Re: BernMetropolisTemplate.R

Kevin McC
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

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

Re: BernMetropolisTemplate.R

John K. Kruschke
Administrator

Sorry for the histInfo error again. One of these days I'll have to do a systematic search-and-convert 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.

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




If you reply to this email, your message will be added to the discussion below:
http://doing-bayesian-data-analysis.12272.x6.nabble.com/BernMetropolisTemplate-R-tp5000491p5000672.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...