我对贝叶斯统计比较陌生,并且正在尝试使用 R2winBUGS 对一些树木放养密度数据应用非线性分层模型。我希望有人可以帮助我找到 R2winBUGS 给我以下错误的原因:
gen.inits 无法执行(灰显)
即使我收到此错误,代码仍会产生输出。对于收敛的参数,模型产生的输出看起来很合理,但是当我使用两条链时,有两个参数(下面代码中的 mean_N0 sigma_N0)不能很好地混合(不收敛)。链似乎从初始值开始(即 mean_N0 从 4800.5 和 4799.5 开始,而 sigma_N0 从 800.5 和 799.5 开始),但与这些值相差不远。两个参数的平均值与初始设定值相差约 0.5。我不确定上述错误是否导致了这个收敛问题。
我已经用尽了对这个问题的调查,现在希望有人能够在下面的 winBUGS 或 R 代码中看到导致我的问题的原因。如果您能提供帮助,我将不胜感激。亲切的问候多姆
WINBUGS 代码
model {
## loop over data for likelihood
for(i in 1:Ntotal){
N[i] <- log(N0[P_ID_Bug[i]] - 25)-(Age[i]/(Beta_0 + Beta_1*Age[i]))
Y[i] ~ dnorm(N[i],tauY)
}
tauY ~ dgamma(1.0E-3, 1.0E-3)
Beta_0 ~ dnorm(9,0.25)
Beta_1 ~ dnorm(0.16,400)
## hierarchical model for each Plots intercept & slope
for (p in 1:P_ID_Length) {
N0[p] ~ dgamma(r_N0, lambda_N0)
}
mean_N0 ~ dnorm(5000,1.0E-6)
sigma_N0 ~ dnorm(5000,0.25E-6)
lambda_N0 <- mean_N0/(sigma_N0*sigma_N0)
r_N0 <- mean_N0 * lambda_N0
}
代码
data <- list(P_ID_Length = length(P_ID),
P_ID_Bug = P_ID_Bug,
Age=Grouped_SDen$Age,
Y =Grouped_SDen$LogSDen_Ha,
Ntotal=nrow(Grouped_SDen))
inits1 <- list(N0= rep(coef(NLS_SDen_Log_1)[[1]], P_ID_Length),
Beta_0 = coef(NLS_SDen_Log_1)[[2]],
Beta_1 = coef(NLS_SDen_Log_1)[[3]],
tauY = 20,
mean_N0 = 4800,
sigma_N0 = 800
)
inits2 <- list(N0= rep(coef(NLS_SDen_Log_1)[[1]], P_ID_Length),
Beta_0 = coef(NLS_SDen_Log_1)[[2]],
Beta_1 = coef(NLS_SDen_Log_1)[[3]],
tauY = 20,
mean_N0 = 4800,
sigma_N0 = 800
)
inits <- list(inits1, inits2)
parameters <- c("N0", "Beta_0", "Beta_1", "tauY", "mean_N0", "sigma_N0")
sims <- bugs(model.file= "C:/WS/Post-Doc/TINNR/PAPER_4_WinBugs/Stocking_Growth.bug",
data = data,
parameters = parameters,
inits = inits,
n.chains = 2,
n.iter = 1000,
n.burnin = 500,
n.thin = 2,
debug=TRUE,
bugs.directory = "C:/Program Files/WinBUGS14/")
我也尝试过使用以下格式进行初始化:
inits <- function(){
list(N0= rep(coef(NLS_SDen_Log_1)[[1]], P_ID_Length),
Beta_0 = coef(NLS_SDen_Log_1)[[2]],
tauY = 1,
mean_N0 = 4800,
sigma_N0 = 800
)
}
请注意,第一个和第二个初始化是相同的,下面我将向您展示第二个列表的外观。再次感谢您抽出时间提供帮助。
[[2]]
[[2]]$N0
[1] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[13] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[25] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[37] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[49] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[61] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[73] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[85] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[97] 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531 4920.531
[[2]]$Beta_0
[1] 8.710965
[[2]]$Beta_1
[1] 0.1623536
[[2]]$tauY
[1] 20
[[2]]$mean_N0
[1] 4800
[[2]]$sigma_N0
[1] 800