2

我正在教一门统计课，让学生通过使用 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 小时。现在代码已经运行了 2 多天，并且没有完成的迹象。我对此感到困惑，因为我只将操作数量增加了一个数量级，这意味着运行时间应该接近 10 小时。现在代码已经运行了 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)
}
4

3 回答 3

3

一个好的经验法则是永远不要for在R. 这是一个替代解决方案：

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
于 2013-10-30T16:24:27.670 回答
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 会话——有没有办法让它更快？

于 2013-10-30T16:43:28.747 回答
2

向量化几乎总是优于任何 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，但错误答案为 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

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

于 2013-12-09T18:46:58.130 回答