Run Jags — извлечь несколько реализаций из объекта mcmc

У меня есть скрипт runjags, который генерирует предсказанную плотность нор для каждой клетки на острове. Я хочу получить несколько рисунков (около 100) от объекта mcmc для каждой ячейки. Мой научный руководитель считает, что я смогу сделать это с помощью пакета 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

Заранее спасибо за любые ответы.


mcmc jags runjags coda
person Bong112    schedule 21.07.2018    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Да, вы можете сделать это, просто извлекая объект MCMC из объекта runjags. Пример модели:

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 итераций, а функция comb.mcmc гарантирует, что они будут равномерно удалены от апостериорной, чтобы минимизировать эффекты автокорреляции внутри цепочек.

Или использовать инструменты в пакете coda, чтобы сделать то же самое:

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

Мэтт

person Matt Denwood    schedule 06.09.2018
comment
Спасибо за это, Мэтт. Мне удалось получить информацию с помощью tidybayes, который показался мне немного более интуитивным, чем coda. - person Bong112; 19.09.2018