LoginSignup
5
3
@tabintone(Komatsu)

不偏分散の平方根は母標準偏差の不偏推定量ではないという話

Last updated at Posted at 2020-11-21

はじめに

製造業に転職してから「バラツキ」に関する話を聞くことが多くなりました。実際にデータから標準偏差$\sigma$を計算して、平均値$\pm3\sigma$などの数値を品質管理に利用することは多いと思います。そんな中で「不偏分散の平方根は母標準偏差の不偏推定量ではない」ということを知ったので自分なりにまとめておくことにします。

不偏推定量

まず、不偏推定量とは何かを簡単に説明します。得られたデータから母数$\theta$を推定するための推定量$\hat{\theta}$を作ったとき、$\hat{\theta}$が$\theta$の不偏推定量であるとは

 E[\hat{\theta}]=\theta

が成立することを言います。$E[\hat{\theta}]$は$\hat{\theta}$の期待値です。この式は$\hat{\theta}$の値が母数$\theta$の周りに散らばることを意味しています。したがって、不偏推定量を用いることで平均的には正しく母数$\theta$を推定できることが保証されます。不偏性は推定量が満たしていてくれると嬉しい性質の一つです。

標準偏差の求め方

$n$個のデータ$ y_1, y_2,\cdots,y_n$があるとき、不偏分散は以下のように計算します。

 \hat{\sigma}^2  = \frac{1}{n-1}\sum_{i=1}^{n}(y_i-\overline{y})^2

$\overline{y}$は$n$個のデータの平均値です。$\hat{\sigma}^2$の期待値を計算すると母分散$\sigma^2$に等しくなりますので、$\hat{\sigma}^2$は母分散$\sigma^2$の不偏推定量になっています。

標準偏差については以下のように不偏分散の平方根で計算することが多いです。

 \hat{\sigma}  = \sqrt{\hat{\sigma}^2}

それでは、$\hat{\sigma}$は母標準偏差$\sigma$の不偏推定量になっているでしょうか？

シミュレーションで確かめてみる

とりあえず、$\hat{\sigma}$が不偏推定量になっているかシミュレーションして確かめてみることにします。正規分布$N(10,5^2)$から乱数を10個発生させて先ほどの計算式で標準偏差$\hat{\sigma}$を計算する、という作業をたくさん繰り返します。$\hat{\sigma}$が不偏推定量であれば、たくさん計算した$\hat{\sigma}$の平均値は母標準偏差である5に近い値になるはずです。

Rを使ってシミュレーションしてみます。比較のため平均と分散についても計算することにします。まずは1000回繰り返してみましょう。

R
#　シミュレーションのための関数の作成
calc_error <- function(n, repeat_num, mu, sigma){
  # 乱数発生
  test_data <- tibble(value = rnorm(n*repeat_num, mean = mu, sd = sigma)) %>% 
    dplyr::mutate(id = as.factor(row_number()%%repeat_num))
  
  # 平均、不偏分散、標準偏差を算出
  estimates <- test_data %>% 
    dplyr::group_by(id) %>% 
    dplyr::summarise(mu = mean(value), 
                     sigma_sq = var(value), 
                     sigma = sd(value))
  
  # 計算した平均、不偏分散、標準偏差の平均を算出
  mu_es <- mean(estimates$mu)
  sigma_sq_es <- mean(estimates$sigma_sq)
  sigma_es <- mean(estimates$sigma)
  
  # 真の値と推定値のパーセント誤差を計算
  mu_error <- ((mu_es-mu)/mu)*100
  sigma_sq_error <- ((sigma_sq_es-sigma^2)/(sigma^2))*100
  sigma_error <- ((sigma_es-sigma)/sigma)*100
  
  # 結果をリストにして出力
  result <- list(mu = mu_es, sigma_sq = sigma_sq_es, sigma = sigma_es, 
                 mu_error = mu_error, sigma_sq_error = sigma_sq_error, 
                 sigma_error = sigma_error)
  return(result)
}
R
# 1000回でシミュレーションの実行
calc_error(n = 10, repeat_num = 1000, mu = 10, sigma = 5)

