# ーーーーーーー【分散推定値の不偏性】ーーーーーーー
#
# 　偏向推定値（var.biased） = 平方和／データ数
# 　不偏推定値（var.unbiased） = 平方和／自由度
# 
# 　標準正規乱数10個の無作為抽出を千回繰り返して，
# 　二つの分散推定値の平均を真値（=1）と比較する．
# ーーーーーーーーーーーーーーーーーーーーーーーーー

# 長さ1000の数値ベクトルをオブジェクトとして定義

   var.biased <- numeric(1000);
   var.unbiased <- numeric(1000);

# 標準正規乱数10個の抽出を千回ループ

   for (i in 1:1000)  {
     data <- rnorm(10, mean=0, sd=1);
     var.biased[i] <- sum((data - mean(data))^2)/10;
     var.unbiased[i] <- sum((data - mean(data))^2)/9
     }

# それぞれの分散推定値の平均値を計算

   cat("平均偏向推定値 =", mean(var.biased), "\n");

## 平均偏向推定値 = 0.8881874

   cat("平均不偏推定値 =", mean(var.unbiased), "\n");

## 平均不偏推定値 = 0.9868748

# ヒストグラム描画と平均推定値の表示
# 緑=真値=1
# 赤=平均偏向推定値
# 青=平均不偏推定値

   hist(var.biased, freq=F, density=15, angle=135);
   hist(var.unbiased, freq=F, density=10, angle=45, add=T);
   abline(v=1, col="green", lwd=5);
   abline(v=mean(var.biased), col="red", lwd=3); 
   abline(v=mean(var.unbiased), col="blue", lwd=3)