估计 n 次掷出 m 个六面骰子的机会

Question

我想生成一个 Monte Carlo 模拟来估计从 n 掷 m 个公平的六面骰子中掷出至少一个的概率。

我的代码：

m<-5000
n<-3
x<-replicate(m, sample(1:6,n,TRUE)==1)
p<-sum(x)/m

p 是估计的概率。在这里我得到值 0.4822.

我的问题：

1) 有没有不用sum的其他方法呢？

2) 我怀疑代码是错误的，因为概率可能太高了。

Answer 1

虽然所陈述的问题有点不清楚，但代码表明您想要估计在 n 个独立骰子中获得至少一个结果为“1”的机会，并且您的目标是通过模拟实验 m 次。

从内到外的程序模拟。从单次迭代开始。你开始得很好，但要非常清楚，让我们使用高度暗示的语法重做它。试试这个：

1 %in% sample(1:6,n,TRUE)

这里使用sample来实现n个独立公平骰子的结果，并检查结果1是否出现在其中任何一个

一旦您对模拟您的实验感到满意（运行多次），那么 replicate 确实会执行模拟：

x <- replicate(m, 1 %in% sample(1:6,n,TRUE))

这会产生 m 个结果。在 1 出现的所有迭代中，每个都将为 TRUE（解释为等于 1），否则将为 FALSE（解释为 0）。因此，1出现的平均次数可以得到

mean(x)

这个 经验频率 是对 理论概率 的良好估计。

作为检查，请注意 1 不会不会出现在单个骰子上的概率为 1-1/6 = 5/6，因此——因为 n 个骰子是独立的—— 不会出现在任何一个骰子上的概率为 (5/6)^n。因此 1 将出现的机会必须是 1 - (5/6)^n。让我们输出这两个值：模拟平均值和理论结果。我们可能还会包括 Z 分数 ，这是衡量均值与理论结果相差多远的指标。通常，-2 和 2 之间的 Z 分数不是任何差异的重要证据。

这是完整的代码。尽管有更快的方法来编写它，但这已经非常快并且尽可能清晰。

m <- 5000     # Number of simulation iterations
n <- 3        # Number of dice per iteration 
set.seed(17)  # For reproducible results
x <- replicate(m, 1 %in% sample(1:6,n,TRUE))

# Compare to a theoretical result.
theory <- 1-(5/6)^n
avg <- mean(x)
Z <- (avg - theory) / sd(x) * sqrt(length(x))
c(Mean=signif(avg, 5), Theoretical=signif(theory, 5), Z.score=signif(Z, 3))

输出为

Mean Theoretical Z.score

0.4132 0.4213 -1.1600

请注意，这两个结果都不接近 n/6，即 1/2 = 0.500。

估计 n 次掷出 m 个六面骰子的机会

Estimate the chance n rolls of m fair six-sided dice

r

dice