Is it possible to pass value to JAGS when running the R program?

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

Is it possible to pass value to JAGS when running the R program?

Frank
Hi John,

Thank you so much for the great book. I'm trying to go through the exercises. I have a question about Exercise 9.1. I converted BernBetaMuKappaJags.R to a function so that I can pass the number of coins into the program. But I got the error message - ""

Here is how I run the program:
    source("BernBetaMuKappaJags.R")
    BernBetaMuKappaJags(5)

Here is the code:

# change this code into function so that I can do Ex9.1 with different number of coins (5 or 50)
BernBetaMuKappaJags = function( numCoins = 50 ) {

rm(list = ls())
graphics.off()
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.

# Specify the model in JAGS language, but save it as a string in R:
modelString = "
# JAGS model specification begins ...
model {
    # Likelihood:
    for ( t in 1:nTrialTotal ) {
        y[t] ~ dbern( theta[ coin[ t ] ] )
    }
    # Prior:
    for ( j in 1:nCoins ) {
        theta[j] ~ dbeta( a , b )T(0.001,0.999)
    }
    a <- mu * kappa
    b <- ( 1.0 - mu ) * kappa
    mu ~ dbeta( Amu , Bmu )
    kappa ~ dgamma( Skappa , Rkappa )
    Amu <- 2.0
    Bmu <- 2.0
    Skappa <- pow(10,2)/pow(10,2)
    Rkappa <- 10/pow(10,2)
}
# ... JAGS model specification ends.
" # close quote to end modelString

# Write the modelString to a file, using R commands:
##### START: comment out the line below for inlining JAGS
# writeLines(modelString,con="model.txt")

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

##### CHANGE: from "ncoins = 50" to "ncoins = numCoins"
ncoins = numCoins ; nflipspercoin = 5
muAct = .7 ; kappaAct = 20
thetaAct = rbeta( ncoins ,muAct*kappaAct ,(1-muAct)*kappaAct )
z = rbinom( n=ncoins ,size=nflipspercoin ,prob=thetaAct )
N = rep( nflipspercoin , ncoins )

# Demo data for various figures in the book:
# N =  c( 5, 5, 5, 5, 5 ) # c( 10, 10, 10 )  # c( 15, 5 )
# z =  c( 1, 1, 1, 1, 5 ) # c(  1,  5,  9 )  # c(  3, 4 )
 
# Therapeutic touch data:
#  z = c(1,2,3,3,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,6,6,7,7,7,8)
#  N = rep(10,length(z))
# Convert z,N to vectors of individual flips.
# coin vector is index of coin for each flip.
# y vector is head or tail for each flip.
# For example,
#  coin = c( 1, 1, 2, 2, 2 )
#  y    = c( 1, 0, 0, 0, 1 )
# means that the first flip was of coin 1 and it was a head, the second flip
# was of coin 1 and it was a tail, the third flip was of coin 2 and it was a
# tail, etc.

coin = NULL ; y = NULL
for ( coinIdx in 1:length(N) ) {
    coin = c( coin , rep(coinIdx,N[coinIdx]) )
    y = c( y , rep(1,z[coinIdx]) , rep(0,N[coinIdx]-z[coinIdx]) )
}
nTrialTotal = length( y )
nCoins = length( unique( coin ) )

dataList = list(
    y = y ,
    coin = coin ,
    nTrialTotal = nTrialTotal ,
    nCoins = nCoins
)

#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.

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

parameters = c( "mu" , "kappa" , "theta" )   # The parameter(s) to be monitored.
adaptSteps = 1000              # Number of steps to "tune" the samplers.
burnInSteps = 5000            # 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( 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.

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 the posterior sample from JAGS:
muSample = mcmcChain[,"mu"]
kappaSample = mcmcChain[,"kappa"] # BRugs gets sample from JAGS
thetaSample = matrix( 0 , nrow=nCoins , ncol=nChains*nIter )
for ( coinIdx in 1:nCoins ) {
    nodeName = paste( "theta[" , coinIdx , "]" , sep="" )
    thetaSample[coinIdx,] = mcmcChain[,nodeName]
}

# Make a graph using R commands:
source("plotPost.R")
if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
  windows(3.2*3,2.5*(1+nCoins))
  layout( matrix( 1:(3*(nCoins+1)) , nrow=(nCoins+1) , byrow=T ) )
  par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1))
  nPtsToPlot = 500
  plotIdx = floor(seq(1,length(muSample),length=nPtsToPlot))
  kPltLim = signif( quantile( kappaSample , p=c(.01,.99) ) , 4 )
  plot( muSample[plotIdx] , kappaSample[plotIdx] , type="p" , ylim=kPltLim ,
        xlim=c(0,1) , xlab=expression(mu) , ylab=expression(kappa) , cex.lab=1.5 ,
        col="skyblue" )
  plotPost( muSample , xlab="mu" , xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
  plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , col="skyblue" ,
            showMode=T )
  for ( coinIdx in 1:nCoins ) {
      plotPost( thetaSample[coinIdx,] , xlab=paste("theta",coinIdx,sep="") ,
                xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , muSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=c(0,1) , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(mu) , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , kappaSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=kPltLim , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(kappa) , col="skyblue")
  }
} # end if ( nCoins <= ...

