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 中的所有采样都出现在开始迭代的部分，即使用函数updatex和updatey其他函数的地方。然而，问题在于这些函数的参数之一——矩阵ok（讽刺，对吧？）。它由updatex并且updateb仅用于重要的位置ok == 1（in crosscheck，crosscheckx）。

在此之前， 的某些值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 回答