ggplot 中的概率热图

I asked this question https://stackoverflow.com/questions/7305803/plot-probability-heatmap-hexbin-with-different-sized-bins a year ago and got code for this "probability heatmap":

numbet <- 32
numtri <- 1e5
prob=5/6
#Fill a matrix 
xcum <- matrix(NA, nrow=numtri, ncol=numbet+1)
for (i in 1:numtri) {
x <- sample(c(0,1), numbet, prob=c(prob, 1-prob), replace = TRUE)
xcum[i, ] <- c(i, cumsum(x)/cumsum(1:numbet))
}
colnames(xcum) <- c("trial", paste("bet", 1:numbet, sep=""))

mxcum <- reshape(data.frame(xcum), varying=1+1:numbet, 
idvar="trial", v.names="outcome", direction="long", timevar="bet")


library(plyr)
mxcum2 <- ddply(mxcum, .(bet, outcome), nrow)
mxcum3 <- ddply(mxcum2, .(bet), summarize, 
            ymin=c(0, head(seq_along(V1)/length(V1), -1)), 
            ymax=seq_along(V1)/length(V1),
            fill=(V1/sum(V1)))
head(mxcum3)

library(ggplot2)

p <- ggplot(mxcum3, aes(xmin=bet-0.5, xmax=bet+0.5, ymin=ymin, ymax=ymax)) + 
geom_rect(aes(fill=fill), colour="grey80") + 
scale_fill_gradient("Outcome", formatter="percent", low="red", high="blue") +
scale_y_continuous(formatter="percent") +
xlab("Bet")

print(p)

（可能需要稍微更改此代码，因为this https://stackoverflow.com/questions/10146109/formatter-argument-in-scale-continuous-throwing-errors-in-r-2-15)

This is almost正是我想要的。除了每个垂直轴应具有不同数量的垃圾箱外，即第一个应有 2 个，第二个应有 3 个，第三个应有 4 个（N+1）。在图表中，轴 6 +7 具有相同数量的箱 (7)，其中 7 应具有 8 (N+1)。

如果我是对的，代码这样做的原因是因为它是观察到的数据，如果我进行更多试验，我们会得到更多的垃圾箱。我不想依靠试验次数来获得正确的垃圾箱数量。

我如何调整此代码以提供正确的垃圾箱数量？

我用过Rdbinom生成头部的频率n=1:32现在进行试验并绘制图表。这将是你所期望的。我读过你之前关于 SO 和 on 的一些帖子math.stackexchange。我还是不明白你为什么想要simulate实验而不是从二项式 R.V. 生成如果你能解释一下，那就太好了！我将尝试使用 @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:

Edit:解释你的旧代码是如何产生的Andrie有效以及为什么它没有达到您的预期。

基本上，安德烈所做的（或者更确切地说是一种看待它的方式）是使用这样的想法：如果你有两个二项式分布，X ~ B(n, p) and Y ~ B(m, p), where n, m = size and 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 and 1's概率=5/6 and 1/6分别。也就是说，如果你有numtri = 1000第834章0's和 1661'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 and size = numtri。如果我们将前两列相加，并用该总和替换第二列，那么，从 (1) 开始，由于概率相等，我们最终将得到二项式分布n = 2。同样，如果您添加了前三列并用该总和替换了第三列，您将获得二项式分布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 进行这种划分的原因，因为这对于生成您需要的图表来说不是必需的。不过，我认为这与frequency你谈到的价值观在先前关于 math.stackexchange 的文章中 https://math.stackexchange.com/questions/37655/calculate-number-of-sequences-in-frequency-matrix

现在谈谈为什么你很难获得我所附的图表（带有n+1 bins):

对于二项式分布n=1:32试验，5/6作为尾部（失败）的概率和1/6作为正面（成功）的概率，k头数由下式给出：

nCk * (5/6)^(k-1) * (1/6)^k # where nCk is n choose k

对于我们生成的测试数据，对于n=7 and n=8（试验），概率k=0:7 and 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。让我们通过使用以下命令直接从二项式分布生成概率来看看这 8 个和 9 个箱中每个箱的概率是多少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 and k=6,7,8对应于n=7 and n=8 are ~ 0。它们的价值非常低。这里的最小值是5.8 * 1e-7实际上 （n=8, k=8）。这意味着如果您模拟的话，您有机会获得 1 值1/5.8 * 1e7次。如果您检查相同的n=32 and k=32，值为1.256493 * 1e-25。因此，您必须模拟那么多值才能获得至少 1 个结果，其中所有值32结果是走向n=32.

这就是为什么您的结果没有某些箱的值，因为对于给定的箱来说，具有该值的概率非常低numtri。出于同样的原因，直接从二项分布生成概率克服了这个问题/限制。

我希望我已经写得足够清楚，以便您能够理解。如果您遇到困难请告诉我。

Edit 2:当我模拟上面刚刚编辑的代码时numtri=1e6，我得到这个n=7 and n=8并计算出头的数量k=0:7 and k=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您将获得更多其他丢失的垃圾箱。但这需要大量的时间/内存（如果有的话）。