7

我正在使用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)

我可以增加迭代次数,但如果 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 上。

4

3 回答 3

7

正如 OP 和@SchaunW 所怀疑的那样,问题出在 C 代码中。“一点”挖掘揭示了一个非常微妙的问题(参见代码,虽然不是最新版本):

Ideal.c 中的所有采样都出现在开始迭代的部分,即使用函数updatexupdatey其他函数的地方。然而,问题在于这些函数的参数之一——矩阵ok(讽刺,对吧?)。它由updatex并且updateb仅用于重要位置ok == 1in crosscheckcrosscheckx)。

在此之前, 的某些值ok被分配为 1 in 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)返回,反之亦然。此外,不同数量的解释不同。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 代码的修复(虽然它可能不适合您的应用程序)。当for时,函数crossxyi,crossxyj被用来代替crosscheck, crosscheckx(坏的)。impute = TRUEideal

于 2013-06-15T13:44:08.913 回答
3

编辑

我无法重现我最初发布的结果。当我第一次得到这些结果时,我关闭了 R,重新启动它,然后再次运行整个过程以确保,我又得到了相同的结果。下面显示的内容完全是从我的 R 控制台复制的。但是,我只是第三次(以及第四次和第五次)尝试了该代码,但它不起作用。我留下我原来的答案,以防万一我遇到了一些事情并且没有意识到它可能对其他人有用,但下面的建议似乎不起作用(至少不是一致的)。

问题似乎确实出在 C 代码中。当我打开ideal函数并逐行运行时,all.equal这行代码中的每个输入都返回 TRUE:

output <- .C("IDEAL", PACKAGE = .package.Name, as.integer(n), 
      as.integer(m), as.integer(d), as.double(yToC), as.integer(maxiter), 
      as.integer(thin), as.integer(impute), as.integer(mda), 
      as.double(xp), as.double(xpv), as.double(bp), as.double(bpv), 
      as.double(xstart), as.double(bstart), xoutput = as.double(rep(0, 
        n * d * numrec)), boutput = as.double(0), as.integer(burnin), 
      as.integer(usefile), as.integer(store.item), as.character(file), 
      as.integer(verbose))

但是,当我多次运行上述代码output$xoutput时,每次返回的结果都略有不同,即使我set.seed(42)在每次运行前立即调用。

sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] SoDA_1.0-5       pscl_1.04.4      vcd_1.2-13       colorspace_1.2-0 gam_1.06.2       coda_0.16-1      lattice_0.20-10  mvtnorm_0.9-9994
[9] MASS_7.3-22     

loaded via a namespace (and not attached):
[1] tools_2.15.2

原始答案

ideal函数有一个startvals参数。该参数的默认值为“eigen”。为了使您的调用set.seed生效,您需要将该参数更改为“随机”。这是您已经尝试过的:

run1 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "eigen"),
   seed = 42)

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

all.equal(run1@firstState, run2@firstState)
[1] TRUE

all.equal(run1@result$xbar, run2@result$xbar)
[1] "Mean relative difference: 0.01832379"

这与startvals设置为“随机”的情况相同:

run1 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "random"),
   seed = 42)

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

all.equal(run1@firstState, run2@firstState)
[1] TRUE    

all.equal(run1@result$xbar, run2@result$xbar)
[1] TRUE

据我所知,startvals为了获得可复制的结果,需要设置为“随机”并没有在包文档中明确指出。在我弄清楚之前,我不得不玩它一段时间。

于 2013-06-12T20:44:34.067 回答
1

它是一个 MCMC 模型,因此它必须使用随机数生成。要获得可重复的结果,您需要通过为随机数生成器设置“种子”来开始分析。这样每次构建模型时,它都使用相同的“随机”数字(只要每次构建模型时重置种子。使用set.seed()函数并为其提供任意值,例如1234.

我不熟悉这个包,看起来你可能已经在你的函数调用中设置了随机数生成的种子seed=42,但我还是建议明确地设置它set.seed()。然后您的代码变为:

set.seed(1234)
run1 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

set.seed(1234)
run2 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)
于 2013-06-07T17:06:50.837 回答