我已经使用 Rdbinom来生成n=1:32试验的正面频率并现在绘制图表。这将是你所期望的。我在 SO 和 上阅读了您之前的一些帖子math.stackexchange。仍然我不明白你为什么要进行simulate实验而不是从二项式 RV 生成如果你能解释它，那就太好了！我将尝试使用@Andrie 的模拟解决方案来检查我是否可以匹配下面显示的输出。现在，这里有一些你可能感兴趣的东西。
set.seed(42)
numbet <- 32
numtri <- 1e5
prob=5/6
require(plyr)
out <- ldply(1:numbet, function(idx) {
outcome <- dbinom(idx:0, size=idx, prob=prob)
bet <- rep(idx, length(outcome))
N <- round(outcome * numtri)
ymin <- c(0, head(seq_along(N)/length(N), -1))
ymax <- seq_along(N)/length(N)
data.frame(bet, fill=outcome, ymin, ymax)
})
require(ggplot2)
p <- ggplot(out, aes(xmin=bet-0.5, xmax=bet+0.5, ymin=ymin, ymax=ymax)) +
geom_rect(aes(fill=fill), colour="grey80") +
scale_fill_gradient("Outcome", low="red", high="blue") +
xlab("Bet")
The plot:
编辑：解释你的旧代码是如何Andrie工作的，以及为什么它没有给出你想要的。
基本上，安德烈所做的（或者更确切地说是一种看待它的方式）是使用这样的想法，即如果你有两个二项分布，X ~ B(n, p)并且Y ~ B(m, p)，其中n, m = size和p = probability of success，那么它们的总和，X + Y = B(n + m, p)（1）。所以，目的xcum是为了得到所有n = 1:32折腾的结果，但为了更好地解释，让我一步一步构造代码。伴随着解释， for 的代码xcum也将非常明显，并且可以立即构建（无需每次都for-loop构建一个。cumsum
如果您到目前为止一直关注我，那么，我们的想法是首先创建一个numtri * numbet矩阵，每列 ( length = numtri) 分别具有0's和1's概率 =5/6和1/6。也就是说，如果你有numtri = 1000，那么，你将有 ~ 8340's和 166 1's* 为每一numbet列（= 32 这里）。让我们先构建它并测试它。
numtri <- 1e3
numbet <- 32
set.seed(45)
xcum <- t(replicate(numtri, sample(0:1, numbet, prob=c(5/6,1/6), replace = TRUE)))
# check for count of 1's
> apply(xcum, 2, sum)
[1] 169 158 166 166 160 182 164 181 168 140 154 142 169 168 159 187 176 155 151 151 166
163 164 176 162 160 177 157 163 166 146 170
# So, the count of 1's are "approximately" what we expect (around 166).
n = 1现在，这些列中的每一列都是具有和的二项分布样本size = numtri。如果我们将前两列相加并用这个总和替换第二列，那么，从 (1) 开始，由于概率相等，我们最终会得到一个二项式分布n = 2。同样，相反，如果您添加了前三列并用这个总和替换了第 3 列，您将获得二项式分布，n = 3依此类推...概念是，如果您cumulatively添加每一列，那么您最终会得到numbet二项分布的数量（此处为 1 到 32）。所以，让我们这样做。
xcum <- t(apply(xcum, 1, cumsum))
# you can verify that the second column has similar probabilities by this:
# calculate the frequency of all values in 2nd column.
> table(xcum[,2])
0 1 2
694 285 21
> round(numtri * dbinom(2:0, 2, prob=5/6))
[1] 694 278 28
# more or less identical, good!
如果你划分xcum，到目前为止，我们cumsum(1:numbet)以这种方式在每一行上生成：
xcum <- xcum/matrix(rep(cumsum(1:numbet), each=numtri), ncol = numbet)
这将与产生的xcum矩阵相同for-loop（如果您使用相同的种子生成它）。但是，我不太明白 Andrie 进行这种划分的原因，因为这对于生成您需要的图表不是必需的。但是，我想这与您在之前关于 math.stackexchange 的帖子中frequency谈到的值有关
现在谈谈为什么你很难获得我附上的图表（带n+1垃圾箱）：
n=1:32对于带有试验的二项分布，5/6作为尾部概率（失败）和正面1/6概率（成功），正面概率由k下式给出：
nCk * (5/6)^(k-1) * (1/6)^k # where nCk is n choose k
对于我们生成的测试数据，对于n=7和（试验），和正面n=8的概率由下式给出：k=0:7k=0:8
# n=7
0 1 2 3 4 5
.278 .394 .233 .077 .016 .002
# n=8
0 1 2 3 4 5
.229 .375 .254 .111 .025 .006
为什么他们都有 6 个箱子而不是 8 个和 9 个箱子？当然，这与 的值有关numtri=1000。让我们通过直接从二项分布生成概率来了解这 8 个和 9 个 bin 中的每一个的概率是多少，dbinom以了解为什么会发生这种情况。
# n = 7
dbinom(7:0, 7, prob=5/6)
# output rounded to 3 decimal places
[1] 0.279 0.391 0.234 0.078 0.016 0.002 0.000 0.000
# n = 8
dbinom(8:0, 8, prob=5/6)
# output rounded to 3 decimal places
[1] 0.233 0.372 0.260 0.104 0.026 0.004 0.000 0.000 0.000
你看到对应k=6,7和k=6,7,8对应n=7的概率n=8是 ~ 0。它们的价值非常低。这里的最小值5.8 * 1e-7实际上是 ( n=8, k=8)。这意味着如果您进行1/5.8 * 1e7多次模拟，您有机会获得 1 个值。如果对 进行相同检查n=32 and k=32，则值为1.256493 * 1e-25。因此，您必须模拟那么多值才能获得至少 1 个结果，其中所有32结果都指向n=32.
这就是为什么您的结果没有某些 bin 的值的原因，因为对于给定的numtri. 出于同样的原因，直接从二项分布生成概率克服了这个问题/限制。
我希望我已经设法写得足够清晰，以便您遵循。让我知道您是否遇到困难。
编辑2：
当我模拟我刚刚在上面编辑过的代码时numtri=1e6，我得到了n=7and并计算了andn=8的正面数：k=0:7k=0:8
# n = 7
0 1 2 3 4 5 6 7
279347 391386 233771 77698 15763 1915 117 3
# n = 8
0 1 2 3 4 5 6 7 8
232835 372466 259856 104116 26041 4271 392 22 1
请注意，对于 n=7 和 n=8，现在有 k=6 和 k=7。此外，对于 n=8，k=8 的值为 1。随着增加numtri，您将获得更多其他丢失的垃圾箱。但这需要大量的时间/内存（如果有的话）。
