【統計】二項分布をR言語で描画してみた
今回は二項分布のヒストグラムをR言語で可視化してみた。二項分布は「コイン投げで表が出る」や、「サイコロを振って1の目が出る」のような事象が一定確率で「成功するか、しないか」の2者択一の試行をそれぞれ独立に繰り返したときに現れる分布(ちなみにこういう試行をベルヌーイ試行という)。パラメータとして、試行回数と事象が成功する確率の2つをとりうるので、今回は様々な値を使って実際にヒストグラムを描いてみることにする。
二項分布の定義
まず二項分布の正確な定義を確認する。成功確率pの試行をn回繰り返したときに、x回成功する確率f(x)は、以下の式で求められる。ただし、x = 0, 1, 2, … , n。
このf(x)の確率分布が二項分布であり、B(n, p)で表す。
二項分布のヒストグラム
上記の定義を踏まえて、早速Rで実装してみる。まず、成功確率pの試行をn回繰り返したときの成功回数を求める関数を作成し、それを10,000回繰り返したときの分布(ヒストグラム)を描く関数を作成する。
# B(n,p)の二項分布に従う試行をしたときの成功回数 bin <- function(n, p = 0.5) { sum(sample(0:1, size = n, replace = TRUE, prob = c(p, 1 - p))) } # B(n,p)のヒストグラムを描画 hist_bin <- function(n = 20, p = 0.5) { x <- sapply(c(1:10000), function(x) { bin(n, p)} ) main <- paste("二項分布 B(", n, ", ", p, ")", seq="") hist(x, breaks=0:n, freq=FALSE, main=main, ylim=c(0,0.30)) }
まずは確率1/2の試行を20回行う場合の分布を描画する。また、例によってR標準の二項分布の関数を使って正しいグラフを重ねてみる。
hist_bin(20,0.5) curve(dbinom(x, 20, 0.5), 0, 20, 21, add=TRUE)
x=10を中心とした、正規分布のように左右対称な形のヒストグラムになった。ちなみに、二項分布の期待値E(X)はnとpの積で求められる。今回の場合だと20*1/2=10となっており、今回の実験結果と矛盾していない。
パラメータを変えながらヒストグラムを描画
続いて、二項分布のパラメータを徐々に変えながらいろいろなヒストグラムを描いてみる。
試行回数を固定し、確率を上げていく
まず、試行回数を20回で固定したまま、確率を0.1から0.9まで0.1ずつ増やしていく。
# B(n=20, p=[0.1,0.9]) for (i in seq(0.1, 0.9, 0.1)) { .file <- paste("out/plot", i, ".png", sep="") hist_bin(20, i) }
以下は描画した画像をgifアニメーションにまとめたもの。
同じ形の山が右から左に流れているのが分かる。
確率を固定し、試行回数を増やしていく
続いて、確率を1/2に固定したまま、試行回数を10~100まで10回ずつ増やしていった場合のヒストグラムを以下のようなコードで描いてみた。
# B(n=[10,100], p=0.5) for (i in seq(10, 100, 10)) { .file <- paste("out/plot", sprintf("%03d", i), ".png", sep="") hist_bin(i, 0.5) }
高い山が次第に小さくなっていくのが分かる。中心は期待値であるnとpの積になっている。
まとめ
二項分布をパラメータの値を変えながらヒストグラムを描いてみた。このような複数パラメータがある分布は教科書とかだと代表的なものを2,3個程度を重ねて描かれる程度だが、Rを使って連続実行することでパラメータの値を少しずつ変えながら表示できるので、動きが分かりやすい。
以上
参考文献
- 作者: 東京大学教養学部統計学教室
- 出版社/メーカー: 東京大学出版会
- 発売日: 1991/07/09
- メディア: 単行本
- 購入: 158人 クリック: 3,604回
- この商品を含むブログ (82件) を見る
- 作者: 山田剛史,杉澤武俊,村井潤一郎
- 出版社/メーカー: オーム社
- 発売日: 2008/01/25
- メディア: 単行本
- 購入: 64人 クリック: 782回
- この商品を含むブログ (68件) を見る
関連する記事