1000回ではシミュレーションの結果は以下のようになりました。

<1000回> 平均 分散 標準偏差
真の値 　 10 25 5
推定値の平均 10.05704 24.77481 4.841726
パーセント誤差 0.57% -0.90% -3.16%

標準偏差の値のズレが一番大きいですが、平均と分散もまだズレているのでさらに繰り返し回数を1万回、10万回に増やしてみることにします。

R
calc_error(n = 10, repeat_num = 10000, mu = 10, sigma = 5)   # 1万回
calc_error(n = 10, repeat_num = 100000, mu = 10, sigma = 5)  # 10万回
<1万回> 平均 分散 標準偏差
真の値 　 10 25 5
推定値の平均 10.00483 25.04079 4.8668
パーセント誤差 0.05% 0.16% -2.66%
<10万回> 平均 分散 標準偏差
真の値 　 10 25 5
推定値の平均 10.0027 24.98578 4.863767
パーセント誤差 0.03% -0.06% -2.72%

繰り返し数を多くすると平均と分散の推定値の平均は母平均・母分散に近い値になっています。しかし、標準偏差に関しては繰り返し数を増やしても$2.7$%ほど小さい値になってしまいます。10万回も繰り返したのに真の値に近づかなかったので、$\hat{\sigma}$は不偏推定量ではない気がしてきました。

期待値を計算してみる1

シミュレーションの結果から$\hat{\sigma}$は不偏推定量ではない気がしてきたので、実際に$\hat{\sigma}$の期待値を求めてみます。$n$個のデータ$ y_1, y_2,\cdots,y_n$がすべて独立に正規分布$N(\mu, \sigma^2)$に従うとします。このとき、

 \hat{\sigma}^2  = \frac{1}{n-1}\sum_{i=1}^{n}(y_i-\overline{y})^2

 \hat{\sigma}  = \sqrt{\hat{\sigma}^2}

として、$E[\hat{\sigma}]$の値を求めることにします。まず、母分散を$\sigma^2$とすると、

 X  = \frac{(n-1)\hat{\sigma}^2}{\sigma^2}

は自由度$n-1$の$\chi^2$分布に従います。$X$を用いて$\hat{\sigma}$は以下のように表すことができます。

 \hat{\sigma}  = \sqrt{\hat{\sigma}^2}=\sqrt{\frac{\sigma^2}{n-1}X}

したがって、$\hat{\sigma}$の期待値は

 E[\hat{\sigma}]=E\Biggl[\sqrt{\frac{\sigma^2}{n-1}X}\Biggl]=\sqrt{\frac{\sigma^2}{n-1}}\cdot E[\sqrt{X}]

と表すことが出来ます。

$X$は自由度$n-1$の$\chi^2$分布に従う確率変数なので、期待値を計算するために$\chi^2$分布の確率密度関数ついておさらいしておきます。2以下の確率密度関数を持つ確率分布をガンマ分布と呼び、$Ga(\alpha, \lambda)$と表します。

f(x) = \frac{\lambda^\alpha}{\Gamma(\alpha)} x^{\alpha-1}e^{-\lambda x}\hspace{4pt}(x\ge 0), \hspace{8pt}0\hspace{4pt}(x\lt 0)

ガンマ分布には$\alpha$と$\lambda$の2つのパラメータがありますが、$\alpha=\frac{n}{2}$、$\lambda=\frac{1}{2}$としたガンマ分布$Ga(\frac{n}{2}, \frac{1}{2})$を特別に自由度$n$の$\chi^2$分布と呼びます。なお、$\Gamma$はガンマ関数という関数ですが、期待値を求めるにあたって詳細を知っている必要はありません。また、実際に計算してガンマ関数の値を求めるときはRの関数にお任せするので、そんな関数があるんだなくらいに思っておいてください。

