カタカタブログ

SIerで働くITエンジニアがカタカタした記録を残す技術ブログ。Java, Oracle Database, Linuxが中心です。たまにRuby on Railsなども。

【統計】カイ二乗分布をR言語で実験してみた

再び、統計の勉強。今回はカイ二乗分布について勉強してみたので、分布をR言語で実際にデータを生成してヒストグラムを描いてみる。

カイ二乗分布の標準正規分布の二乗和として定義され、自由度dfごとに分布がある。数式での定義は他の文献に譲るとして、要するに標準正規分布からdf個のデータをサンプリングして、それらの二乗したものを足し合わせる、という操作を十分な回数だけ行えば、カイ二乗分布が得られるはず!ということで早速実験してみた。

※前回同様、今回の実験も以下の書籍に感化されてやっています。

統計学が最強の学問である[実践編]---データ分析のための思想と方法

統計学が最強の学問である[実践編]---データ分析のための思想と方法

カイ二乗分布の定義

まず教科書的なカイ二乗分布の定義を確認しておく。
Z1, Z2, … , Zkを独立な標準正規分布N(0,1)に従う確率変数とする。
このとき、
f:id:osn_th:20160409101016j:plain
となる確率変数χ2が従う確率分布を自由度kのカイ二乗分布と呼ぶ。

カイ二乗分布のヒストグラム

まず、自由度df、つまり標準正規分布からサンプリングするデータの数を引数として、その二乗和を計算するkai2という関数を定義する。

# 標準正規分布から自由度dfの数だけサンプルを抽出し、それらの二乗和を求める
kai2 <- function(df) {
  sum(rnorm(df, mean=0, sd=1)**2)
}

> kai2(1)
[1] 0.04196556
> kai2(1)
[1] 0.8841481
> kai2(2)
[1] 1.8568
> kai2(2)
[1] 0.8412345
> kai2(3)
[1] 7.518514
> kai2(3)
[1] 0.7282372

自由度1

続いて、このkai2を自由度1で10,000回呼び出して得られたデータから、ヒストグラムを描いてみる。なお、縦軸は確率密度(つまり、その値をとる確率)を表している。

# カイ二乗分布
x <- sapply(c(1:10000), function(x) { kai2(1)} )
hist(x, freq=FALSE, breaks=30)

すると以下のヒストグラムが得られる。
f:id:osn_th:20160402102211p:plain
さらに、比較のためR言語標準のカイ二乗分布を表す関数をこのヒストグラムの上に被せてみる。

curve(dchisq(x, 1), 0, 8, add=TRUE)

確かに、一致している!
f:id:osn_th:20160402102213p:plain

自由度もいろいろ変えてみる。

自由度2

# カイ二乗分布(自由度2)
df <- 2
x <- sapply(c(1:10000), function(x) { kai2(df)} )
hist(x, freq=FALSE, breaks=30, ylim=c(0, 0.5))
curve(dchisq(x, df), 0, 8, add=TRUE)

f:id:osn_th:20160402102216p:plain

自由度3

f:id:osn_th:20160402102219p:plain

自由度8

f:id:osn_th:20160402102221p:plain

自由度100

f:id:osn_th:20160402102224p:plain
自由度も大きくすると、山が次第に分布の中心に寄って行き、
100ほどになるとほぼ正規分布に近づいていることが分かる!

まとめ

カイ二乗分布が標準正規分布からサンプリングした値の二乗和の分布になる、ということが実験でより感覚的に確かめられた!また、自由度、つまりサンプリングするデータ数を増やしていくと正規分布に近づいていくことも確認できた!

以上。

※統計は今このあたりの本で勉強してます。

統計学が最強の学問である

統計学が最強の学問である

統計学が最強の学問である[実践編]---データ分析のための思想と方法

統計学が最強の学問である[実践編]---データ分析のための思想と方法

Rによるやさしい統計学

Rによるやさしい統計学

関連記事