R2Jags :: autojags（）＆＃x27;函数不会收敛

Question

我正在运行以下代码：

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（）函数的模型），所以它应该不很重要。到底是怎么回事？

Answer 1

这就是autojags（）的代码。我们将仔细查看其工作原理：

1. 使用您在功能调用中指示的稀疏和迭代更新对象。

   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, ...)

2. 检查是否有任何R-HAT值大于函数调用中指定的R-hat值。

    check <- any(object$BUGSoutput$summary[, "Rhat"] > Rhat)

3. 如果R-hat值都小于函数调用中的值，则停止，否则继续进行。

    if (check) {

4. 初始化计数器，虽然该计数器小于n.update（默认为2），而至少一个R-hat值大于函数调用中指示的一个值，请执行操作在牙套中。

        count <- 1
        while (check & (count < n.update)) {

5. 再次更新先前估计的模型，并使用相同的迭代，变薄等...

            object <- update(object, n.iter = n.iter, n.thin = n.thin, 
                refresh = refresh, progress.bar = progress.bar, 
                ...)

6. 递增count计数器并再次检查收敛。

            count <- count + 1
            check <- any(object$BUGSoutput$summary[, "Rhat"] > 
                Rhat)

7. 继续循环直到满足条件，然后

        }
    }
    return(object)

将调用中的结果返回autojags（），您正在使用默认n.iter（1000）和默认代码> n.update （2）。这意味着该函数最多将停止3000次迭代（第一次运行1000个，而对于2个更新中的每个更新中的每一个）都将停止。贝叶斯模型仅被保证会融合，因为无限距离的绘制数量。停止规则存在，因此您的计算机不会尝试在无限的时间内运行。还值得注意的是，检查不是按迭代进行迭代，而是n.Iter -by- n.iter块。