$\chi^2$分布の確率密度関数がわかったので期待値を計算していきます。$X$は自由度$n-1$の$\chi^2$分布に従うので、確率密度関数は$\alpha=\frac{n-1}{2}$、$\lambda=\frac{1}{2}$として、

 E[\hat{\sigma}]=\sqrt{\frac{\sigma^2}{n-1}}\cdot E[\sqrt{X}]

=\sqrt{\frac{\sigma^2}{n-1}}\int_{-\infty}^{\infty}\sqrt{x}\hspace{2pt}f(x)dx

 =\sqrt{\frac{\sigma^2}{n-1}}\int_{0}^{\infty}\sqrt{x}\cdot\frac{\bigl(\frac{1}{2}\bigl)^\frac{n-1}{2}}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot x^{(\frac{n-1}{2}-1)}\cdot e^{-\frac{1}{2}x} dx

$\sqrt{x}=x^\frac{1}{2}$を利用して式を変形していきます。

 =\sqrt{\frac{\sigma^2}{n-1}}\int_{0}^{\infty}\frac{\bigl(\frac{1}{2}\bigl)^\frac{n-1}{2}}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot x^{(\frac{n}{2}-1)}\cdot e^{-\frac{1}{2}x} dx

 =\sqrt{\frac{\sigma^2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n}{2}\bigl)}\cdot\frac{\bigl(\frac{1}{2}\bigl)^\frac{n}{2}}{\bigl(\frac{1}{2}\bigl)^\frac{n}{2}}\int_{0}^{\infty}\frac{\bigl(\frac{1}{2}\bigl)^\frac{n-1}{2}}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot x^{(\frac{n}{2}-1)}\cdot e^{-\frac{1}{2}x} dx

