不偏分散を計算する際になぜnではなくn-1で割るのか。理由はともかくnで割ってマズイのは割ってみれば分かります。
やってみよう
- サイコロ: 1から6までが等しい確率で出る理想的なサイコロ様を用意
- 分散の計算: サイコロをn回振るごとに分散と不偏分散を計算
- nはとりあえず2から50まで
- それぞれのnについて1000回分散を計算し、その平均をそのn数における分散とする
var <- function(data) sum((mean(data) - data)^2)/length(data)
var2 <- function(data) sum((mean(data) - data)^2)/(length(data)-1)
dice <- 1:6
true.var <- var(dice)
n.limit <- 50
rep.limit <- 1000
mean.var <- mean.var2 <- numeric(n.limit-1)
for(n in 2:n.limit){
v <- v2 <- numeric(rep.limit)
for(i in 1:rep.limit){
v[i] <- var(sample(dice, n, rep=TRUE))
v2[i] <- var2(sample(dice, n, rep=TRUE))
}
mean.var[n-1] <- mean(v)
mean.var2[n-1] <- mean(v2)
}
plot(2:n.limit, ylim=c(0, 4), type="n",
xlab="サイコロを振った回数", ylab="分散の平均値")
points(2:n.limit, mean.var, col=2)
points(2:n.limit, mean.var2, col=3)
abline(h = true.var, col="skyblue")
legend("bottomright", inset=0.05,
legend=c("分散(nで割ったやつ)の平均",
"不偏分散(n-1で割ったやつ)の平均",
"真の分散"),
pch=c(1, 1, NA),
lty=c(NA, NA, 1),
col=c(2, 3, "skyblue"),
box.lty=0)
結果
不偏分散マジ不偏。それで結局なんでn-1なんですか。