第７章 二項分布（授業ノート）

茶柱が立つと縁起がいいと言われます。統計学を学び始めたP子さんは茶柱の立つ確率を調べてみることにしました。毎朝一番茶を入れ、茶柱が立つかどうかをみるのです。カレンダーに丸をつけてもいいのですが、統計学の記号を使うことにしました。

• 正月、茶柱が立ちました　\(X_1=1\).
• 大晦日、茶柱は立ちませんでした　\(X_{365}=0\).

\(X=X_1+X_2+\cdots+X_{365}\)　は１年間で茶柱のたった日数を表します。
このように統計学では、０と１の値をとる変数\(X_i\)を用いて回数をカウントします。

1 ベルヌーイ分布.

茶柱の立つお茶はよく売れます。
あるお茶のメーカーのパックには「茶柱の立つ確率 \(\dfrac{1}{6}\) と記載されています。
P子さんは早速買ってきました。

物事を一般化するのが数学の方法です。そのほうが法則を見つけやすいからです。
茶柱の立つ確率を \(p\) とします。

　　　\(X_1=0,1\) 　確率変数は２通りの値しか取りません. 　　　

\(X_1\)	1	0
確率	\(p\)	\(1-p\)

このような分布をベルヌーイ分布といいます.

• 期待値　\(E(X_1)=p\).
  
• 分散　　\(V(X_1)=(1-p)^2p+(-p)^2(1-p)=p(1-p)\)

2 二項分布.

\(X_i\) は \(i\)日めに茶柱が立ったかどうかをカウントする確率変数です.

• \(X_1\) は１日目に茶柱がたった回数　(0 or 1).

• \(X_2\) は２日目に茶柱がたった回数　(0 or 1).
  ‥.

• \(X_n\) は n 日目に茶柱がたった回数　(0 or 1).

これらは互いに影響しない独立な確率変数です.
（確率変数の性質より、独立のときは分散は加算されます）.
\(X=X_1+X_2+\cdots+X_{n}\)　はn日間で茶柱の立った合計日数を表します。

• 期待値　\(E(X)=E(X_1)+E(X_2)+\cdots+E(X_{n})=np\).
  
• 分散　　\(V(X)=V(X_1)+V(X_2)+\cdots+V(X_{n})=np(1-p)\)

具体的な例でシミュレーションします.

2.1 コインを4回投げるとき、表の出る回数X.

• 表が出たら１.
• 裏が出たら０.

とし、１回目から４回目までの結果を x1,x2,x3,x4 とします。
起こりうる全ての場合を樹形図の考えで書き出すと１６通りになります。

x1 <- c(0,1)
x2 <- c(0,1)
x3 <- c(0,1)
x4 <- c(0,1)
data <- crossing(x1,x2,x3,x4)
data

x1	x2	x3	x4
0	0	0	0
0	0	0	1
0	0	1	0
0	0	1	1
0	1	0	0
0	1	0	1
0	1	1	0
0	1	1	1
1	0	0	0
1	0	0	1
1	0	1	0
1	0	1	1
1	1	0	0
1	1	0	1
1	1	1	0
1	1	1	1

表の出る回数 x1+x2+x3+x4 を表す列を付け加えると.

dat <- mutate(data,回数=x1+x2+x3+x4)
dat

x1	x2	x3	x4	回数
0	0	0	0	0
0	0	0	1	1
0	0	1	0	1
0	0	1	1	2
0	1	0	0	1
0	1	0	1	2
0	1	1	0	2
0	1	1	1	3
1	0	0	0	1
1	0	0	1	2
1	0	1	0	2
1	0	1	1	3
1	1	0	0	2
1	1	0	1	3
1	1	1	0	3
1	1	1	1	4

表の出る回数の昇順で並び替えます.

arrange(dat,回数)

x1	x2	x3	x4	回数
0	0	0	0	0
0	0	0	1	1
0	0	1	0	1
0	1	0	0	1
1	0	0	0	1
0	0	1	1	2
0	1	0	1	2
0	1	1	0	2
1	0	0	1	2
1	0	1	0	2
1	1	0	0	2
0	1	1	1	3
1	0	1	1	3
1	1	0	1	3
1	1	1	0	3
1	1	1	1	4

それぞれの回数が何通りあるかを集計すると、

table(dat$回数)

## 
## 0 1 2 3 4 
## 1 4 6 4 1

#dat %>% arrange(dat$回数)

この数値は　\({}_4 C_0,~ {}_4 C_1,~{}_4 C_2,~{}_4 C_3,~{}_4 C_4~\) であり二項係数と呼ばれます.
　　　　\({}_4 C_0+~ {}_4 C_1+~{}_4 C_2+~{}_4 C_3+~{}_4 C_4~=16\)　となっています.

\({}_4 C_i\) は 4個から \(i\) 個選ぶ組み合わせであり、Rでは関数chooseで計算することができます.

i <- 0:4
cho <- choose(4,i)
cho

## [1] 1 4 6 4 1

これがコインを4回投げるとき、表の出る回数 X の理論分布です。
　　　\(E(X)=4\cdot \dfrac{1}{2}=2~~~~~V(X)=4\cdot \dfrac{1}{2}\left(1-\dfrac{1}{2}\right)=1\).

赤は、同じ平均と分散を持つ正規分布\(N(2,1^2)\)です。
n=4 なのでヒストグラムは曲線とは言えませんが、二項分布が正規分布で近似できるのではないかと考えたのが ラプラスです.

2.1.1 半整数値補正.

正規分布の確率は標準正規分布表を利用できます。
二項分布に正規分布の近似を適用する場合は補正が必要です。
二項分布の確率はヒストグラムの長方形ですが、正規分布の確率はグラフとx軸で囲まれた面積だからです。
\(N(2,1^2)\) の確率密度関数を\(f(x)\) として、

