母分散の信頼区間推定のシミュレーション

第109回Tokyo.Rの書籍抽選で、『数値シミュレーションで読み解く統計のしくみ～Rでためしてわかる心理統計』(小杉考司・紀ノ定保礼・清水裕士著, 技術評論社)をいただきました。統計的推測や分析をシミュレーションで理解するという本で、たいへん参考になります。

この本では直接はあつかわれていない題材で、母分散の信頼区間の推定をやってみます。

問題

正規分布\(N(\mu, \sigma^2)\)から抽出された標本\(x_1, x_2, x_3, \dots, x_n\)から母分散の95%信頼区間を推定します。

シミュレーション

理論的には、標本の偏差平方和を\(T^2 = \sum_{i=1}^{n}(x_n - \bar{x})^2\)とすると、\(\frac{T^2}{\sigma^2}\)は自由度\(n-1\)のカイ二乗分布にしたがいます。これを確かめてみます。

設定

その前に、ライブラリの読み込みと、偏差平方和の関数の定義です。

library(ggplot2)

sim_sumsq <- function(n = 10, mu = 0, sigma = 1) {
  x <- rnorm(n, mu, sigma)
  x_bar <- mean(x)
  sum((x - x_bar)^2)
}

サンプルサイズn = 10, 平均mu = 0, 標準偏差sigma = 3として、シミュレートしてみます。

n <- 10
mu <- 0
sigma <- 3
R <- 2000

set.seed(123)
T2 <- replicate(R, sim_sumsq(n, mu, sigma))
X2 <- T2 / sigma^2

結果

シミュレーションで生成した、偏差平方和を分散で割った値をヒストグラムに、理論的なカイ二乗分布の確率密度関数の値を曲線でプロットします。

ggplot(data.frame(X = X2)) +
  geom_histogram(aes(x = X, y = after_stat(density)),
                 binwidth = 1) +
  geom_function(fun = ~ dchisq(.x, df = n - 1), colour = "red") 

プロットしてみると、両者はうまくかさなっていることがわかりました。

母分散の信頼区間

したがって、自由度\(n-1\)のカイ二乗分布の2.5%点を\(\chi^2_{0.025}\)、97.5%点を\(\chi^2_{0.975}\)とすると、以下のようになります。

\[ \chi^2_{0.025} < \frac{T^2}{\sigma^2} < \chi^2_{0.975} \]

これを変形すると、母分散の95%信頼区間は以下のようになります。

\[ \frac{T^2}{\chi^2_{0.975}} < \sigma^2 < \frac{T^2}{\chi^2_{0.025}} \]

シミュレーションで計算された値

シミュレーションで計算された偏差平方和を\(\chi^2_{0.975}\)および\(\chi^2_{0.025}\)で割った値はそれぞれ以下のようになりました。

summary(T2 / qchisq(0.975, df = n - 1))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5213  2.8342  3.9727  4.2649  5.4681 14.7835
summary(T2 / qchisq(0.025, df = n - 1))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.672  19.966  27.985  30.044  38.519 104.142