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

我正在使用R中的pscl包并试图让它产生可测试/可重现的结果.我已经看了底层的C代码,看起来似乎在正确的位置调用GetRNGstate()和PutRNGstate(),但似乎无法重复MCMC模型的输出.

我已经在SoDA包中的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)

我可以增加迭代次数,但如果RNG状态得到传播则这并不重要.我错过了一些非常简单的事吗?谢谢.

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

EDIT2

我还应该在OS X 10.8.4上使用pscl 1.04.4添加我的R 3.0.1.

正如OP和@SchaunW所怀疑的那样,问题在于C代码. “挖掘一点”揭示了一个相当微妙的问题(参见 source代码,但不是最新版本):

ideal.c中的所有采样都出现在开始迭代的部分中,即使用函数updatex,updatey和其他函数.然而,问题在于这些函数的一个参数 – 矩阵确定(具有讽刺意味,对吧?).它由updatex和updateb使用,只有ok == 1的位置很重要(in crosscheck,crosscheckx).

在此之前,某些ok值被指定为1(y,ok,n,m).

但是,在最开始时,ok的初始值表示为

ok = imatrix(n,m);

它分配一个整数矩阵(参见ima.c的util.c).问题是那么ok包含各种数字,即不仅是零,有时是零.看起来它们与R的RNG状态无关,这解释了@SchaunW注意到的行为:all.equal(run1 @ result $xbar,run2 @ result $xbar)如果!any(ok == 1)则返回TRUE反之亦然.另外,不同数量的解释了不同的lastState.

我不是C的专家,我不确定代码中是否存在逻辑错误,或者是否应该纠正imatrix函数,但是直接修复可能是在初始化后立即用零填充:

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时理想.

相关文章
相关标签/搜索