• \(\displaystyle P(X=3)=\int_{2.5}^{3.5}f(x)dx\)　　　　\(\displaystyle P(2<X<3)=\int_{1.5}^{3.5}f(x)dx\).
  
• 一般には、　\(\displaystyle P(a<X<b)=\int_{a-0.5}^{b+0.5}f(x)dx\).

つまり、二項確率\(P(a<X<b)\) を求めるには、
　　　a-0.5<X<b+0.5　　として標準化し標準正規分布表を用いればよいことになります。
これを半整数値補正といいます。

2.2 コインを16回投げるとき、表の出る回数X.　

回数n が大きくなれば近似の精度が上がります。
（一般に、\(np>10\) が望ましいとされたいます）

i <- 0:16
x=seq(0,16,by=1)
X <- rep(x,choose(16,i))
table(X)

## X
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
##     1    16   120   560  1820  4368  8008 11440 12870 11440  8008  4368  1820 
##    13    14    15    16 
##   560   120    16     1

hist(X,breaks=seq(-0.5,16.5,by=1),prob=T)
curve(dnorm(x,8,2),add=T,co="red")

var(X)*(2^16-1)/2^16

## [1] 4

X は平均 8 、標準偏差 ２ の確率変数になります.
赤は、同じ平均と分散を持つ正規分布\(N(8,2^2)\)です。
近似精度が上がっていることが分かります。

ここまでは二項分布の理論確率であり、正規分布に近い分布になっています。
６章「自然現象、社会現象の多くは正規分布に従う」の例と考えることができます。

実際の実験では分布はどうなるのでしょうか。乱数を用いて確かめてみましょう。

3 研究.

二項分布は、「n回の試行でイベントが何回起こるか」という回数の確率分布です.
実際にコインを投げたりサイコロを降ったりするのは大変です。
乱数によるシミュレーションを行います.

rbinom(n , size , prob) は二項分布に従う乱数を生成します。

• n : 生成する乱数の数.
• size ：試行回数.
• prob : １回の試行でイベントが起こる確率.

3.1 　ベルヌーイ乱数　

rbinom(n=1 , size=1 , prob) は、
１回の試行で.
　　イベントが起きれば 1. 　　起きなければ 0.
　という回数をカウントする乱数を生成します.
この乱数をベルヌーイ乱数といいます.

3.1.1 コイン乱数.

乱数を用いて確率変数、
　　　「コインの表が出れば１、裏が出れば０」　を出力しましょう.

x <- rbinom(n = 1, size = 1, prob = 0.5)
x

## [1] 1

コインを４回投げる場合は.

x <- rbinom(n = 4, size = 1, prob = 0.5)
x

## [1] 0 0 1 1

これが4回、１セットの結果です。
このセットを手作業で繰り返すのは大変ですが、プログラミング「for 構文」なら容易に実行できます。

10セット繰り返してみます.

data=data.frame(x1=c(),x2=c(),x3=c(),x4=c(),計=c())
for (i in 1:30) {
  x <- rbinom(n = 4, size = 1, prob = 0.5)
  x <- c(x,sum(x))
  data <- rbind(data,x)
}
colnames(data) <- c("１回目","２回目","３回目","４回目","合計")
data

１回目	２回目	３回目	４回目	合計
0	0	0	0	0
1	1	1	1	4
0	1	1	1	3
1	1	0	1	3
1	0	1	1	3
1	0	1	1	3
1	0	0	0	1
0	1	0	0	1
1	0	0	0	1
0	1	1	0	2
1	1	1	1	4
0	1	1	1	3
0	0	1	0	1
1	0	0	0	1
0	1	0	0	1
1	0	0	0	1
0	0	0	1	1
1	1	1	1	4
1	0	0	1	2
1	1	1	0	3
1	0	0	1	2
0	0	1	0	1
1	1	0	0	2
1	1	1	0	3
1	0	0	1	2
0	1	1	0	2
0	0	0	1	1
1	1	1	1	4
0	0	1	0	1
1	1	1	0	3

par(family = "HiraKakuProN-W3") #日本語フォントの指定
hist(data$合計,prob=T,breaks=seq(-0.5,5.5,by=1) , xlab="x")

項目2.1は理論分布なので対称なヒストグラムになりましたが、実験では理想的なヒストグラムにはなりません。

3.2 適合度.

実験結果が理論分布に適合しているかどうかは§14の適合度検定で学習します。少し先走りますが chisq.test で確かめてみましょう。

chisq.test(table(data$合計), p = cho/16)

## Warning in chisq.test(table(data$合計), p = cho/16): Chi-squared approximation
## may be incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  table(data$合計)
## X-squared = 6.9333, df = 4, p-value = 0.1395

誤差変動で適合していないように見えますが、
帰無仮説「適合している」は棄却されません。
これくらいの誤差はありうるという判断です。

回数を増やしてヒストグラムを表示します。

data=data.frame(x1=c(),x2=c(),x3=c(),x4=c(),計=c())
for (i in 1:500) {
  x <- rbinom(n = 4, size = 1, prob = 0.5)
  x <- c(x,sum(x))
  data <- rbind(data,x)
}
colnames(data) <- c("１回目","２回目","３回目","４回目","合計")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
hist(data$合計,breaks=seq(-0.5,5.5,by=1),prob=T)
curve(dnorm(x,2,1),add=T)

回数が増えればほぼ正規分布に近くなります。この性質はラプラスの定理として知られています。

二項分布の正規近似条件について

4 例題７.

5 演習７.