母集団(成人男性全体)では、平均身長が170cm、身長の標準偏差が6cm であり、身長は正規分布に従うことを知っているとする。つまり、母分散は \(6^2 = 36\) である。
母集団の身長分布は以下のようになる。
以下では、様々な標本サイズで標本を抽出し、その標本の分散を求める作業を 繰り返すことで、分散に対する理解を深めよう。
身長の観測値を \(x\) で現すことにしよう。統計学1で(おそらく)学んだとおり、\(x\) の分散 \(s^2\) は、以下の式で定義される。 \[s^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2\] ただし、\(n\) は標本サイズである。
\(n=10\) として母集団から標本を1つ抽出し、この分散を計算してみよう。
抽出した標本について、1つひとつの値 (\(x_1, x_2, \dots, x_{10}\)) を確認する。
## [1] 175.7723 163.8309 173.8970 176.5091 161.1755 172.1482 169.9440
## [8] 164.7393 165.6629 164.9124
次に、偏差 (deviation) \(x - \bar{x}\) を計算する。
## [1] 6.913121 -5.028296 5.037844 7.649967 -7.683695 3.289069 1.084874
## [8] -4.119859 -3.196253 -3.946773
偏差の二乗 (squares) を計算する。
## [1] 47.791240 25.283762 25.379874 58.521988 59.039169 10.817977 1.176952
## [8] 16.973235 10.216032 15.577014
偏差の二乗を合計する(これを偏差平方和と呼ぶ)。
## [1] 270.7772
これを標本サイズ \(n\) で割ったものが、標本の分散 \(s^2\) である。
## [1] 27.07772
毎回この作業をするのは面倒なので、標本分散を計算する関数を作ってしまおう。
get_var <- function(x) {
## 標本の分散を計算する関数
## x は1つの標本(値が複数入ったベクトル)
n <- length(x) # 標本サイズを求める
dev <- x - mean(x) # 偏差
ssd <- sum(dev^2) # 偏差平方和
return(ssd / n) # 標本の分散を返す
}
この関数を使って、標本の分散を求めてみよう。
## [1] 27.07772
先ほどと同じ値が得られた。
しかし、Rには分散を求める var()
という関数がある。この関数を使って分散を求めてみよう。
## [1] 30.08636
答えが一致しない。なぜだろうか。
実は、Rのvar()
が計算しているのは、不偏分散と呼ばれる分散である。この分散は、母分散の不偏推定量である。不偏分散 \(u^2\) は、次のように定義される。 \[u^2 = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})^2.\] つまり、分散の式の分母を \(n\) から \(n-1\) に変えたものが不偏分散である。 分母を小さくしているので、分散自体は大きく推定される。このように、あえて分散が大きくなるようにするのは、標本分散 \(s^2\) が母分散を過小推定してしまうためである。これをシミュレーションで確かめてみよう。
ここでは、上で定義した母集団から \(n=10\) の標本を1万回抽出し、それぞれの標本について、標本分散 \(s^2\) と不偏分散 \(u^2\) を計算し、分散の分布を確認してみよう。
まず、シミュレーションの設定を行う。
次に、結果を保存するための行列を用意する。保存する対象は2つ(2種類の分散)で、シミュレーションの繰り返し回数が n_sim なので、n_sims行2列の行列を用意する。
1列目に標本分散を、2列目に不偏分散を保存する。
forループを使ってシミュレーションを実行する。
for (i in 1:n_sims) {
x <- sample(pop, size = n, replace = FALSE)
res_sim_1[i, 1] <- get_var(x)
res_sim_1[i, 2] <- var(x)
}
結果を確認する。まず、標本分散の分布を確認してみよう。
hist(res_sim_1[, 1], col = "gray",
xlab = "標本分散", ylab = "度数", main = "標本分散の標本分布")
abline(v = get_var(pop), col = "red")
標本分散の平均値を確認する。
## [1] 32.25378
母分散(母集団の分散)は、
## [1] 36.01808
だから、標本分散の平均値は母分散より小さい。つまり、平均すると、標本分散は母分散を過小推定する。言い換えると、標本分散は偏った推定量である。
同様に、不偏分散を確認してみよう。
hist(res_sim_1[, 2], col = "gray",
xlab = "不偏分散", ylab = "度数", main = "不偏分散の標本分布")
abline(v = get_var(pop), col = "red")
不偏分散の平均値を確認する。
## [1] 35.83753
母分散(母集団の分散)は、
## [1] 36.01808
だから、不偏分散の平均値は、母分散と(ほぼ)同じである。つまり、平均すると、不偏分散は母分散を正しく推定することができる。言い換えると、不偏分散は母分散の不偏推定量である。
カイ二乗 (\(\chi^2\)) 分布は、自由度 \(df \geq 1\) によってその形を変える。
たとえば、自由度1, 3, 5, 10のカイ二乗分布は、以下のように分布する。
x <- seq(0, 20, length = 1000)
chi1 <- dchisq(x, df = 1)
chi3 <- dchisq(x, df = 3)
chi5 <- dchisq(x, df = 5)
chi10 <- dchisq(x, df = 10)
plot(x, chi1, type = "l", lwd = 2, xlim = c(0, 20), ylim = c(0, 0.6),
ylab = "確率密度", main = expression(paste(chi^2, "分布")))
abline(h = 0)
lines(x, chi3, lwd = 2, col = "red")
lines(x, chi5, lwd = 2, col = "skyblue")
lines(x, chi10, lwd = 2, col = "orange")
legend("topright", lty = 1, lwd = 2, col = c("black","red","skyblue","orange"),
legend = c("1", "3", "5", "10"), title = "自由度")
カイ二乗分布の形状を確認するための関数を用意したので、これを使って色々なカイ二乗分布の形状を確認してみよう(この関数の中身を理解する必要はない)。
plot_chisq <- function(df) {
upper <- max(c(df + 10, df*2))
x <- seq(from = 0, to = upper, length = 1000)
chisq <- dchisq(x, df = df)
plot(x, chisq, type = "l", col = "royalblue", lwd = 1.5,
ylab = "確率密度",
main = paste("自由度", df, "のカイ二乗分布"))
}
自由度1のカイ二乗分布は、
自由度3のカイ二乗分布は、
実習課題:自由度 (df) の値をによって、カイ二乗分布がどのように変化するか確かめてみよう。
上で定義した母集団から、標本サイズ \(n=10\) の標本を1つ取り出す。
この標本から母集団の分散を推定したい場合、点推定値は不偏分散 \(u^2\) なので、
## [1] 51.23188
である。
母分散の区間推定には、自由度 \(n-1\) のカイ二乗分布を利用する。 母分散の信頼区間は、 \[\frac{n-1}{\chi_U^2} u^2 \leq \sigma^2 \leq \frac{n-1}{\chi_L^2} u^2\] と表される。\(\chi_L^2\) と\(\chi_U^2\)には、自由度と信頼区間の信頼度(パーセンテージ)に応じて異なる数字を入れる。
例として、95パーセント信頼区間を求めることにしよう。私たちの標本の大きさは \(n=10\)なので、自由度は \(n-1=9\) である。カイ二乗分布の中心部分95%を推定に利用したい。言い換えると、分布の端にある5%分は除外したい。両端から2.5%ずつ除外することにすると、下側の点である\(\chi_L^2\)は、
## [1] 2.700389
であり、上側の点である\(\chi_U^2\)は、
## [1] 19.02277
であることがわかる。
これらの値を使うと、母分散 \(\sigma^2\) の95%信頼区間の下限値は、
## [1] 24.23869
となる。同様に、上限値は、
## [1] 170.7483
となる。したがって、今回の標本から得られる母分散の95%信頼区間は、[24.24, 170.75] である。
毎回この作業を繰り返すのは面倒なので、母分散の信頼区間を計算する関数を作ってしまおう。 信頼度は選べるようにする。
confint4var <- function(x, level = 0.95) {
## 母分散の信頼区間を計算する関数
## 引数: x = 標本(観測値のベクトル)
## level = 信頼度(既定値は0.95)
## 返り値:不偏分散、信頼区間の下限値、信頼区間の上限値のベクトル
n <- length(x)
u_sq <- var(x)
chisq_l <- qchisq(p = (1 - level) / 2, df = n - 1, lower.tail = TRUE)
chisq_u <- qchisq(p = (1 - level) / 2, df = n - 1, lower.tail = FALSE)
lb <- (n - 1) * u_sq / chisq_u
ub <- (n - 1) * u_sq / chisq_l
return(c(var = u_sq, lb = lb, ub = ub))
}
試しにこの関数を使ってみよう。
## var lb ub
## 51.23188 24.23869 170.74829
上で計算したのと同じ結果が得られた。
この標本から得られる母分散の87%信頼区間を求めてみよう。
## var lb ub
## 51.23188 28.65448 127.62425
求める区間は、[28.65, 127.62] である。
日本の成人女性の平均身長が、平均 160cm, 標準偏差 6cm の正規分布にしたがっていると仮定する。 このとき、標本サイズ \(n = 100\) の標本を3つ抽出し、それぞれの標本について以下の推定値を求めなさい。