JAGS/rjags에서 여러 그룹에 대한 베이 즈 독립 매개 변수를 별도로 할당합니다.

Question

Kruschke의 Doing Bayesian Data Analysis (제 9 장)에서 외삽하여 JAGS에서 계층 적 분석을 수행하려고합니다. 4 개의 헤드 비율에 대한 사후 매개 변수 추정을 얻고 싶습니다. 두 박하에서 나오는 동전 (세타 1, 2, 3, 4)과 각 박하에서 나온 동전의 평균 편견에 대한 견적 (민트 바이어스 : 오메가). 각 박하의 편차 인 kappa의 변동성을 상수로 유지했습니다. 문제는 두 번째 박하에서 사후 평가를 얻을 수 없다는 것입니다. 이전 샘플링을 샘플링하는 것처럼 보입니다. 누구든지 두 번째 민트에 대한 사후 평가를 생성하기 위해 모델 문자열 텍스트 (아래 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") 

Answer 1

문제 아래의 분석을위한

전체 스크립트는 theta에 이전을 설정할 때 인덱스 동전의 민트을 시도했던 방법이었다. 이 경우에만 theta이 4 개 있으며 nFlips이 아닙니다. 중첩 된 인덱싱 mints[coins]은 mints 데이터 벡터에 액세스하는 것이지 각 동전이 어떤 동전인지를 나타내는 벡터가 아닙니다. 아래 수정 버전을 만들었습니다. 각 동전이 어느 동전에 속하는지를 나타내는 벡터의 명시적인 구성에 주목하십시오. 또한 모델 사양에서 각 for-loop 인덱스는 데이터 이름과 다른 고유 한 인덱스 이름을가집니다.

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