R を使用して、試行dbinomのヘッドの頻度を生成しn=1:32、グラフをプロットしました。それはあなたが期待するものになります。SO および に関する以前の投稿をいくつか読みましたmath.stackexchange。それでも、二項式 RV から生成するのではなく、実験を行う理由がわかりません。simulate説明していただければ幸いです。@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 が行ったこと (または、むしろそれを見る 1 つの方法) は、2 つの二項分布がある場合、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-loopcumsum
ここまで説明してきたのであれば、最初にnumtri * numbet、各列 ( length = numtri)に確率 =0'sおよびそれぞれ確率 =をもつ行列を作成するという考えです。つまり、 がある場合、列ごとに~ 834と 166 * になります (ここでは =32)。これを構築して、最初にこれをテストしましょう。1's5/61/6numtri = 10000's1'snumbet
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。最初の 2 列を追加し、2 列目をこの合計に置き換えると、(1) から、確率が等しいため、 の二項分布になりn = 2ます。同様に、代わりに、最初の 3 つの列を追加し、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: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
両方とも 8 と 9 のビンではなく、6 つのビンを持っているのはなぜですか? もちろん、これは の値と関係がありますnumtri=1000。dbinomなぜこれが起こるのかを理解するためにを使用して、二項分布から直接確率を生成することにより、これらの 8 と 9 の各ビンの確率を見てみましょう。
# 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に対応する確率は ~であることがわかります。それらは値が非常に低いです。ここでの最小値は実際には ( , ) です。これは、何度もシミュレートした場合、1 つの値を取得する可能性があることを意味します。についても同様にチェックすると、値はです。したがって、その多くの値をシミュレートして、すべての結果が.n=7n=805.8 * 1e-7n=8k=81/5.8 * 1e7n=32 and k=321.256493 * 1e-2532n=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
現在、n=7 と n=8 に対して k=6 と k=7 があることに注意してください。また、n=8 の場合、k=8 の値は 1 です。増加numtriすると、他の不足しているビンをより多く取得できます。ただし、膨大な時間/メモリが必要になります (必要な場合)。