#
windows(width=7,height=5)
layout( matrix( 1:4 , nrow=2 , byrow=T ) )
par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1) )
plotPost( muSample , xlab="mu" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , HDItextPlace=.1 , col="skyblue")
plotPost( thetaSample[1,] , xlab="theta1" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
# Uncomment the following if using therapeutic touch data:
#plotPost( thetaSample[28,] , xlab="theta28" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
} # end of function


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

Re: Is it possible to pass value to JAGS when running the R program?

John K. Kruschke
Administrator


What is the error message?



On Tue, Mar 26, 2013 at 10:03 PM, Frank [via Doing Bayesian Data Analysis] <[hidden email]> wrote:
Hi John,

Thank you so much for the great book. I'm trying to go through the exercises. I have a question about Exercise 9.1. I converted BernBetaMuKappaJags.R to a function so that I can pass the number of coins into the program. But I got the error message - ""

Here is how I run the program:
    source("BernBetaMuKappaJags.R")
    BernBetaMuKappaJags(5)

Here is the code:

# change this code into function so that I can do Ex9.1 with different number of coins (5 or 50)
BernBetaMuKappaJags = function( numCoins = 50 ) {

rm(list = ls())
graphics.off()
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.

# Specify the model in JAGS language, but save it as a string in R:
modelString = "
# JAGS model specification begins ...
model {
    # Likelihood:
    for ( t in 1:nTrialTotal ) {
        y[t] ~ dbern( theta[ coin[ t ] ] )
    }
    # Prior:
    for ( j in 1:nCoins ) {
        theta[j] ~ dbeta( a , b )T(0.001,0.999)
    }
    a <- mu * kappa
    b <- ( 1.0 - mu ) * kappa
    mu ~ dbeta( Amu , Bmu )
    kappa ~ dgamma( Skappa , Rkappa )
    Amu <- 2.0
    Bmu <- 2.0
    Skappa <- pow(10,2)/pow(10,2)
    Rkappa <- 10/pow(10,2)
}
# ... JAGS model specification ends.
" # close quote to end modelString

# Write the modelString to a file, using R commands:
##### START: comment out the line below for inlining JAGS
# writeLines(modelString,con="model.txt")

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

##### CHANGE: from "ncoins = 50" to "ncoins = numCoins"
ncoins = numCoins ; nflipspercoin = 5
muAct = .7 ; kappaAct = 20
thetaAct = rbeta( ncoins ,muAct*kappaAct ,(1-muAct)*kappaAct )
z = rbinom( n=ncoins ,size=nflipspercoin ,prob=thetaAct )
N = rep( nflipspercoin , ncoins )

# Demo data for various figures in the book:
# N =  c( 5, 5, 5, 5, 5 ) # c( 10, 10, 10 )  # c( 15, 5 )
# z =  c( 1, 1, 1, 1, 5 ) # c(  1,  5,  9 )  # c(  3, 4 )
 
# Therapeutic touch data:
#  z = c(1,2,3,3,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,6,6,7,7,7,8)
#  N = rep(10,length(z))
# Convert z,N to vectors of individual flips.
# coin vector is index of coin for each flip.
# y vector is head or tail for each flip.
# For example,
#  coin = c( 1, 1, 2, 2, 2 )
#  y    = c( 1, 0, 0, 0, 1 )
# means that the first flip was of coin 1 and it was a head, the second flip
# was of coin 1 and it was a tail, the third flip was of coin 2 and it was a
# tail, etc.

coin = NULL ; y = NULL
for ( coinIdx in 1:length(N) ) {
    coin = c( coin , rep(coinIdx,N[coinIdx]) )
    y = c( y , rep(1,z[coinIdx]) , rep(0,N[coinIdx]-z[coinIdx]) )
}
nTrialTotal = length( y )
nCoins = length( unique( coin ) )

dataList = list(
    y = y ,
    coin = coin ,
    nTrialTotal = nTrialTotal ,
    nCoins = nCoins
)

#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.

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

parameters = c( "mu" , "kappa" , "theta" )   # The parameter(s) to be monitored.
adaptSteps = 1000              # Number of steps to "tune" the samplers.
burnInSteps = 5000            # 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( 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.

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 the posterior sample from JAGS:
muSample = mcmcChain[,"mu"]
kappaSample = mcmcChain[,"kappa"] # BRugs gets sample from JAGS
thetaSample = matrix( 0 , nrow=nCoins , ncol=nChains*nIter )
for ( coinIdx in 1:nCoins ) {
    nodeName = paste( "theta[" , coinIdx , "]" , sep="" )
    thetaSample[coinIdx,] = mcmcChain[,nodeName]
}

# Make a graph using R commands:
source("plotPost.R")
if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
  windows(3.2*3,2.5*(1+nCoins))
  layout( matrix( 1:(3*(nCoins+1)) , nrow=(nCoins+1) , byrow=T ) )
  par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1))
  nPtsToPlot = 500
  plotIdx = floor(seq(1,length(muSample),length=nPtsToPlot))
  kPltLim = signif( quantile( kappaSample , p=c(.01,.99) ) , 4 )
  plot( muSample[plotIdx] , kappaSample[plotIdx] , type="p" , ylim=kPltLim ,
        xlim=c(0,1) , xlab=expression(mu) , ylab=expression(kappa) , cex.lab=1.5 ,
        col="skyblue" )
  plotPost( muSample , xlab="mu" , xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
  plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , col="skyblue" ,
            showMode=T )
  for ( coinIdx in 1:nCoins ) {
      plotPost( thetaSample[coinIdx,] , xlab=paste("theta",coinIdx,sep="") ,
                xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , muSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=c(0,1) , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(mu) , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , kappaSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=kPltLim , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(kappa) , col="skyblue")
  }
} # end if ( nCoins <= ...

