r - R 包 pscl 中的 Ideal() 不会产生可重复的结果

Question

我正在使用 pscl R 中的包并试图让它产生可测试/可重复的结果。我查看了底层的 C 代码，看起来好像 GetRNGstate()和 PutRNGstate()在正确的地方被调用，但似乎不可能重复来自 MCMC 模型的输出。

我已经把simulationResult中的功能打包了来自 苏打水 包，以便我可以验证 R 端每个模拟 R 的开始状态。

library(pscl)
library(SoDA)
run1 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

run2 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

我们可以验证至少在 R 端的起始状态是相同的:

all.equal(run1@firstState, run2@firstState)

但输出不同:

all.equal(run1@result$xbar, run2@result$xbar)

我可以增加迭代次数，但如果 R​​NG 状态得到传播，这应该无关紧要。我错过了一些非常简单的东西吗？谢谢。

编辑 : 我还应该注意到 all.equal(run1@lastState, run2@lastState) (每次运行的结束状态)应该相同，但它们最终不同。我的猜测是，由 C 调用的 R RNG 函数之外的某些意外源正在影响调用这些 RNG 函数的次数。好奇的。

编辑 2

我还应该在 OS X 10.8.4 上添加我在 R 3.0.1 和 pscl 1.04.4 上。

Answer 1

正如 OP 和@SchaunW 所怀疑的那样，问题在于 C 代码。 “一点点”挖掘揭示了一个非常微妙的问题(请参阅 source 代码，但不是最新版本):

Ideal.c 中的所有采样都出现在开始迭代的部分，即函数 updatex 的位置。 , updatey和其他人使用。然而，问题在于这些函数的参数之一——矩阵 ok (讽刺，对吧？)。它由 updatex 使用和 updateb和 只有 ok == 1所在的位置很重要( in crosscheck ， crosscheckx )。

在此之前，ok 的一些值在 check(y,ok,n,m) 中被分配为 1 .

然而，一开始 ok 的初始值表示为

ok = imatrix(n,m);

它分配一个整数矩阵(请参阅 util.c for imatrix )。问题是，然后 ok包含各种数字，即不仅是零，有时是一。似乎它们与 R 的 RNG 状态无关，这解释了@SchaunW 指出的行为: all.equal(run1@result$xbar, run2@result$xbar)返回 TRUE如果 !any(ok == 1)反之亦然。此外，不同数量的解释不同 lastState .

我不是 C 专家，我不确定代码中是否存在逻辑错误或者 imatrix应该更正函数，但一个简单的解决方法是填充 ok初始化后立即为零:

ok = imatrix(n,m);
for(a=0; a<n; a++) {
    for(aa=0; aa<m; aa++) {
      ok[a][aa] = 0;
    }
}

最后，还有一个不包括修改 C 代码的修复程序(尽管它可能不适合您的应用程序)。函数 crossxyi , crossxyj用于代替 crosscheck , crosscheckx (坏的)当 impute = TRUE为 ideal .