R概率模拟不会终止？答案

【问题标题】：R probability simulation that won't terminate?R概率模拟不会终止？
【发布时间】：2013-11-10 11:18:52
【问题描述】：

我正在教一门统计课，让学生通过使用 R 进行模拟来探索概率和统计方面的问题。最近，对于掷 5 个骰子时恰好得到两个 6 的概率有些困惑。答案是choose(5,2)*5^3/6^5，但也有同学认为“顺序不重要”；即答案应该是choose(5,2)*choose(25,3)/choose(30,5)。我认为让他们模拟滚动 5 个骰子数千次，跟踪每个实验的经验概率，然后重复实验多次会很有趣。问题是上面的两个数字足够接近，以至于很难通过模拟以统计上显着的方式梳理出差异（当然我可能做错了）。我试着掷 5 个骰子 100000 次，然后重复实验 10000 次。这需要一个小时左右才能在我的 i7 linux 机器上运行，并且仍然有 25% 的机会正确答案是选择（5,2）*选择（25,3）/选择（30,5）。所以我将每个实验的掷骰数增加到 10^6。现在代码已经运行了 2 多天，并且没有完成的迹象。我对此感到困惑，因为我只是将操作数量增加了一个数量级，这意味着运行时间应该接近 10 小时。

第二个问题：有没有更好的方法来做到这一点？请参阅下面发布的代码：

probdist = rep(0,10000)

for (j in 1:length(probdist))
{
   outcome = rep(0,1000000)
   for (k in 1:1000000)
   {
      rolls = sample(1:6, 5, replace=T)
      if (length(rolls[rolls == 6]) == 2) outcome[k] = 1
   }

   probdist[j] = sum(outcome)/length(outcome)
}

【问题讨论】：

标签： r probability simulation

【解决方案1】：

一个好的经验法则是永远不要在R 中编写for 循环。这是一个替代解决方案：

doSample <- function()
{
   sum(sample(1:6,size=5,replace=TRUE)==6)==2
}

> system.time(samples <- replicate(n=10000,expr=doSample()))
user  system elapsed 
0.06    0.00    0.06 
> mean(samples)
[1] 0.1588
> choose(5,2)*5^3/6^5
[1] 0.160751

10,000 美元的样本似乎不太准确。 100,000 美元更好：

> system.time(samples <- replicate(n=100000,expr=doSample()))
user  system elapsed 
0.61    0.02    0.61 
> mean(samples)
[1] 0.16135

【讨论】：

谢谢，这非常有帮助。显然需要超过 100000 次掷骰子：> choose(5,2)*5^3/6^5 = 0.160751, choose(5,2)*choose(25,3)/choose(30,5) = 0.1613967，但这比我做的要快数千倍，而且代码更简单易懂。
知道为什么 for 循环在 R 中这么慢，而 replicate() 这么快吗？
此答案中的代码进行了 1e5 次迭代。您问题中的代码执行 1e4*1e6=1e10 次迭代。

【解决方案2】：

我最初授予 M. Berk 一个正确答案检查，因为他/她建议使用 R replicate() 函数。进一步的调查迫使我撤销我之前的认可。事实证明，replicate() 只是 sapply() 的一个包装器，它实际上并没有为 for 循环提供任何性能优势（这似乎是一个常见的误解）。无论如何，我准备了 3 个版本的模拟，2 个使用 for 循环，一个使用复制，如建议的那样，并一个接一个地运行它们，每次从一个新的 R 会话开始，以比较执行时间：

# dice26dist1.r: For () loop version with unnecessary array allocation
probdist = rep(0,100)

for (j in 1:length(probdist))
{
  outcome = rep(0,1000000)
  for (k in 1:1000000)
  {
    rolls = sample(1:6, 5, replace=T)
    if (length(rolls[rolls == 6]) == 2) outcome[k] = 1
  }
  probdist[j] = sum(outcome)/length(outcome)
}

system.time(source('dice26dist1.r'))
用户系统已过
596.365 0.240 598.614

# dice26dist2.r: For () loop version
probdist = rep(0,100)

for (j in 1:length(probdist))
{
  outcomes = 0
  for (k in 1:1000000)
  {
    rolls = sample(1:6, 5, replace=T)
    if (length(rolls[rolls == 6]) == 2) outcomes = outcomes + 1
  }
  probdist[j] = outcomes/1000000
}

system.time(source('dice26dist2.r'))
用户系统已过
506.331 0.076 508.104

# dice26dist3.r:  replicate() version
doSample <- function()
{
   sum(sample(1:6,size=5,replace=TRUE)==6)==2
}

probdist = rep(0,100)

for (j in 1:length(probdist))
{
  samples = replicate(n=1000000,expr=doSample())
  probdist[j] = mean(samples)
}

system.time(source('dice26dist3.r'))
用户系统已过
804.042 0.472 807.250

从这里可以看出，从任何 system.time 指标来看，replicate() 版本都比任何一个 for 循环版本都慢。我原本以为我的问题主要是通过分配百万字符的结果[]数组导致缓存未命中，但是比较 dice26dist1.r 和 dice26dist2.r 的时间表明这对性能只有名义上的影响（尽管对系统的影响时间相当可观：>300% 差异。

有人可能会争辩说，我在所有三个模拟中仍然使用 for 循环，但据我所知，在模拟随机过程时这是完全不可避免的；我每次都必须模拟实际经历随机过程（在这种情况下，滚动 5 个骰子）。我很想知道任何可以让我避免使用 for 循环的技术（当然，以提高性能的方式）。我知道这个问题非常适合并行化，但我说的是使用单个 R 会话——有没有办法让它更快？

【讨论】：

【解决方案3】：

向量化几乎总是优于任何 for 循环。在这种情况下，通过首先生成所有掷骰子，然后检查每组 5 个等于 6 的骰子数量，您应该会看到显着的加速。

set.seed(5)
N <- 1e6
foo <- matrix(sample(1:6, 5*N, replace=TRUE), ncol=5)
p <- mean(rowSums(foo==6)==2)
se <- sqrt(p*(1-p)/N)
p
## [1] 0.160382

这是一个 95% 的置信区间：

p + se*qnorm(0.975)*c(-1,1)
## [1] 0.1596628 0.1611012

我们可以看到正确答案（ans1）在区间内，而错误答案（ans2）不在，或者我们可以进行显着性检验；测试正确答案时的 p 值为 0.31，但测试错误答案时的 p 值为 0.0057。

(ans1 <- choose(5,2)*5^3/6^5)
## [1] 0.160751
pnorm(abs((ans1-p)/se), lower=FALSE)*2
## [1] 0.3145898

ans2 <- choose(5,2)*choose(25,3)/choose(30,5)
## [1] 0.1613967
pnorm(abs((ans2-p)/se), lower=FALSE)*2
## [1] 0.005689008

请注意，我一次生成所有掷骰子；如果内存是个问题，您可以将其拆分并组合，就像您在原始帖子中所做的那样。这可能是导致您意外加速的原因；如果有必要使用交换内存，这将大大减慢它。如果是这样，最好增加运行循环的次数，而不是循环内的滚动次数。

【讨论】：