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

Question

我正在使用psclR中的软件包并尝试使其生成可测试/可重现的结果.我已经看到了在底层的C代码,它看起来好像GetRNGstate()和PutRNGstate()被称为在正确的地方,但它似乎是不可能从MCMC模型重复输出.

我已经simulationResult从SoDA包中打包了函数,因此我可以验证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.

Answer 1

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

在ideal.c所有抽样中出现的部分动工迭代,即其中的功能updatex,updatey和其他人使用.然而,问题在于这些函数的一个参数 - 矩阵ok(具有讽刺意味,对吧？).使用它updatex,并updateb与只在那里的位置ok == 1是重要的(在 crosscheck,crosscheckx).

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

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

ok = imatrix(n,m);

它分配一个整数矩阵(参见util.c imatrix).问题是,然后ok包含各种数字,即不仅是零,有时是零.看起来它们与R的RNG状态无关,这解释了@SchaunW注意到的行为:all.equal(run1@result$xbar, run2@result$xbar)返回TRUEif !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.