0

Kruschke の Doing Bayesian Data Analysis の第 9 章から推定して、JAGS で階層分析を実行しようとしています。 2 つのミント、および各ミントからのコインの平均バイアスの推定値 (ミント バイアス: オメガ)。各ミントのバイアス、カッパの変動性を一定に保ちました。問題は、2 番目の造幣局から事後推定を取得できないことです。それは、前のものをサンプリングしているだけのようです。2 番目のミントの事後推定値を生成するために、モデル文字列テキストを修正する方法を知っている人はいますか (以下のステップ 3 を参照)。

以下の分析用スクリプト全体

library(rjags)
library(runjags)
library(coda)


############### 1. Generate the data 

flips <- c(sample(c(rep(1,3), rep(0,9))), # coin 1, mint 1, 12 flips total
           sample(c(rep(1,1), rep(0,4))), # coin 2, mint 1, 5 flips total
           sample(c(rep(1,10), rep(0,5))), # coin 1, mint 2, 15 flips
           sample(c(rep(1,17), rep(0,6)))) # coin 2, mint 2, 23 flips

coins <- factor(c(rep(1,12), rep(2,5), rep(3, 15), rep(4, 23)))

mints <- factor(c(rep(1,17), rep(2,38)))

nFlips <- length(flips) 
nCoins <- length(unique(coins))
nMints <- length(unique(mints))


#################### 2. Pass data into a list 

dataList <- list(
  flips = flips,
  coins = coins,
  mints = mints,
  nFlips = nFlips,
  nCoins = nCoins,
  nMints = nMints)


################### 3. specify and save the model 

modelString <- "
model{

      # start with nested likelihood function
      for (i in 1:nFlips) {

              flips[i] ~ dbern(theta[coins[i]])
      } 

      # next the prior on theta
      for (coins in 1:nCoins) {

              theta[coins] ~ dbeta(omega[mints[coins]]*(kappa - 2) + 1, (1 - omega[mints[coins]])*(kappa - 2) + 1) 
      }

      # next we specify the prior for the higher-level parameters on the mint, omega and kappa
      for (mints in 1:nMints) {

              omega[mints] ~ dbeta(2,2)

      }

      kappa <- 5
}
"


writeLines(modelString, "tempModelHier4CoinTwoMint.txt")

############################### Step 4: Initialise Chains 

initsList <- list(theta1 = mean(flips[coins==1]),
                  theta2 = mean(flips[coins==2]),
                  theta3 = mean(flips[coins==3]),
                  theta4 = mean(flips[coins==4]),
                  omega1 = mean(c(mean(flips[coins==1]),
                                 mean(flips[coins==2]))),
                  omega2 = mean(c(mean(flips[coins==3]),
                                 mean(flips[coins==4]))))

initsList


############################### Step 5: Generate Chains 

runJagsOut <- run.jags(method = "simple",
                       model = "tempModelHier4CoinTwoMint.txt",
                       monitor = c("theta[1]", "theta[2]", "theta[3]", "theta[4]", "omega[1]", "omega[2]"),
                       data = dataList,
                       inits = initsList,
                       n.chains = 1,
                       adapt = 500,
                       burnin = 1000,
                       sample = 50000,
                       thin = 1,
                       summarise = FALSE,
                       plots = FALSE)



############################### Step 6: Convert to Coda Object 

codaSamples <- as.mcmc.list(runJagsOut)

head(codaSamples)


############################### Step 7: Make Graphs 

df <- data.frame(as.matrix(codaSamples))

theta1 <- ggplot(df, aes(x = df$theta.1.)) + geom_density()
theta2 <- ggplot(df, aes(x = df$theta.2.)) + geom_density()
theta3 <- ggplot(df, aes(x = df$theta.3.)) + geom_density()
theta4 <- ggplot(df, aes(x = df$theta.4.)) + geom_density()
omega1 <- ggplot(df, aes(x = df$omega.1.)) + geom_density()
omega2 <- ggplot(df, aes(x = df$omega.2.)) + geom_density()

require(gridExtra)

ggsave("coinsAndMintsHier/hierPropFourCoinsTwoMints.pdf", grid.arrange(theta1, theta2, theta3, theta4, omega1, omega2, ncol = 2), device = "pdf", height = 30, width = 10, units = "cm")
4

1 に答える 1

1

