なんでn-1なんですか

不偏分散を計算する際になぜnではなくn-1で割るのか。理由はともかくnで割ってマズイのは割ってみれば分かります。

やってみよう

• サイコロ: 1から6までが等しい確率で出る理想的なサイコロ様を用意
• 分散の計算: サイコロをn回振るごとに分散と不偏分散を計算
  • nはとりあえず2から50まで
  • それぞれのnについて1000回分散を計算し、その平均をそのn数における分散とする

## 分散(nで割るやつ)
var <- function(data) sum((mean(data) - data)^2)/length(data)
## 不偏分散(n-1で割るやつ) 
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なんですか。