我已经使用R的dbinom函数生成了n=1:32次试验中正面朝上的频率，并绘制了图表。这是您所期望的结果。我阅读了您在SO和math.stackexchange上的一些早期帖子，但仍然不明白为什么您要“模拟”实验而不是从二项式随机变量中生成。如果您能解释一下，那就太好了！我将尝试使用@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")
情节：
编辑：解释你的旧代码为什么不能达成你的意图以及Andrie的代码是如何工作的。
基本上，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 次投掷的结果，但为了更好地解释它，让我逐步构造代码。随着解释的进行，xcum 的代码也将非常明显，可以在很短的时间内构建出来（无需任何必要的 for-loop 和每次构建 cumsum）。
如果你到目前为止都跟上了我的思路，那么我们的想法首先是创建一个 numtri * numbet 矩阵，每个列（length = numtri）都有 0's 和 1's，概率分别为 5/6 和 1/6。也就是说，如果你有 numtri = 1000，那么你将会有大约 834 个 0'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)))
> 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
现在，每个列都是二项分布的样本，其中
n = 1，
size = numtri。如果我们将前两列相加并用这个总和替换第二列，则根据(1)，由于概率相等，我们最终将得到一个
n = 2的二项分布。同样地，如果您将前三列相加并将第3列替换为此总和，那么您将获得一个
n = 3的二项分布，以此类推... 概念是，如果您
累计添加每列，那么您将得到
numbet个二项分布（这里是1到32）。所以，让我们这样做。
xcum <- t(apply(xcum, 1, cumsum))
> table(xcum[,2])
0 1 2
694 285 21
> round(numtri * dbinom(2:0, 2, prob=5/6))
[1] 694 278 28
如果您按照以下方式对每行使用 cumsum(1:numbet)，将我们迄今为止生成的 xcum 进行除法：
xcum <- xcum/matrix(rep(cumsum(1:numbet), each=numtri), ncol = numbet)
这将与for-loop生成的xcum矩阵完全相同（如果您使用相同的种子生成它）。然而，我不太理解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:7和k=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的值有关。通过使用
dbinom从二项分布直接生成概率，让我们看看每个这些8和9个箱子的概率，以了解为什么会发生这种情况。
dbinom(7:0, 7, prob=5/6)
[1] 0.279 0.391 0.234 0.078 0.016 0.002 0.000 0.000
dbinom(8:0, 8, prob=5/6)
[1] 0.233 0.372 0.260 0.104 0.026 0.004 0.000 0.000 0.000
您会发现，与
n=7 和
n=8 相关的
k=6,7 和
k=6,7,8 的概率约为
0。它们的值非常低。这里的最小值实际上是
5.8 * 1e-7（
n=8，
k=8）。这意味着，如果您模拟了
1/5.8 * 1e7 次，您就有机会获得 1 个值。如果您检查
n=32 和
k=32 的情况，则该值为
1.256493 * 1e-25。因此，您将需要模拟那么多的值，才能获得至少一个结果，其中所有
32 种结果都是正面朝上的，适用于
n=32。
这就是为什么您的结果在某些箱子中没有值的原因，因为对于给定的
numtri，其发生的概率非常低。出于同样的原因，直接从二项分布生成概率可以克服这个问题/限制。
我希望我写得足够清楚，让您能够理解。如果您阅读有困难，请告诉我。
编辑 2： 当我用
numtri=1e6 模拟上面编辑的代码时，对于
n=7 和
n=8，并计算
k=0:7 和
k=0:8 的正面朝上的次数，我得到了以下结果：
0 1 2 3 4 5 6 7
279347 391386 233771 77698 15763 1915 117 3
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的增加，您将获得更多其他缺失的bin。但这将需要大量的时间/内存（如果可能的话）。
