乱数によるシミュレーションで中心極限定理を確かめる
統計の初学者としての感想なのですが、統計モデルを「乱数発生器」と見なす考え方は学習の初期でしっかり教えてほしかったなと思っております。私が最初に勉強した『心理統計学の基礎』という本にも、よくみたら「確率モデル」とは何なのかを説明する節で、
母集団を一種のデータ発生装置とみなし、そこから確率的に発生したデータによってサンプルが構成されると考える
と書いてあります。
今にして思えばべつにふつうのことを言っているのですが、最初はこの「データ発生装置」「乱数発生器」のイメージをしっかり持つというのが、統計学を理解する第一歩になるというのをあまり分かってなかった気がします。
たとえばRだと、モデル式を書いて乱数を発生させるコマンドと組み合わせれば、そのモデルを「データ発生装置」として自分でサンプルデータを作り、数値実験を行ったりできます。「乱数」を実際にコンピュータで発生させていろいろシミュレーションしてみるというのにな慣れてくると、理解が進むというか、「イメージ」がつきやすくなってきますね。いろんな面で。
さて今日は、そういう簡単なシミュレーションの練習として、「モデル」をつくってどうのこうのとは違うのですが、中心極限定理は果たしてほんとなのかを自作データで確認してみようと思いました。これって初学者の定番ネタのひとつなんじゃないでしょか。
中心極限定理というのは、数学的な導出については理解してないのですが（汗）、
正規分布と中心極限定理
http://www.kisc.meiji.ac.jp/~nino/2010/1020.html
「平均
で標準偏差
の母集団から十分大きい標本数
の標本を抽出して平均
を求める。同じ標本数の標本を幾つも抽出し、それぞれの標本平均をデータとして集めると、母集団の分布とは関係なくそのデータは、平均
の正規分布に近づく。これを中心極限定理という」
というやつです。とにかく「平均値」を議論している限り、母集団の分布がどんなへんな形であろうが、ランダムサンプリングされた標本の平均値は何回もサンプリングを繰り返すと正規分布にしたがって変動するという話です。
Rでやってみる
まず、ガンマ分布を２つ重ねて、ヘンな分布の乱数を100万個つくってみます。
平均と標準偏差とヒストグラムも表示させます。
> pop <- c((rgamma(n=500000, shape=3, scale=20)*3), (1000-(rgamma(n=500000, shape=12, scale=10))*2)) > mean(pop) [1] 470.1164 > sd(pop) [1] 303.0091 > > hist(pop)
平均470.1164、標準偏差303.0091のヘンな乱数セットをつくることができました。*1
この母集団から、サンプリングを繰り返して平均値を記録していきたいのですが、そのためにまず関数をつくっておきます。
> cen.lim1 <- function (population, size, repeatation) { + # populationには乱数で作った母集団ベクトルを入れる + # sizeには母集団から何個サンプリングするかを指定 + # repeatationはサンプリングを何回繰り返すかを指定 + means <- c() + for (i in 1:repeatation) { + samp <- sample(pop, size) + means <- c(means, mean(samp)) + } + + return(means) + }
ここで、さっきのヘンな分布の母集団からのサンプリングを行ってみます。
サンプルサイズ500のランダムサンプリングを1000回繰り返して、1000回分の平均値がどのように分布するかをみてみます。
> result1 <- cen.lim1(population=pop, size=500, repeatation=1000) > > mean(result1) [1] 470.3386 > sd(result1) [1] 13.55627 > sd(pop)/sqrt(500) # 中心極限定理から予想されるsd [1] 13.55098
平均値はほぼ母集団平均の470.1164に近く、標準偏差も、中心極限定理から予想される値になってますね。
ヒストグラムみてみましょう。ついでに正規分布のグラフを重ねておきます。表示を整えるために作図関数にごちゃごちゃ書いてます。
> x <- seq(from=min(result1), to=max(result1)) > par(mar = c(5, 5, 5, 5)) > hist(result1, breaks=20, xlim=c(min(result1), max(result1)), xlab="Mean", main="変な分布からサンプリング", cex.main=0.9) > par(new=TRUE) > plot(x, dnorm(x, mean=mean(result1), sd=sd(result1)), type="l", lwd=1.5, col="blue", xlim=c(min(result1), max(result1)), xlab="", ylab="", axes=FALSE) > axis(4) > mtext("density", side = 4, line=3)
たしかに、もとのデータの分布がかなりヘンな割に、平均値のサンプリング結果は綺麗な正規分布になってますね。
ついで
次に、ついでなので、サンプルサイズをいろいろ変えた場合の結果を一覧にしてみます。
まぁ、関数つくる練習って感じです……。
> cen.lim2 <- function (population, size, repeatation) { + # populationには乱数で作った母集団ベクトルを入れる + # sizeには母集団から何個サンプリングするかをベクトルで指定 + # ※c(100, 500, 1000)とすると、サンプル10個、500個、1000個の場合の結果が出る + # repeatationにはサンプリングを繰り返す回数を指定 + + mmeans <- c() + sdmeans <- c() + normality <- c() + for (i in size) { # 指定したサンプルサイズごとに実行 + + means <- c() + for (j in 1:repeatation) { # サンプリングを指定した回数繰り返し + samp <- sample(population, i) + means <- c(means, mean(samp)) + } + + mmeans <- c(mmeans, mean(means)) + sdmeans <- c(sdmeans, sd(means)) + normality <- c(normality, ks.test(means, "pnorm", mean=mean(means), sd=sd(means))$p.value) + } + + # 結果の出力テーブル + result <- data.frame( + Size=size, # 指定したさまざまなサンプルサイズ + Mean=mmeans, # 繰り返し取得したサンプル平均の平均値 + SD=sdmeans, # 繰り返し取得したサンプル平均の標準偏差 + Normality=normality, # 正規性検定結果（p値） + PopMean=rep(mean(population), length(size)), # 母集団平均 + PopSD=rep(sd(population), length(size)), # 母集団のsd + EstimatedSD=(sd(population)/sqrt(size)) # sdの理論的な推定値 + ) + + return(result) + }
こういう関数をいったんつくっといて……
> result2 <- cen.lim2(population=pop, size=c(5, 10, 100, 500, 1000), repeatation=1000) > result2 Size Mean SD Normality PopMean PopSD EstimatedSD 1 5 473.41 134.3656 0.073836 470.12 303.01 135.510 2 10 470.93 93.3876 0.862887 470.12 303.01 95.820 3 100 472.54 30.4350 0.505786 470.12 303.01 30.301 4 500 470.20 13.8713 0.691897 470.12 303.01 13.551 5 1000 470.44 9.3685 0.945892 470.12 303.01 9.582
なんか、options(digits=5)としているのですが、表示桁数が項目によって異なってしまい、正規性の検定結果のところがすごい少数になってますが（解決策不明）。
平均値はまぁサンプルが大きい方が微妙に真の値に近づいていて、標準偏差もほぼ予想通りの値になっています。
ちなみに、正規性の検定は、データ数がそこそこ大きいのでシャピロ・ウィルク検定ではなくコルモゴロフ・スミルノフ検定の結果を（あまり理解せずに）表示しています。
正規性の検定についてはあとでまとめて勉強して理解しておかないとやばいなと思ってるのですが、さしあたり、
「統計学関連なんでもあり」の過去ログ--- 043
統計勉強会 - 東邦大学・理論生態 ←ここの「正規性の検定（宮内）.pptx」
Rによるコルモゴロフ・スミルノフ検定
このへんを参考にしました。１つめのリンク先は注意喚起として大事で、２つめのリンクはそれぞれの検定方式の数式が解説されています。３つめのリンクはタイトルのとおりRでの実行方法。
参考に、次のエントリで、一様乱数と正規乱数に対してシャピロ・ウィルク検定とコルモゴロフ・スミルノフ検定を実行したときの結果をグラフにしておきます。
さらについでに、サンプルサイズではなくサンプリング回数を変えた場合の結果も一覧にしておきます。
だから何という感じではあります。
> cen.lim3 <- function (population, size, repeatation) { + # populationには乱数で作った母集団ベクトルを入れる + # sizeには母集団から何個サンプリングするかを指定 + # repeatationにはサンプリングを繰り返す回数をベクトルで入れる + # ※c(10, 100)とすると、繰り返し数10の場合と100の場合の結果が出る + # サンプリングを繰り返してそれぞれの平均値を取得 + # 平均値ベクトルの平均、標準偏差を得る + # 参考として、母集団のmean、SD、中心極限定理による推定SDを並べる + + mmeans <- c() + sdmeans <- c() + normality <- c() + for (i in repeatation) { # 指定した複数の繰り返し回数ごとに実行 + + means <- c() + for (j in 1:i) { # サンプリングを指定した回数繰り返し + samp <- sample(population, size) + means <- c(means, mean(samp)) + } + + mmeans <- c(mmeans, mean(means)) + sdmeans <- c(sdmeans, sd(means)) + normality <- c(normality, ks.test(means, "pnorm", mean=mean(means), sd=sd(means))$p.value) + } + + + # 結果の出力テーブル + result <- data.frame( + Trial=repeatation, # 指定したさまざまなサンプリング回数 + Mean=mmeans, # 繰り返し取得したサンプル平均の平均値 + SD=sdmeans, # 繰り返し取得したサンプル平均の標準偏差 + Normality=normality, # 正規性検定結果（p値） + PopMean=rep(mean(population), length(repeatation)), # 母集団平均 + PopSD=rep(sd(population), length(repeatation)), # 母集団のsd + EstimatedSD=(sd(population)/sqrt(size)) # sdの理論的な推定値 + ) + + return(result) + } > > result3 <- cen.lim3(population=pop, size=500, repeatation=c(10, 100, 1000, 10000)) > result3 Trial Mean SD Normality PopMean PopSD EstimatedSD 1 10 465.98 13.821 0.66877 470.12 303.01 13.551 2 100 470.60 14.018 0.77609 470.12 303.01 13.551 3 1000 469.98 13.768 0.94606 470.12 303.01 13.551 4 10000 470.12 13.537 0.82283 470.12 303.01 13.551 >
*1:sd()はほんとは不偏標準偏差を求める関数であり、ここは何かを推定する場面ではないのでおかしいかもですが、データ数が100万もあるのだからどうでもいいでしょう。
