Run Jags -从mcmc对象中提取多个实现

Question

我有一个runjags脚本，可以为岛上的每个单元生成预测的洞穴密度。我希望从每个单元格的mcmc对象中获得多个绘制(大约100个)。我的论文导师认为我应该能够使用coda包来做到这一点，但我只能提取每个单元格的平均值，而不是多个实现。

用于运行模型并提取平均值的代码：

runjags.options(force.summary=TRUE)
print(runjags.options())
S2VS1_best_fit_result <- run.jags(model=S2VS1_best_fit_model, burnin=100000, sample=1000, n.chains=3, modules="glm", thin = 100)
S2_result <- as.mcmc(S2VS1_best_fit_result, vars = "S2")
S2_result_list <- as.mcmc.list(S2VS1_best_fit_result, vars = "S2")
S1_summary <- summary(S2_result_list)
S1_stats <- S2_summary$statistics

谁能告诉我如何为每个单元格获取多个值？

模型：

S2VS1_best_fit_model <- "model{
for(i in 1:K) { # Cells loop
S2[i]~dpois(lambda1[i])
lambda1[i]<- exp(a0+a1*normalise_DEM_aspect[i]+a2*normalise_DEM_elevation[i]+a3*normalise_DEM_slope[i]+
a4*normalise_DEM_elevation[i]*normalise_DEM_slope[i]+
a5*normalise_sentinel5[i]+a6*normalise_sentinel10[i]+
a8*S1[i]+
a9*Tussac[i])

muLogit_tussac[i]<-b0+b1*normalise_sentinel1[i]+b2*normalise_sentinel7[i]+b3*normalise_sentinel8[i]+
               b4*normalise_sentinel9[i]+b5*normalise_DEM_slope[i]

Logit_tussac[i]~dnorm(muLogit_tussac[i], tau) # tau = precision (1/variance or 1/sd^2) - see Lecture 5, Slide 17
Tussac[i]<-exp(Logit_tussac[i])/(1+exp(Logit_tussac[i]))

S1[i]~dpois(lambda2[i])
lambda2[i]<-exp(c0)

}

# Priors

a0~dnorm(0, 10)
a1~dnorm(0, 10)
a2~dnorm(0, 10)
a3~dnorm(0, 10)
a4~dnorm(0, 10)
a5~dnorm(0, 10)
a6~dnorm(0, 10)
a7~dnorm(0, 10)
a8~dnorm(0, 10)
a9~dnorm(0, 10)

b0~dnorm(0, 10)
b1~dnorm(0, 10)
b2~dnorm(0, 10)
b3~dnorm(0, 10)
b4~dnorm(0, 10)
b5~dnorm(0, 10)

c0~dnorm(0, 10)

tau~dgamma(0.001, 0.001)

#data# S1, S2, K
#data# normalise_sentinel1, normalise_sentinel5, normalise_sentinel7
#data# normalise_sentinel9, normalise_sentinel8, normalise_sentinel10
#data# normalise_DEM_aspect, normalise_DEM_elevation, normalise_DEM_slope
#inits# a0, a1, a2, a3, a4, a5
#inits# b0, b1, b2, b3, b4, b5
#inits# c0
#monitor# a0, a1, a2, a3, a4, a5, b0
#monitor# b0, b1, b2, b3, b4, b5
#monitor# c0
#monitor# ped, dic
#monitor# S1, S2
}"

前5行数据集：

S1 S2 Logit_tussac moisture DEM_slope DEM_aspect DEM_elevation sentinel1 sentinel2 sentinel3 sentinel4 sentinel5 sentinel6 sentinel7 sentinel8 sentinel9 sentinel10
NA NA        NA        NA   14.917334   256.1612      12.24432    0.0513    0.0588    0.0541    0.1145    0.1676    0.1988    0.1977    0.1658    0.1566     0.0770
0  0  -9.210240         1   23.803741   225.1231      16.88028    0.1058    0.1370    0.2139    0.2387    0.2654    0.2933    0.3235    0.2928    0.3093     0.1601
NA NA        NA        NA   20.789165   306.0945      18.52480    0.0287    0.0279    0.0271    0.0276    0.0290    0.0321    0.0346    0.0452    0.0475     0.0219
NA NA -9.210240         1    6.689442   287.9641      36.08975    0.0462    0.0679    0.1274    0.1535    0.1797    0.2201    0.2982    0.2545    0.4170     0.2252
0  0  -9.210240         1   25.476444   203.0659      23.59964    0.0758    0.1041    0.1326    0.1571    0.2143    0.2486    0.2939    0.2536    0.3336     0.1937
1  0  -1.385919         3    1.672511   270.0000      39.55215    0.0466    0.0716    0.1227    0.1482    0.2215    0.2715    0.3334    0.2903    0.3577     0.1957

提前感谢您的回复。

Answer 1

可以，您只需从runjags对象中提取MCMC对象即可完成此操作。示例模型：

X <- 1:100
Y <- rnorm(length(X), 2*X + 10, 1)

model <- "model { 
for(i in 1 : N){ 
    Y[i] ~ dnorm(true.y[i], precision);
    true.y[i] <- (m * X[i]) + c
} 
m ~ dunif(-1000,1000)
c ~ dunif(-1000,1000) 
precision ~ dexp(1)
}"

data <- list(X=X, Y=Y, N=length(X))

results <- run.jags(model=model, monitor=c("m", "c", "precision"), 
data=data, n.chains=2)

从中我们可以获得矩阵形式的汇总统计信息：

summary(results)

或者来自后验的给定迭代次数作为MCMC矩阵：

combine.mcmc(results, return.samples=10)

在这种情况下，我们要求10次迭代，combine.mcmc函数确保它们与后验均匀间隔，以便最小化链内自相关的影响。

或者使用coda包中的工具来做同样的事情：

allmcmc <- coda::as.mcmc(results)
window(allmcmc, thin=1000)

哑光