ideal()在R包中不能产生可重复的结果

Question

我正在使用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调用的RNG函数之外的一些偶然性来源正在影响这些RNG函数被调用的次数。好奇。

Edit2

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

Answer 1

OP和@SchaunW怀疑问题出在C代码上。"一点点"挖掘揭示了一个相当微妙的问题(请参阅源代码，而不是最新版本):

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

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

然而，一开始ok的初始值用

表示

ok = imatrix(n,m);

分配一个整数矩阵(参见util.c中的imatrix)。问题是ok包含各种数字，即不仅是0，有时还包含1。它们似乎与R的RNG状态无关，这解释了@SchaunW注意到的行为:如果!any(ok == 1),all.equal(run1@result$xbar, run2@result$xbar)返回TRUE，反之亦然。另外，不同的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代码的修复(尽管它可能不适合您的应用程序)。当impute = TRUE用于ideal时，用crossxyi、crossxyj代替crosscheck、crosscheckx(坏的)。

Answer 2

EDIT

我无法复制我最初发布的结果。当我第一次得到这些结果时，我关闭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的调用产生效果，您需要将该参数更改为"random"。这是你已经尝试过的:

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设置为"random"的同样的事情:

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设置为"随机"，在包文档中没有明确指出。我摆弄了一会儿才弄明白。

Answer 3

这是一个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)