#
windows(width=7,height=5)
layout( matrix( 1:4 , nrow=2 , byrow=T ) )
par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1) )
plotPost( muSample , xlab="mu" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , HDItextPlace=.1 , col="skyblue")
plotPost( thetaSample[1,] , xlab="theta1" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
# Uncomment the following if using therapeutic touch data:
#plotPost( thetaSample[28,] , xlab="theta28" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
} # end of function


thanks,
Frank


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: Is it possible to pass value to JAGS when running the R program?

John K. Kruschke
Administrator
In reply to this post by Frank

I just noticed that you still have
rm(list = ls())
at the beginning of the function. It's basically erasing itself.


On Tue, Mar 26, 2013 at 10:03 PM, Frank [via Doing Bayesian Data Analysis] <[hidden email]> wrote:
Hi John,

Thank you so much for the great book. I'm trying to go through the exercises. I have a question about Exercise 9.1. I converted BernBetaMuKappaJags.R to a function so that I can pass the number of coins into the program. But I got the error message - ""

Here is how I run the program:
    source("BernBetaMuKappaJags.R")
    BernBetaMuKappaJags(5)

Here is the code:

# change this code into function so that I can do Ex9.1 with different number of coins (5 or 50)
BernBetaMuKappaJags = function( numCoins = 50 ) {

rm(list = ls())
graphics.off()
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.

# Specify the model in JAGS language, but save it as a string in R:
modelString = "
# JAGS model specification begins ...
model {
    # Likelihood:
    for ( t in 1:nTrialTotal ) {
        y[t] ~ dbern( theta[ coin[ t ] ] )
    }
    # Prior:
    for ( j in 1:nCoins ) {
        theta[j] ~ dbeta( a , b )T(0.001,0.999)
    }
    a <- mu * kappa
    b <- ( 1.0 - mu ) * kappa
    mu ~ dbeta( Amu , Bmu )
    kappa ~ dgamma( Skappa , Rkappa )
    Amu <- 2.0
    Bmu <- 2.0
    Skappa <- pow(10,2)/pow(10,2)
    Rkappa <- 10/pow(10,2)
}
# ... JAGS model specification ends.
" # close quote to end modelString

# Write the modelString to a file, using R commands:
##### START: comment out the line below for inlining JAGS
# writeLines(modelString,con="model.txt")

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

##### CHANGE: from "ncoins = 50" to "ncoins = numCoins"
ncoins = numCoins ; nflipspercoin = 5
muAct = .7 ; kappaAct = 20
thetaAct = rbeta( ncoins ,muAct*kappaAct ,(1-muAct)*kappaAct )
z = rbinom( n=ncoins ,size=nflipspercoin ,prob=thetaAct )
N = rep( nflipspercoin , ncoins )

# Demo data for various figures in the book:
# N =  c( 5, 5, 5, 5, 5 ) # c( 10, 10, 10 )  # c( 15, 5 )
# z =  c( 1, 1, 1, 1, 5 ) # c(  1,  5,  9 )  # c(  3, 4 )
 
# Therapeutic touch data:
#  z = c(1,2,3,3,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,6,6,7,7,7,8)
#  N = rep(10,length(z))
# Convert z,N to vectors of individual flips.
# coin vector is index of coin for each flip.
# y vector is head or tail for each flip.
# For example,
#  coin = c( 1, 1, 2, 2, 2 )
#  y    = c( 1, 0, 0, 0, 1 )
# means that the first flip was of coin 1 and it was a head, the second flip
# was of coin 1 and it was a tail, the third flip was of coin 2 and it was a
# tail, etc.

coin = NULL ; y = NULL
for ( coinIdx in 1:length(N) ) {
    coin = c( coin , rep(coinIdx,N[coinIdx]) )
    y = c( y , rep(1,z[coinIdx]) , rep(0,N[coinIdx]-z[coinIdx]) )
}
nTrialTotal = length( y )
nCoins = length( unique( coin ) )

dataList = list(
    y = y ,
    coin = coin ,
    nTrialTotal = nTrialTotal ,
    nCoins = nCoins
)

#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.

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

parameters = c( "mu" , "kappa" , "theta" )   # The parameter(s) to be monitored.
adaptSteps = 1000              # Number of steps to "tune" the samplers.
burnInSteps = 5000            # 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( 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.

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 the posterior sample from JAGS:
muSample = mcmcChain[,"mu"]
kappaSample = mcmcChain[,"kappa"] # BRugs gets sample from JAGS
thetaSample = matrix( 0 , nrow=nCoins , ncol=nChains*nIter )
for ( coinIdx in 1:nCoins ) {
    nodeName = paste( "theta[" , coinIdx , "]" , sep="" )
    thetaSample[coinIdx,] = mcmcChain[,nodeName]
}

# Make a graph using R commands:
source("plotPost.R")
if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
  windows(3.2*3,2.5*(1+nCoins))
  layout( matrix( 1:(3*(nCoins+1)) , nrow=(nCoins+1) , byrow=T ) )
  par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1))
  nPtsToPlot = 500
  plotIdx = floor(seq(1,length(muSample),length=nPtsToPlot))
  kPltLim = signif( quantile( kappaSample , p=c(.01,.99) ) , 4 )
  plot( muSample[plotIdx] , kappaSample[plotIdx] , type="p" , ylim=kPltLim ,
        xlim=c(0,1) , xlab=expression(mu) , ylab=expression(kappa) , cex.lab=1.5 ,
        col="skyblue" )
  plotPost( muSample , xlab="mu" , xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
  plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , col="skyblue" ,
            showMode=T )
  for ( coinIdx in 1:nCoins ) {
      plotPost( thetaSample[coinIdx,] , xlab=paste("theta",coinIdx,sep="") ,
                xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , muSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=c(0,1) , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(mu) , col="skyblue")
      plot( thetaSample[coinIdx,plotIdx] , kappaSample[plotIdx] , type="p" ,
            xlim=c(0,1) , ylim=kPltLim , cex.lab=1.5 ,
            xlab=bquote(theta[.(coinIdx)]) , ylab=expression(kappa) , col="skyblue")
  }
} # end if ( nCoins <= ...

#
windows(width=7,height=5)
layout( matrix( 1:4 , nrow=2 , byrow=T ) )
par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1) )
plotPost( muSample , xlab="mu" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , HDItextPlace=.1 , col="skyblue")
plotPost( thetaSample[1,] , xlab="theta1" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
# Uncomment the following if using therapeutic touch data:
#plotPost( thetaSample[28,] , xlab="theta28" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
} # end of function


thanks,
Frank


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: Is it possible to pass value to JAGS when running the R program?

Frank
In reply to this post by John K. Kruschke
I'm sorry I forgot paste the error message. Here it is:
  ## Error: object ’numCoins’ not found
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Is it possible to pass value to JAGS when running the R program?

Frank
In reply to this post by John K. Kruschke
Thank you so much for your quick response, John! I just comment out the line of "rm(list = ls())". It works now.

thanks,
Frank
Loading...