読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第14回 推測統計学(標準誤差、中心極限定理)

単変量解析

第14回は推測統計学について書きます。

 

本来はSPSSで行っていた分析をRで再現するためのコードを書くという趣旨のはずだったのですが、Rを使っていると本当にいろんなことができるので、忘れないうちにいろいろ書いておきます。SPSSでは中心極限定理のシュミレーションはできませんし...できるのかな... 

a. 標本抽出

簡単におさらい。母集団から標本抽出して標本を得る。次に、標本から標本平均、標本分散、標本相関係数などの標本統計量を算出し、それを用いて母数(母平均、母分散など)を推定します。つまり、母集団→標本抽出→標本統計量算出→母数推定、という流れ。

※サンプル数:

標本抽出を行った回数のことで、サンプルサイズとは違うので、混同しないよう注意。サンプルサイズは1標本におけるnの数。

例えばn=1000などと表記されている。サンプル数=1000とは表記しない。

※母数:

母集団の特徴を表す数値であり、時々、分数の分母を「母数」という人がいますが、統計学的に母数は、母集団の確率分布を特徴付ける数なので注意。

b. 標本抽出誤差

抽出することで重要なことは「どのくらい母数に近い値を得られるか」「そのくらい信用できるか」ということですが、標本抽出によって得られた標本統計量と母数はほぼ一致しません(幸運の持ち主は別)。つまり、標本抽出で得られた値には誤差が伴います。そして、繰り返し標本抽出していると標本誤差も分布します。

ここからは、山田、杉澤、村井(2013)「Rによるやさしい統計学 」、オーム社のコードを参考に、正規母集団の母平均を推定していきます。

まず、母平均を推定するといっても、母集団は神様でない限りわかりません。そのため、母集団を仮定してから実験してみましょう。言い換えると「もし母集団がこのような正規分布のとき、このくらいあてになるという推定値がえられる」ということをみていきましょう。

sample <- rnorm(10,50,10) 
#average50、sd10の正規分布から10個の乱数を発生
sample
[1] 67.06154 59.00746 51.77259 50.16063 55.37940 44.02646 34.46429 34.92906 59.44837 55.32392
mean(sample)
[1] 51.15737 #およそ1.1の標本誤差
シュミレーション
sample.average <- numeric(length=10000)#推定値を入れる箱を10000個準備
> for (i in 1:10000){ #中括弧内の処理を10000回繰り返す
+  sample <- rnorm(10,50,10)
+   sample.average[i] <- mean(sample)
+ }

モンテカルロシュミレーション部分の解説、中括弧の所です。

まず、10000回分の推定値を入れる箱を作成。それ以降では、各番号が割り付けられた箱の中に、平均50、標準偏差10の正規分布から発生する乱数の平均を求めて、それを1〜10000の箱に入れていく。

hist(sample.average)

f:id:teruaki-sugiura:20150712051715p:plain

標本抽出を10000回繰り返したときの標本平均の分布

そして、誤差の接待値が5以下の回数を調べる。つまり、母平均の推定値が45以上、55以下の回数を調べましょう。

absolute.value.under5<- ifelse(abs(sampl.average-50)<=5,1,0)
table(absolute.value.under5)
absolute.value.under5
   0   1
1092 8908

89%は誤差5以内(=母数)に収まっていることがわかりました。 

最後に標本分布の平均と分散はどうでしょうか。

mean(sample.average)
[1] 50.02495
mean((sample.average-mean(sample.average))^2)
[1] 9.938753

 仮定した正規分布と近い値が得られていることがわかります。

 

おまけ

中心極限定理のシュミレーションによる証明。

※一様分布
for(i in 1:10000){
+ 標本 <- runif(100,0,1)
+ 標本平均[i] <- mean(標本)
+ }
hist(標本平均)
f:id:teruaki-sugiura:20150712053426p:plainポアソン分布
for(i in 1:10000){
+ 標本 <- rpois(100,3)
+ 標本平均[i] <- mean(標本)
+ }
f:id:teruaki-sugiura:20150712053526p:plain

※Γ分布
for(i in 1:10000){
+ 標本 <- rgamma(100,5)
+ 標本平均[i] <- mean(標本)
+ }
hist(標本平均)
f:id:teruaki-sugiura:20150712053708p:plain

※対数正規分布
for(i in 1:10000){
+ 標本 <- rlnorm(100,0,1)
+ 標本平均[i] <- mean(標本)
+ }
hist(標本平均)
f:id:teruaki-sugiura:20150712053734p:plain

 

広告を非表示にする