積分の中身の分数部分を入れ替えます。

 =\sqrt{\frac{\sigma^2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot\frac{\bigl(\frac{1}{2}\bigl)^\frac{n-1}{2}}{\bigl(\frac{1}{2}\bigl)^\frac{n}{2}}\int_{0}^{\infty}\frac{\bigl(\frac{1}{2}\bigl)^\frac{n}{2}}{\Gamma\bigl(\frac{n}{2}\bigl)}\cdot x^{(\frac{n}{2}-1)}\cdot e^{-\frac{1}{2}x} dx

すると、積分の中身は自由度$n$の$\chi^2$分布の確率密度関数になります。確率密度関数を定義域全体で積分すると1になるので、積分を消去することができて、

 =\sqrt{\frac{\sigma^2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot\frac{\bigl(\frac{1}{2}\bigl)^\frac{n-1}{2}}{\bigl(\frac{1}{2}\bigl)^\frac{n}{2}}

 =\sqrt{\frac{\sigma^2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n-1}{2}\bigl)}\cdot \sqrt{2}

$\sigma^2$をルートの外に出して式を整理すると、

 =\sigma\cdot\sqrt{\frac{2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n-1}{2}\bigl)}

と表すことが出来ます。上の式は母標準偏差$\sigma$とデータ数$n$に依存する定数の積の形になっているので、

 C_n=\sqrt{\frac{2}{n-1}}\cdot\frac{\Gamma\bigl(\frac{n}{2}\bigl)}{\Gamma\bigl(\frac{n-1}{2}\bigl)}

と置くと、$\hat{\sigma}$の期待値は

 E[\hat{\sigma}]=C_n\sigma

となります。

$E[\hat{\sigma}]=\sigma$とはならなかったので、$\hat{\sigma}$は母標準偏差$\sigma$の不偏推定量ではありませんでした。シミュレーションから示唆された通りですね。

ちなみに、$\hat{\sigma}$の期待値は母標準偏差$\sigma$の$C_n$倍なので、不偏推定量にするには$C_n$で割ってズレを補正する必要があります。したがって、母標準偏差の不偏推定量は

 \frac{\hat{\sigma}}{C_n}

になります。また、$C_n$を標準偏差のcorrection factorと呼びます。3

シミュレーションの結果を確認してみる

10万回のシミュレーションで得た標準偏差の推定値の平均は$4.863767$でした。この値は$\hat{\sigma}$の期待値である、$C_n\sigma$と一致しているでしょうか？念の為、確認しておくことにします。

R
# Cnを計算するための関数の作成
Cn <- function(n){
  sqrt(2/(n-1))*(gamma(n/2)/gamma((n-1)/2))
}

シミュレーションでは10個の乱数を発生させて標準偏差を計算したので$n=10$です。$C_{10}$の値を以下のように計算します。

R
# C10の値を計算する
> Cn(10)
[1] 0.9726593

母標準偏差は5としてシミュレーションを行ったので、$\hat{\sigma}$の期待値は、

 E[\hat{\sigma}]=C_{10}\sigma=0.9726593\hspace{1pt}\times\hspace{1pt}5=4.863296

となります。シミュレーションで得た$4.863767$とほとんど同じ値になっています。これでシミュレーションの結果を理論的に説明することができました。

Correction factor Cn

correction factor $C_n$が$1$から離れた値になるほど$\hat{\sigma}$の期待値と母標準偏差の乖離が大きくなり、$\hat{\sigma}$は偏りの大きな推定量になっていきます。$C_n$の値はデータ数に依存するため、$C_n$とデータ数の間にどのような関係があるかを確認しておくことにします。$x$軸をデータ数、$y$軸をcorrection factor $C_n$としてプロットすると以下のようになります。赤線は$C_n=1$の線です。

plot.png

$C_n$は常に$1$より小さいため、$\hat{\sigma}$の期待値は常に母標準偏差よりも小さくなることが分かります。また、データ数が多くなるほど$C_n$が$1$に近づいて行くのが確認できます。ある程度データ数が多いときは単純に不偏分散の平方根を計算することで、母標準偏差を不偏に推定できていると考えても特に問題なさそうです。しかし、データ数が少ないときは$C_n$は$1$から離れた値になり、$\hat{\sigma}$の期待値と母標準偏差のズレは大きくなります。そのような場合は、標準偏差を小さく推定してしまう可能性が高くなるのでズレを補正した$\hat{\sigma}/{C_n}$を母標準偏差の推定量とすることを検討しても良いかもしれません。4

ただし、$C_n$は「$n$個のデータ$ y_1, y_2,\cdots,y_n$がすべて独立に正規分布$N(\mu, \sigma^2)$に従う」という仮定の下で導出した値であることに注意する必要があります。5データが正規分布に従わない場合は$C_n$で補正しても不偏推定量にはなりません。6

おわりに

ここまで不偏分散の平方根が母標準偏差の不偏推定量にならないこと、データが正規分布に従う場合の母標準偏差の不偏推定量の計算方法をまとめました。今まであまり意識したことはありませんでしたが、単に標準偏差と言っても奥が深いですね。

  1. こちらのサイトを参考に自分が理解しやすい形に直して期待値を計算しています。https://bellcurve.jp/statistics/blog/13645.html ↩

  2. ガンマ分布と$\chi^2$分布の確率密度関数については、東京大学出版会の「基礎統計学I 統計学入門」を参考にしています。 ↩

  3. おそらく$C_n$の逆数をcorrection factorとして定義する方が一般的だと思いますが、自分の中での分かりやすさを優先してここでは$C_n$そのままをcorrection factorと呼ぶことにしています。$C_n$の逆数をcorrection factorとした場合は、correction factorを不偏分散の平方根に掛けることになります。 ↩

  4. correction factorによる補正を行う場合は文章中に記載したコードで計算しても良いですが、計算済みのcorrection factorの表を活用した方が楽かもしれません。検索するとたくさん出てきます。 ↩

  5. 一般に標準偏差の不偏推定量は分布の仮定なしで求めることは出来ないらしいです。これも脚注1のサイトに書いてありました。 ↩

  6. 一般的に小さめに推定することになるので、補正を行うことで不偏推定量にはならなくても多少ましな推定量にはなるかもしれません。また、こちらのサイトではデータが指数分布に従うときにcorrection factorで補正しても不偏推定量にならないことをシミュレーションで示しています。https://biolab.sakura.ne.jp/unbiased-standard-deviation.html ↩

5
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
3