我有一个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
}"
前五行数据集:
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
提前感谢您的回复。
是的,您可以通过从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)
马特