在ggplot Probabilty热图(Probabilty heatmap in ggplot)

我问这个问题在一年前拿到代码为这个“概率热图”：

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)

（可能需要稍微改变，因为这个代码本）

这几乎正是我想要的。除了每个垂直轴应具有仓的不同数目，即，在第一应具有2秒3，第三4（N + 1）。在该曲线图轴6 7具有相同数量的二进制位（7），其中7应具有8（N + 1）。

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

我怎样才能适应这个代码给仓的是否正确？

我已经使用的r dbinom产生的磁头的频率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工作的，为什么它不给你打算什么。

基本上，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)))

# 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的原因，因为这是没有必要为你生成所需的图形。不过，我想这事做与frequency你谈到值在math.stackexchange先前的帖子

现在到为什么你有获得我已经连接了的图（有困难的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 。让我们来看看直接使用二项式分布概率产生什么这些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和k=6,7,8对应于n=7和n=8是〜 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 。出于同样的原因，直接从二项式分布产生的概率克服了这个问题/限制。

我希望我已经成功了足够的清晰度为你写跟随。让我知道，如果你的麻烦经历。

编辑2：当我模拟的代码我刚刚上面编辑numtri=1e6 ，我得到这个对于n=7和n=8和计数的用于头的数量k=0:7和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

需要注意的是，有k = 6的且k = 7现在对于n = 7和n = 8。此外，对于n = 8，则有对于k = 8的值为1。随着numtri你会获得更多的其他失踪箱。但它会需要巨大的时间/内存容量（如果有的话）。