0

我有一个 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

提前感谢您的任何回复。

4

1 回答 1

0

是的，您可以通过简单地从 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)

马特

于 2018-09-06T13:12:47.373 回答