二項分布
二項分布 (binomial distibution) とは
「質的データで二値（仮にAとBとする）の事象をN回おこなったとき、そのうち何回がAなのか」という状況を表した確率分布を二項分布といいます。例えば、バスケットボールのゴールにボールを10回投げたとき、そのうち何回入ったか」みたいな状況の確率分布ですね。半々の確率（0.5の確率）で成功する人も、ときには全く入らないかもしれないし（10回中0回）、全部入るかもしれません。でも、10回中5回入る確率が一番高そうですね。
二項分布は、数式としては、Aとなる割合を母比率pとしたときに、一度もAにならない (Aが0回の) 確率、Aが1回の確率、、、、N回すべてでAとなる確率によって描かれる確率を表現した確率関数です。そして、その確率関数を、X軸を回数、Y軸を確率という図に示すと「二項分布の山」が描かれます。
例えば、表が出るという母比率0.5のコインを３回投げたとき、表が０回、１回、２回、３回でる確率は、母比率 p = 0.5、データ数 N = 3の二項分布です。表が０回となる確率は、組み合わせが１通りしかないので、0.53ですが、表が１回となる確率は、３回投げたうちのどこか１回が表になるわけですから、組み合わせが３通りあります。そのため、確率は 3 x 0.53です。同じように、表が２回となる確率も3 x 0.53となり、３回とも表である確率は0.53となります。
これらの計算結果をグラフにして、例えばX軸の下限を-1、上限を4にすると、-1や4になる確率は０なので、この二項分布が「山」のように見えることが納得できることでしょう。ちなみに、どんな確率分布も、すべての取りうる確率を足し合わせたら１になります。つまり「山」の面積は１です。これはベルヌーイ分布でも同じです。
二項分布の確率関数
二項分布の確率関数は以下の形で表されます。
f(x) = nCxpx(1-p)n-x
ここでnCxとは組み合わせの数です。上の例ではp = 0.5だったため、1-p = 0.5となり、0.53が常に含まれる形となっていました。
Rで二項分布を扱う
では、二項分布をRで描いてみましょう。まずは前知識です。Rにおける確率分布のコマンドは頭文字１字で「Rがすべき作業内容」を表し、２文字目から確率分布の略称となっていることが多いです。二項分布の英名はbinomial distributionなので「binom」が二項分布の略称です。
頭文字1字には、d、p、rの３通りがあります。
dbinom: 想定した二項分布におけるX軸の値ごとの「確率」を求める。densityのd
pbinom: 想定した二項分布において、X軸が、ある値以下となる(orより大きな値となる)「確率」を求める。probabilityのp
rbinom: 想定した二項分布から、X軸の値をランダムに選択する。randomのr
ある値となる確率
ある二項分布における各値の確率を計算したい場合はdbinomを用います。
母比率probが0.5で試行回数sizeが10の二項分布で、xが3になる確率密度を求めましょう。
dbinom(3, size=10, prob=0.5)
別の値も試してみましょう。
dbinom(5, size=10, prob=0.5)
dbinom(7, size=10, prob=0.5)
dbinom(10, size=10, prob=0.5)
ある値以下 (よりも大きな値) となる確率
では、値が3以下になる確率はなんでしょうか。
pbinom(3, size=10, prob=0.5)
?pbinom #このヘルプ画面を読むと、pbinomではlower.tail=TRUEになっています。
つまり、ここでは、x=3, 2, 1, 0となる確率の和を計算しています。
dbinom(3, size=10, prob=0.5)+dbinom(2, size=10, prob=0.5)+dbinom(1, size=10, prob=0.5)+dbinom(0, size=10, prob=0.5)
lower.tail=FALSEにすると、Xが3よりも大きな値となる確率を計算します。
pbinom(3, size=10, prob=0.5, lower.tail=F)
pbinom(3, size=10, prob=0.5) + pbinom(3, size=10, prob=0.5, lower.tail=F)
両方を足すと１になりましたね。どんな確率分布でも、確率分布の総面積は１です。
二項分布の作図
まずはbarplotを使って、二項分布を描いてみましょう。
母比率がデータ数（試行回数, size）が5と10の二項分布も描きます。
x <- 0:5 #まず0, 1, 2, ...5の数列を作り、それをxに代入します。
y <- dbinom(x, size=5, prob=0.5) #xの各値の確率を求めて、yに代入します。
barplot(y, col="darkgreen")
今度はデータ数が10の場合です。
x <- 0:10
barplot(dbinom(x, size=10, prob=0.5), col="coral", ylim=c(0, 0.25))
ylimを付けたのは、付けない呪文だと0.25を表示してくれなかったからです。
次にplotを使って母比率0.5でデータ数10の二項分布を折れ線グラフで描きます。
barplotと比べて、折れ線グラフのほうが複数の分布を比べやすいからです。
x <- 0:10 #0から10までの整数の数列をつくります。
このxの値ごとで確率を計算して、それを線で結びます (type="l")。
なお、type="l"のlはイチではなくエルの小文字です。
plot(x, dbinom(x, size=10, prob=0.5), type="l", ylim=c(0, 0.5))
いちど作図されたあとで線を描き足すにはlinesを使います。
母比率を変えて重ね書きしてみましょう。
lines(x, dbinom(x, size=10, prob=0.2), col="red")
lines(x, dbinom(x, size=10, prob=0.9), col="blue")
二項分布は母比率が0あるいは1に近いほど左右非対称な形になります。
今度はデータ数を増やして、母比率probが0.5でデータ数（試行回数, size）が100の二項分布を描きましょう。
x <- 0:100 #0から100までの整数の数列をつくります。
plot(x, dbinom(x, size=100, prob=0.5), type="l", ylim=c(0, 0.2))
lines(x, dbinom(x, size=100, prob=0.2), col="red")
lines(x, dbinom(x, size=100, prob=0.9), col="blue")
このように、二項分布はデータ数が多いほど左右対称な形になります。
二項分布からのランダムなサンプリング
データ数が10、母比率が0.5の二項分布からランダムに10個の数字をサンプリングします。
rbinom(10, size=10, prob=0.5) -> x
hist(x)
母集団の母比率は0.5でも、サンプリングした数が少ないと、必ずしも左右対称な図になりません。
もし幸運にも（？）キレイな左右対称な山になっている人は、上の呪文をもう一度唱えてみてください。
サンプリングの数が少ないと、母集団の分布の「山」が必ずしもうまく反映されないのです。
今度はランダムに1000個の数字をサンプリングします。rbinom(1000, size=10, prob=0.5) -> x1
hist(x1, , breaks=0:10)rbinom(1000, size=10, prob=0.5) -> x1
hist(x1, , breaks=0:10) #ほぼ左右対称な二項分布になったはずです。
それでは、母比率を0.1にしてみましょう。
rbinom(1000, size=10, prob=0.1) -> x2
hist(x2, breaks=0:10)
今度は左右非対称になりましたね。
このように、二項分布は母比率とデータ数（試行回数）で形が決まります。
二項分布の最尤推定
二項分布では、データ数は明らかなので、母集団推定するのは母比率だけです。母比率を推定するには、いろいろな母比率の二項分布で、データとまったく同じものが得られる確率（尤度）が最大になる母比率を選びます。ここで、じつは二項分布の最大尤度の求め方はベルヌーイ分布と同じになります。実際にやってみれば分かります。母比率は表が出る割合 p で、データ数が5の二項分布を考えてみましょう。母比率 p を推定します。
５回のうち、表、表、表、裏、裏となった場合を考えてみましょう。これら各データが得られる確率は順に、p, p, p, (1-p), (1-p)ですから、尤度はp3(1-p)2となります。このように、二項分布の最尤推定では、組み合わせの数が登場しないため、ベルヌーイ分布と同じになるのです。
結論として、二項分布の母比率を最尤推定すると、そのデータにおける母比率の対象となる事象の割合が、母比率の最尤推定値となります。それはx/n、上の場合であれば3/5 = 0.6が母比率です。