'R2jags::autojags()' 函数不收敛
The 'R2jags::autojags()' Function Doesn't Converge
我正在运行宁以下代码：
Data_Frame <- structure(list(Predictor_Variable = c(20.5273283761926, 11.4039835403673, 5.31423215288669, 12.5192134582903, 15.8487716107629, 7.16369490255602, 3.91227747895755, 8.90797802712768, 14.1797202872112, 4.82832778943703, 14.3634092528373, 18.3354590029921, 1.56595673179254, 15.27424925589, 2.01471757609397, 17.8340036130976, 24.4356162613258, 24.7056286665611, 8.86376699781977, 3.95776731311344, 1.70528281596489, 24.1929906886071, 7.84694050671533, 2.11047385819256, 24.8220994370058, 14.8339046107139, 19.9926545028575, 19.0412578289397, 4.78973534191027, 18.7796145037282, 19.3868586677127, 14.6739382355008, 20.0653701729607, 19.6233123773709, 7.19542927108705, 16.7649424169213, 9.41006114007905, 19.2996966594364, 21.0871973482426, 13.0715743871406, 17.3938042367809, 18.3959823509213, 9.12836091592908, 2.27788038901053, 10.3573848318774, 4.74680115585215, 9.6620995842386, 16.9264123018365, 11.6498990275431, 20.4607627529185, 24.9407043796964, 21.8109163164627, 11.807474057423, 0.530748104210943, 9.57337225554511, 16.5205421217252, 2.07126556197181, 14.2894889519084, 5.47687077778392, 17.44882510975, 14.0162174822763, 11.7680913361255, 1.04056366835721, 1.95715780719183, 21.3338757341262, 19.7388443630189, 18.9714497013483, 19.5131896412931, 14.1728786868043, 10.4211826983374, 0.625250476878136, 17.7159008861054, 17.5371950492263, 3.91660461900756, 20.3449436288793, 1.5704334480688, 13.397057261318, 20.7767494546715, 7.27919469354674, 23.8504881446715, 19.7164187382441, 5.95191353349946, 20.7321806927212, 10.5640658119228, 16.616114001954, 15.9614237083588, 3.23146634036675, 16.6699483001139, 14.8473161039874, 18.7451402307488, 3.54583212174475, 5.97735904157162, 8.23308551334776, 3.3043225936126, 9.93713270290755, 5.86233736248687, 2.55345033365302, 11.350765300449, 7.88409785600379, 16.5159953408875), Response_Variable = c(1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L)), class = "data.frame", row.names = c(NA, -100L))
sink("Generalised Logistic Function Model.txt")
cat("model {
# Priors
Intercept ~ dnorm(0, 0.001)
Slope ~ dnorm(0, 0.001)
Exponent ~ dlnorm(0, 1)
Sigma ~ dlnorm(0, 1)
Tau <- (1 / (Sigma ^ 2))
# Likelihood and Model Fit
for (i in 1:Number_of_Observations) {
Response[i] ~ dnorm(Mean[i], Tau)
Mean[i] <- ((1 + exp(-(Intercept + (Slope * Predictor[i])))) ^ (-Exponent))
Actual_Squared_Residual[i] <- ((Response[i] - Mean[i]) ^ 2)
Simulated_Response[i] ~ dbern(Mean[i])
Simulated_Squared_Residual[i] <- ((Simulated_Response[i] - Mean[i]) ^ 2)
}
Bayesian_p_Value <- step((sum((Simulated_Squared_Residual[]) ^ 2)) / (sum((Actual_Squared_Residual[]) ^ 2)) - 1)
}", fill = T)
sink()
Parameters <- c("Intercept", "Slope", "Exponent", "Bayesian_p_Value")
Data <- list(Predictor = Data_Frame$Predictor_Variable, Response = Data_Frame$Response_Variable, Number_of_Observations = 100)
Initial_Values <- function () {
list()
}
(Output_1 <- R2jags::jags(Data, Initial_Values, Parameters, "Generalised Logistic Function Model.txt", n.chains = 3, n.thin = 1, n.iter = 1000, n.burnin = 100, working.directory = getwd()))
当我运行这个贝叶斯模型时，我得到以下输出：
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 100 discarded)
n.sims = 2700 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.436 0.496 0.000 0.000 0.000 1.000 1.000 1.003 870
Exponent 1.122 1.002 0.149 0.431 0.811 1.412 3.900 1.189 16
Intercept -3.550 3.022 -11.786 -4.783 -2.785 -1.452 0.148 1.299 14
Slope 0.259 0.124 0.117 0.173 0.221 0.300 0.597 1.197 18
deviance 112.166 2.789 108.793 110.218 111.477 113.362 119.313 1.008 270
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 3.9 and DIC = 116.0
DIC is an estimate of expected predictive error (lower deviance is better).
然后当我尝试使用 R2jags::autojags() 函数时，我得到以下输出：
(Output_2 <- R2jags::autojags(Generalised_Logistic_Function_Model_Output, Rhat = 1.05))
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 0 discarded)
n.sims = 3000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.377 0.485 0.000 0.000 0.000 1.000 1.000 1.020 110
Exponent 0.889 0.831 0.037 0.339 0.670 1.186 3.195 1.757 6
Intercept -8.189 12.251 -47.264 -5.649 -3.372 -1.833 -0.135 1.831 6
Slope 0.473 0.560 0.126 0.191 0.241 0.350 2.130 1.871 6
deviance 112.520 2.919 108.919 110.427 111.930 113.896 119.283 1.103 28
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 4.0 and DIC = 116.5
DIC is an estimate of expected predictive error (lower deviance is better).
我的印象是这个函数会继续 运行 直到它收敛（直到 Rhat 达到 1.1）。事实并非如此 - 第二个模型实际上具有更高的 Rhat 值。我试过在 autojags() 函数中使用 Rhat 参数，但这没有帮助。我可以看到第二个模型没有丢弃任何初始迭代，但由于它是从一个已经不错的模型开始的（我使用 jags() 函数得到的模型），所以这无关紧要。怎么回事？
这就是 autojags() 的代码。我们将通过它来了解它是如何工作的：
- 使用您在函数调用中指定的细化和迭代更新对象。
n.burnin = object$n.iter
n.thin.auto <- max(1, floor((n.iter - n.burnin)/1000))
n.thin <- ifelse(n.thin > n.thin.auto, n.thin, n.thin.auto)
if (any(!class(object) %in% c("rjags", "rjags.parallel")))
stop("model must be a rjags object")
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar, ...)
- 检查是否有任何 R-hat 值大于函数调用中指定的 R-hat 值。
check <- any(object$BUGSoutput$summary[, "Rhat"] > Rhat)
- 如果R-hat值都小于函数调用中的值，则停止，否则继续。
if (check) {
- 初始化一个计数器，当该计数器小于
n.update（默认为 2）并且至少有一个 R-hat 值大于函数调用中指示的值时，做大括号里的事。
count <- 1
while (check & (count < n.update)) {
- 再次更新先前估计的模型，使用相同的迭代、细化等...
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar,
...)
- 增加
count 计数器并再次检查收敛。
count <- count + 1
check <- any(object$BUGSoutput$summary[, "Rhat"] >
Rhat)
- 继续循环直到满足条件return结果
}
}
return(object)
在您对 autojags() 的调用中，您使用的是默认值 n.iter (1000) 和默认值 n.update (2)。这意味着该函数将在最多 3000 次迭代时停止（第一个 运行 到 1000 次，两次更新中的每一次更新 1000 次），无论是否收敛。贝叶斯模型只能保证在从后验抽取的次数接近无穷大时收敛。存在停止规则，因此您的计算机不会无限期地尝试 运行。还值得注意的是，检查并未完成 iteration-by-iteration，而是 n.iter-by-n.iter 个块。
我正在运行宁以下代码：
Data_Frame <- structure(list(Predictor_Variable = c(20.5273283761926, 11.4039835403673, 5.31423215288669, 12.5192134582903, 15.8487716107629, 7.16369490255602, 3.91227747895755, 8.90797802712768, 14.1797202872112, 4.82832778943703, 14.3634092528373, 18.3354590029921, 1.56595673179254, 15.27424925589, 2.01471757609397, 17.8340036130976, 24.4356162613258, 24.7056286665611, 8.86376699781977, 3.95776731311344, 1.70528281596489, 24.1929906886071, 7.84694050671533, 2.11047385819256, 24.8220994370058, 14.8339046107139, 19.9926545028575, 19.0412578289397, 4.78973534191027, 18.7796145037282, 19.3868586677127, 14.6739382355008, 20.0653701729607, 19.6233123773709, 7.19542927108705, 16.7649424169213, 9.41006114007905, 19.2996966594364, 21.0871973482426, 13.0715743871406, 17.3938042367809, 18.3959823509213, 9.12836091592908, 2.27788038901053, 10.3573848318774, 4.74680115585215, 9.6620995842386, 16.9264123018365, 11.6498990275431, 20.4607627529185, 24.9407043796964, 21.8109163164627, 11.807474057423, 0.530748104210943, 9.57337225554511, 16.5205421217252, 2.07126556197181, 14.2894889519084, 5.47687077778392, 17.44882510975, 14.0162174822763, 11.7680913361255, 1.04056366835721, 1.95715780719183, 21.3338757341262, 19.7388443630189, 18.9714497013483, 19.5131896412931, 14.1728786868043, 10.4211826983374, 0.625250476878136, 17.7159008861054, 17.5371950492263, 3.91660461900756, 20.3449436288793, 1.5704334480688, 13.397057261318, 20.7767494546715, 7.27919469354674, 23.8504881446715, 19.7164187382441, 5.95191353349946, 20.7321806927212, 10.5640658119228, 16.616114001954, 15.9614237083588, 3.23146634036675, 16.6699483001139, 14.8473161039874, 18.7451402307488, 3.54583212174475, 5.97735904157162, 8.23308551334776, 3.3043225936126, 9.93713270290755, 5.86233736248687, 2.55345033365302, 11.350765300449, 7.88409785600379, 16.5159953408875), Response_Variable = c(1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L)), class = "data.frame", row.names = c(NA, -100L))
sink("Generalised Logistic Function Model.txt")
cat("model {
# Priors
Intercept ~ dnorm(0, 0.001)
Slope ~ dnorm(0, 0.001)
Exponent ~ dlnorm(0, 1)
Sigma ~ dlnorm(0, 1)
Tau <- (1 / (Sigma ^ 2))
# Likelihood and Model Fit
for (i in 1:Number_of_Observations) {
Response[i] ~ dnorm(Mean[i], Tau)
Mean[i] <- ((1 + exp(-(Intercept + (Slope * Predictor[i])))) ^ (-Exponent))
Actual_Squared_Residual[i] <- ((Response[i] - Mean[i]) ^ 2)
Simulated_Response[i] ~ dbern(Mean[i])
Simulated_Squared_Residual[i] <- ((Simulated_Response[i] - Mean[i]) ^ 2)
}
Bayesian_p_Value <- step((sum((Simulated_Squared_Residual[]) ^ 2)) / (sum((Actual_Squared_Residual[]) ^ 2)) - 1)
}", fill = T)
sink()
Parameters <- c("Intercept", "Slope", "Exponent", "Bayesian_p_Value")
Data <- list(Predictor = Data_Frame$Predictor_Variable, Response = Data_Frame$Response_Variable, Number_of_Observations = 100)
Initial_Values <- function () {
list()
}
(Output_1 <- R2jags::jags(Data, Initial_Values, Parameters, "Generalised Logistic Function Model.txt", n.chains = 3, n.thin = 1, n.iter = 1000, n.burnin = 100, working.directory = getwd()))
当我运行这个贝叶斯模型时，我得到以下输出：
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 100 discarded)
n.sims = 2700 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.436 0.496 0.000 0.000 0.000 1.000 1.000 1.003 870
Exponent 1.122 1.002 0.149 0.431 0.811 1.412 3.900 1.189 16
Intercept -3.550 3.022 -11.786 -4.783 -2.785 -1.452 0.148 1.299 14
Slope 0.259 0.124 0.117 0.173 0.221 0.300 0.597 1.197 18
deviance 112.166 2.789 108.793 110.218 111.477 113.362 119.313 1.008 270
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 3.9 and DIC = 116.0
DIC is an estimate of expected predictive error (lower deviance is better).
然后当我尝试使用 R2jags::autojags() 函数时，我得到以下输出：
(Output_2 <- R2jags::autojags(Generalised_Logistic_Function_Model_Output, Rhat = 1.05))
Inference for Bugs model at "Generalised Logistic Function Model.txt", fit using jags,
3 chains, each with 1000 iterations (first 0 discarded)
n.sims = 3000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
Bayesian_p_Value 0.377 0.485 0.000 0.000 0.000 1.000 1.000 1.020 110
Exponent 0.889 0.831 0.037 0.339 0.670 1.186 3.195 1.757 6
Intercept -8.189 12.251 -47.264 -5.649 -3.372 -1.833 -0.135 1.831 6
Slope 0.473 0.560 0.126 0.191 0.241 0.350 2.130 1.871 6
deviance 112.520 2.919 108.919 110.427 111.930 113.896 119.283 1.103 28
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 4.0 and DIC = 116.5
DIC is an estimate of expected predictive error (lower deviance is better).
我的印象是这个函数会继续 运行 直到它收敛（直到 Rhat 达到 1.1）。事实并非如此 - 第二个模型实际上具有更高的 Rhat 值。我试过在 autojags() 函数中使用 Rhat 参数，但这没有帮助。我可以看到第二个模型没有丢弃任何初始迭代，但由于它是从一个已经不错的模型开始的（我使用 jags() 函数得到的模型），所以这无关紧要。怎么回事？
这就是 autojags() 的代码。我们将通过它来了解它是如何工作的：
- 使用您在函数调用中指定的细化和迭代更新对象。
n.burnin = object$n.iter
n.thin.auto <- max(1, floor((n.iter - n.burnin)/1000))
n.thin <- ifelse(n.thin > n.thin.auto, n.thin, n.thin.auto)
if (any(!class(object) %in% c("rjags", "rjags.parallel")))
stop("model must be a rjags object")
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar, ...)
- 检查是否有任何 R-hat 值大于函数调用中指定的 R-hat 值。
check <- any(object$BUGSoutput$summary[, "Rhat"] > Rhat)
- 如果R-hat值都小于函数调用中的值，则停止，否则继续。
if (check) {
- 初始化一个计数器，当该计数器小于
n.update（默认为 2）并且至少有一个 R-hat 值大于函数调用中指示的值时，做大括号里的事。
count <- 1
while (check & (count < n.update)) {
- 再次更新先前估计的模型，使用相同的迭代、细化等...
object <- update(object, n.iter = n.iter, n.thin = n.thin,
refresh = refresh, progress.bar = progress.bar,
...)
- 增加
count计数器并再次检查收敛。
count <- count + 1
check <- any(object$BUGSoutput$summary[, "Rhat"] >
Rhat)
- 继续循环直到满足条件return结果
}
}
return(object)
在您对 autojags() 的调用中，您使用的是默认值 n.iter (1000) 和默认值 n.update (2)。这意味着该函数将在最多 3000 次迭代时停止（第一个 运行 到 1000 次，两次更新中的每一次更新 1000 次），无论是否收敛。贝叶斯模型只能保证在从后验抽取的次数接近无穷大时收敛。存在停止规则，因此您的计算机不会无限期地尝试 运行。还值得注意的是，检查并未完成 iteration-by-iteration，而是 n.iter-by-n.iter 个块。