統計学を勉強していくと、テキストによって、標準偏差や分散の分母がnだったり、n-1だったりして混乱しますよね。結論を言うと
- 目の前にある標本がどの程度ばらついているか表現したい場合はnで割る
- 標本から母集団の分散を「推定」する場合はn-1で割った方が推定精度が高い
ということです。この理由や使い分けについてモヤモヤしている方は、まずこちらをご覧ください。
このページでは、R言語を用いて、乱数を用いたシミュレーションを行い、その結果を交えながら、nなのかn-1なのかの説明をしています。
実験内容
平均0、標準偏差1の正規分布に従う乱数からn個のデータをサンプリング（\(x_i, i=1…n\)）して、以下の2つの分散（標本分散と不偏分散）を計算します。この実験をN回実施して両者の違いを確認してみます。
\(x\sim N(0, 1^2)\)
\(\bar{x}=\frac{1}{n}\Sigma_{i=1}^n{x_i}\)
標本分散： \(\sigma_1^2 = \frac{1}{n}\Sigma_{i = 1}^n(x_i-\bar{x})^2\)
不偏分散： \(\sigma_2^2 = \frac{1}{n-1}\Sigma_{i = 1}^n(x_i-\bar{x})^2\)
nについては、n = 3（違いが現れやすいケース）とn = 100（違いが現れにくいケース）の2ケース、繰り返し数Nは両者とも10万回としました。
結果
n = 3の場合
3個の乱数から標本分散、不偏分散を計算する作業を10万回行うプログラムをR言語で書いて試してみます。
ライブラリとしてtidyverseをロードします。tidyverseはデータの整形、グラフ化、文字列操作など統計学、データサイエンスの分野で必要な操作を使いやすくまとめたパッケージです。
library(tidyverse)
次にtibble表示の桁数を設定しておきます。tibble型はdata.frame型の拡張版でtidyverseパッケージでよく使われるデータ型です。デフォルトでは表示の桁数が有効数字3桁なので、ここでは5桁に修正しておきます。
# tibble表示の桁数設定
tibble_opt <- list(
"pillar.sigfig" = 5
)
options(tibble_opt)
getOption("pillar.sigfig")
次に乱数を発生させます。ここでは、平均(mean) = 0, 分散（\(sd^2 = 1^2\)）の正規分布に従う乱数を3×10万個発生させて、そのデータを3列×10万行のdara.frameに代入しています。
n <- 3 N <- 100000 # 乱数の発生 set.seed(0) x <- rnorm(n * N, mean = 0, sd = 1 ) %>% matrix(ncol = n) %>% as.data.frame()
標本分散、不偏分散のいずれを計算するにしても、偏差平方和を計算する必要がありますので、関数fun()として準備します。
\(fun(x) = \Sigma_{i = 1}^n(x_i – \bar{x})^2\)
# 偏差平方和計算用の関数（xx : データ）
fun <- function(xx){
sum((xx - mean(xx))^2)
}
xの各行について、偏差平方和yを求め、標本分散（y1）と不偏分散（y2）を求めます。
# xの各行について分散を計算 y <- apply(x, MARGIN = 1, FUN = fun) y1 <- y / n y2 <- y / (n - 1)
標本分散（y1）と不偏分散（y2）はサンプリングしたxの値に応じてばらつきます。どのようにばらついているのか確認するため、ヒストグラムを書いてみます。ヒストグラムを書くためにxとy1,y2を一つにまとめて、整形しておきます。
# グラフ化しやすいように整形
x$y1 <- y1
x$y2 <- y2
x <- x %>%
pivot_longer(cols =c(y1, y2),
names_to = "method",
values_to = "y") %>%
mutate(method = ifelse(method == "y1", "標本分散", "不偏分散"))
整形後のデータはこんな感じ。V1〜V3が乱数で、methodの方法で計算した結果がyにまとめられています。
# 計算された分散をグラフ化
ggplot(data = x) +
geom_histogram(aes(x = y, fill = method),
alpha = 0.5,
position = "identity",
binwidth = 0.2, boundary = 0) +
labs(x = "標本分散 or 不偏分散", y = "個数", fill = "計算方法") +
scale_x_continuous(breaks = seq(0, 6, by = 1)) +
coord_cartesian(xlim = c(0, 6)) +
theme_bw(base_family = "HiraKakuProN-W3") +
scale_fill_manual(values = c("blue", "red"))
xの真の分散は1であるのに（標準偏差は1であるのに）、標本分散は0〜4くらいまでばらついています。不偏分散は標本分散よりやや高く0〜6くらいまでばらついています。真の分散は1であるのに、かなり広くばらつきますね。
10万個の標本分散、不偏分散のデータの平均値を求めてみます。標本分散と不偏分散の両方がyに入っているので、一旦、method（=標本分散 or 不偏分散）でグループ分けをしてから、平均の計算をしています。
# 計算された分散（N個）の平均
x %>%
group_by(method) %>%
summarise(mean = mean(y))
# A tibble: 2 × 2 method mea <chr> <dbl> 1 標本分散 0.66509 2 不偏分散 0.99764
手元にあるデータから、母集団の分散（標準偏差）を推定するにあたり、推定である以上、真の値には一致はしないのですが、不偏分散の方が近しい値を算出する可能性が高いのです。不偏分散を使用するのはここに理由があります。
なお、標本分散、不偏分散とも確率的にばらつくので、データによっては標本分散の方が、真の値に近い場合もあり得ます。あくまでも、この推定方法で何回も行った場合に、不偏分散の方が平均的に真の値に近い値が得られるということです。
n = 100の場合
ソースコードはn <- 3となっていた部分のみ書き換えるだけなので、省略します。
n <- 100 # n <- 3の部分を書き換える
グラフは以下のようになります。
また、100個のデータから標本分散、不偏分散を計算する作業を10万回繰り返した場合のそれぞれの値の平均値は以下のようになります。
# A tibble: 2 × 2 method mean <chr> <dbl> 1 標本分散 0.99024 2 不偏分散 1.0002
不偏分散は真の値である1に近く、標本分散は真の値の\(\frac{n-1}{n} = \frac{99}{100}\)倍に近い値になっています。
\(\chi^2\)分布との関係
標本分散にしても、不偏分散にしても分母が違うだけで分子は同じ値（偏差平方和）です。
\(x\sim N(0, 1^2)\)
\(\bar{x}=\frac{1}{n}\Sigma_{i=1}^n{x_i}\)
分子（偏差平方和）： \(Z = \Sigma_{i = 1}^n(x_i-\bar{x})^2\)
偏差平方和（Z）の分布について考えてみます。
ヒストグラムだと高さがNによって変わってきてしまうので、縦軸を密度で表します（ヒストグラムの棒の総面積が1になるように無次元化します）。n = 3,100の時のヒストグラム（密度）と自由度2, 99の\(\chi^2\)分布の確率密度関数を重ねて描くと以下のようになります。
分散の分子である偏差平方和と\(\chi^2\)分布が密接に関係があるということです。また、\(\chi^2\)分布の自由度が和をとった個数nではなく、n-1であるところも興味深いところです。不偏分散の分母の値と偏差平方和が従う分布の自由度が同じn-1であることは、決して偶然ではなく、根本は同じ理由で説明されます。このため、不偏分散の分母がn-1である理由を、自由度が一つ減るからと説明されたりもします。
\(\chi^2\)分布については、こちらでも記事を書いていますのでぜひ参照ください。
まとめ
標本分散、不偏分散の意味をRによるシミュレーションで説明してみました。また、分散計算で出てくる偏差平方和と\(\chi^2\)分布の関係についてもシミュレーションで試してみました。
今回、Rを用いてシミュレーションを行いました。Rは統計計算に特化していると言われることも多いプログラミング言語ですが、グラフを作るのも得意ですし、今回のようにモンテカルロシミュレーションを行ったりもできます。また、公式サイトに登録されているパッケージを利用すれば、機械学習、ディープラーニング、微分方程式などの数値計算を行ったり、Webスクレイピングなどもできます。Rはデータに関わるさまざまな問題に対応できるプログラミング言語です。しかも無料（笑）。これを機にRに触れてみてはいかがでしょうか？