問題は、事前設定を on に設定するときに、コインのミントにインデックスを付けようとしていた方法でしたtheta。thetaこの場合、 ではなく4 しかありませんnFlips。ネストされたインデックスは、各コインが属するミントのベクトルではなく、データ ベクトルmints[coins]にアクセスしていました。mints以下に修正版を作成しました。各コインがどのミントに属しているかを示すベクトルの明示的な構成に注意してください。モデル仕様でも、各 for ループ インデックスには、データ名とは異なる独自のインデックス名があることに注意してください。

graphics.off() # This closes all of R's graphics windows.
rm(list=ls())  # Careful! This clears all of R's memory!

library(runjags)
library(coda)

#library(rjags)

############### 1. Generate the data 

flips <- c(sample(c(rep(1,3), rep(0,9))), # coin 1, mint 1, 12 flips total
           sample(c(rep(1,1), rep(0,4))), # coin 2, mint 1, 5 flips total
           sample(c(rep(1,10), rep(0,5))), # coin 1, mint 2, 15 flips
           sample(c(rep(1,17), rep(0,6)))) # coin 2, mint 2, 23 flips

# NOTE: I got rid of `factor` because it was unneeded and got in the way
coins <- c(rep(1,12), rep(2,5), rep(3, 15), rep(4, 23))

# NOTE: I got rid of `factor` because it was unneeded and got in the way
mints <- c(rep(1,17), rep(2,38))

nFlips <- length(flips) 
nCoins <- length(unique(coins))
nMints <- length(unique(mints))

# NEW: Create vector that specifies the mint of each coin. There must be a     more 
# elegant way to do this, but here is a logical brute-force approach. This
# assumes that coins are consecutively numbered from 1 to nCoins.
mintOfCoin = NULL
for ( cIdx in 1:nCoins ) {
  mintOfCoin = c( mintOfCoin , unique(mints[coins==cIdx]) )
}

#################### 2. Pass data into a list 

dataList <- list(
  flips = flips,
  coins = coins,
  mints = mints,
  nFlips = nFlips,
  nCoins = nCoins,
  nMints = nMints,
  mintOfCoin = mintOfCoin # NOTE
  )


################### 3. specify and save the model 

modelString <- "
model{
  # start with nested likelihood function
  for (fIdx in 1:nFlips) {
    flips[fIdx] ~ dbern( theta[coins[fIdx]] )
  } 
  # next the prior on theta
  # NOTE: Here we use the mintOfCoin index.
  for (cIdx in 1:nCoins) {
    theta[cIdx] ~ dbeta( omega[mintOfCoin[cIdx]]*(kappa - 2) + 1 ,
                          ( 1 - omega[mintOfCoin[cIdx]])*(kappa - 2) + 1 ) 
  }
  # next we specify the prior for the higher-level parameters on the mint, 
  # omega and kappa
  # NOTE: I changed the name of the mint index so it doesn't conflict with 
  # mints data vector.
  for (mIdx in 1:nMints) {
    omega[mIdx] ~ dbeta(2,2)
  }
  kappa <- 5
}
"


writeLines(modelString, "tempModelHier4CoinTwoMint.txt")

############################### Step 4: Initialise Chains 

initsList <- list(theta1 = mean(flips[coins==1]),
                  theta2 = mean(flips[coins==2]),
                  theta3 = mean(flips[coins==3]),
                  theta4 = mean(flips[coins==4]),
                  omega1 = mean(c(mean(flips[coins==1]),
                                  mean(flips[coins==2]))),
                  omega2 = mean(c(mean(flips[coins==3]),
                                  mean(flips[coins==4]))))

initsList


############################### Step 5: Generate Chains 

runJagsOut <- run.jags(method = "parallel",
                       model = "tempModelHier4CoinTwoMint.txt",
                       # NOTE: theta and omega are vectors:
                       monitor = c( "theta", "omega" , "kappa" ),
                       data = dataList,
                       #inits = initsList, # NOTE: Let JAGS initialize.
                       n.chains = 4, # NOTE: Not only 1 chain.
                       adapt = 500,
                       burnin = 1000,
                       sample = 10000,
                       thin = 1,
                       summarise = FALSE,
                       plots = FALSE)



############################### Step 6: Convert to Coda Object 

codaSamples <- as.mcmc.list(runJagsOut)

head(codaSamples)

########################################
## NOTE: Important step --- Check MCMC diagnostics 

# Display diagnostics of chain, for specified parameters:
source("DBDA2E-utilities.R") # For function diagMCMC()
parameterNames = varnames(codaSamples) # from coda package
for ( parName in parameterNames ) {
  diagMCMC( codaObject=codaSamples , parName=parName )
}



############################### Step 7: Make Graphs 
# ...
于 2016-12-14T14:39:51.033 に答える