#本文
統計で「自由度」とか「n-1で割る」などの言葉を聞いたことありませんか？
どんな時に使われるのか、どんな意味かご存知ですか？
また、「nが十分に大きければnで割っても問題ない」という言葉も聞いたことありませんか？
今日はそのあたりをまとめて理解してみます。
内容はある程度基本的な統計を知っている前提で書き込みます。
修正・指摘歓迎です。
#分散って何だっけ
分散とは偏差の二乗和の平均です。
偏差はデータとデータの平均の差です。
データすべての平均を使って各データとの偏差を計算して二乗した平均を求めます。
母集団の分散は母分散と呼びます。
variance = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu)^2
母集団の分散に対して、母集団からn個抽出した「標本データ内での分散」を標本分散と呼びます。
sample variance = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2
母分散は母平均からのすべてのデータのばらつき具合を表す値であり、
標本分散は抽出してきたいくつかの標本データを使って求めた平均と、標本データのばらつき具合を表す値である。
#推測統計の話
統計の分野「推測統計」は「測定できない母集団という集団の性質を知りたい」という事を数学的に実現しようとする分野です。
たとえば全人類の平均身長がどのくらいで、身長の上限から下限までどのような分布で存在しているか説明できますか？
全人類を調べることはできない。
なので調べられる範囲で調べて、データの性質から全人類の身長についてできる限り正確に推定しようとするのです。
母集団の分散が分からないので、標本の分散をなるべく母集団に近づけたいと考えるわけです。
そこでn-1で割る必要がでてきます。
今回は母集団が正規分布という例で考えます。
データは平均150分散81(標準偏差9)を創り出して使います。
n=100000
norm_dist<-rnorm(n,150,sqrt(81))
hist(norm_dist)
正規分布といえば上図のような釣り鐘型の分布です。
上図のようなヒストグラムを見ると分かりますが、山の頂点に近いデータの密度が高くなっています。
正規分布の母集団から無作為に抽出してきた時、密度の高い部分から多くデータが出てきて、密度の低い両端からはデータが取り出しにくいことは想像できると思います。
つまり標本データを抽出すると、母平均の周辺に近い値が集まりやすいと言えます。
従って標本分散は母分散よりも小さくなってしまうわけです。
ではその標本分散の性質を修正してやりたいわけです。
乱暴な言い方をしますが、nで割るより少し値が大きくなるn-1で割った方が母分散に近くなりそうです。
後ほどもうすこし正確に理解していきますが、この修正した標本分散の事を不偏分散と呼びます。
途中で分からなくなるという方は不偏標本分散と読み替えてください。
#数値として確認する
確認の前に今回はR言語を使っていくのでR言語の分散を求める関数について小話
var(norm_dist)
[1] 81.04635
じつはRの分散を求める関数は既に不偏分散を求めるようになっています。
(復習：不偏分散は標本データの範囲内で母分散に近くなるように修正したものでした。)
#標本分散
sum((norm_dist-mean(norm_dist))^2)/n
[1] 81.04554
#不偏分散
sum((norm_dist-mean(norm_dist))^2)/(n-1)
[1] 81.04635
#実際に試して図で理解する
標本を取り出した時の標本分散と不偏分散のどちらが母分散に近くなるかを確認してみます。
標本5個の時
sample_n <- 5
sample_point <- norm_dist[1:sample_n]
#標本分散
sum((sample_point-mean(sample_point))^2)/sample_n
[1] 31.80251
#不偏分散
sum((sample_point-mean(sample_point))^2)/(sample_n-1)
[1] 39.75313
81という母分散からはまだまだかけ離れた値ですが、
不偏分散のほうが8ほど大きい値をとって標本分散よりも母分散に近くなっています。
10個とってきたら?
sample_n <- 10
sample_point <- norm_dist[1:sample_n]
#標本分散
sum((sample_point-mean(sample_point))^2)/sample_n
[1] 91.1652
#不偏分散
sum((sample_point-mean(sample_point))^2)/(sample_n-1)
[1] 101.2947
今度は不偏分散のほうが母分散から離れた値になってしまいました。
というのもランダムに生成したデータなので時にそのようなデータを抽出してしまうこともあります。
#じゃあ不偏分散のほうが母分散に近いとは必ずしも言えないのでは？？
これに関してはランダム抽出という試行を繰り返して、不偏分散と標本分散のどちらが比較的母分散に近いのかを確かめてみたいと思います。
標本のランダム抽出を1万回ほど行って、それぞれで標本分散と不偏分散を求めてみたいと思います。
抽出する標本の数も動かしてみます。
library("animation")
saveGIF({
for(i in 5:150){
sample_size<-i
norm_dist<-rnorm(sample_size*10000,150,sqrt(81))
sample_point <- data.frame(matrix(norm_dist,ncol=sample_size,nrow=10000))
sample_var_n <- rowSums((sample_point-rowMeans(sample_point))^2) / sample_size
sample_var_n_1 <- rowSums((sample_point-rowMeans(sample_point))^2) / (sample_size-1)
par(mfrow=c(1,2))
mode_n<-round(mean(sample_var_n))
plot(density(sample_var_n), xlim=c(-10,350), ylim=c(0,0.02),main = paste0("var divided by n mean = ", mode_n,"\n size = ",sample_size),xlab = "N = 10000 true var = 81")
abline(v=mode_n,col="red")
mode_n_1<-round(mean(sample_var_n_1))
plot(density(sample_var_n_1), xlim=c(-10,350), ylim=c(0,0.02),main = paste0("var divided by n-1 mean = ", mode_n_1,"\n size = ",sample_size),xlab = "N = 10000 true var = 81")
abline(v=mode_n_1,col="red")
}
}, interval = 0.5, movie.name = "sample_n_minus_1_5to100.gif")
右側がn-1で割った不偏分散
左側がnで割った標本分散
赤色の線が真の母分散の81です
不偏分散は最初から81に近い値で、標本分散は小さい値を示しています。
標本分散の方も40や50あたりから不偏分散と変わらないような値になってきています。
これが「nが十分に大きければ標本分散を使っても不偏分散と大きく変わらない」という言葉の意味になります。
#まだ納得いかないので証明する
上記で直観的にn-1で割った方が真の母集団の値に近くなることが分かっていただけたかと思いますが、n-1であってn-2やn-3でない理由について説明していこうと思います。
x_i = 標本i番目のデータ \\\
\bar{x} = 標本の平均 \\\
S^2 = 標本分散 \\\
\mu = 母平均 \\\
\sigma^2 = 母分散 \\\
とおいておく。
標本分散を書き出すと、
S^2 = \frac{1}{n}
\Bigl(
\bigl(
x_1 - \bar{x}
\bigr)^2
+ ・・・ +
\bigl(
x_n - \bar{x}
\bigr)^2
\Bigr)
となる。
ここでいささか突然だが、標本の平均について考える。
そもそも無作為抽出されたデータから求めた標本平均で、標本分散は計算されているのだから、標本平均も母平均とは異なっているはずである。
標本平均が、どのくらい母平均から離れているのかを求めておく。
E[(\bar{x} - \mu)^2] \\\
= \frac{1}{n^2}\sum_{i=1}^{n}(E[(x_i - \mu)^2]) \\\
= \frac{1}{n^2}\sum_{i=1}^{n}(\sigma^2) \\\
= \frac{\sigma}{n}
各標本データを母平均を使って分散を求めると、母分散に等しくなる。
シグマ記号の母分散の部分はiと関係なく、しかし母分散はn個あるのでn個の母分散となる。
分母のnの二乗と打ち消して最後の値が得られる。
話を戻して標本分散の式中に母分散をいれてやる。
というのも、標本分散を母分散を使って表現できれば、標本分散が母分散の何倍なのかを知ることができる。
s^2 = \frac{1}{n}
\Bigl(
\bigl(
x_1-\bar{x}-\mu+\mu
\bigr)^2
+ ・・・ +
\bigl(
x_n-\bar{x}-\mu+\mu
\bigr)^2
\Bigr)
書き換えて
s^2 = \frac{1}{n}
\Bigl(
\bigl(
(x_1-\mu)-(\bar{x}-\mu)
\bigr)^2
+ ・・・ +
\bigl(
(x_n-\mu)-(\bar{x}-\mu)
\bigr)^2
\Bigr)
これを展開公式を使って展開する。
前半の項と後半の項をひとつの塊として(A-B)^2として展開する。
s^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i-\mu)^2 - 2(\bar{x}-\mu)\frac{1}{n}\sum_{i=1}^{n}(x_i-\mu) + \frac{1}{n}n(\bar{x}-\mu)^2
2項目のシグマ記号は母分散からの標本平均の偏差である
1項目は二乗なのでおなじように修正はできない。
\frac{1}{n}\sum_{i=1}^{n}(x_i-\mu) = (\bar{x} - \mu)
よって2項目を書き直すと3項目と計算できるので
s^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i-\mu)^2 - (\bar{x}-\mu)^2
この式の期待値を求めると
s^2 = \frac{1}{n}\sum_{i=1}^{n}E[(x_i-\mu)^2] - E[(\bar{x}-\mu)^2]
1項目の値は母分散の式であり、
2項目の値は先ほど求めた母平均と標本平均のばらつきの期待値である。
置き換えてやると
s^2 = \sigma^2 - \frac{\sigma^2}{n} \\\
= \frac{n-1}{n}\sigma^2
ようやくここまできました。
標本分散は母分散のn分のn-1倍ということが分かりました。
よって標本分散を母分散と等しくするためには、標本分散にn-1分のn倍してやればいいのです。
\sigma^2
= \frac{n}{n-1}s^2 \\\
= \frac{n}{n-1} \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x}) \\\
= \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})
#おまけ
最頻値で表示してみました。
library("animation")
saveGIF({
for(i in 5:150){
sample_size<-i
norm_dist<-rnorm(sample_size*10000,150,sqrt(81))
sample_point <- data.frame(matrix(norm_dist,ncol=sample_size,nrow=10000))
sample_var_n <- rowSums((sample_point-rowMeans(sample_point))^2) / sample_size
sample_var_n_1 <- rowSums((sample_point-rowMeans(sample_point))^2) / (sample_size-1)
par(mfrow=c(1,2))
mode_n<-names(which.max(table(round(sample_var_n))))
plot(density(sample_var_n), xlim=c(-10,350), ylim=c(0,0.02),main = paste0("var divided by n mode = ", mode_n,"\n size = ",sample_size),xlab = "N = 10000 true var = 81")
abline(v=mode_n,col="red")
mode_n_1<-names(which.max(table(round(sample_var_n_1))))
plot(density(sample_var_n_1), xlim=c(-10,350), ylim=c(0,0.02),main = paste0("var divided by n-1 mode = ", mode_n_1,"\n size = ",sample_size),xlab = "N = 10000 true var = 81")
abline(v=mode_n_1,col="red")
}
}, interval = 0.5, movie.name = "sample_n_minus_1_5to100.gif")
