真値と標準偏差の推定

(886d) 更新

はじめての誤差論

平均値および標準偏差の期待値 †

x の確率分布を p(x) とし、その期待値を x^* 、分散を (\sigma_x^*)^2 とする。

n 解の測定で、あるデータセット x_1,x_2,\dots,x_n が得られる確率は、

　 \prod_{k=1}^n p(x_k)

である。

このデータセットに対して平均値は \bar x=\frac{1}{n}\sum_{k=1}^n x_k だから、 「平均値の期待値」を求めるにはこれに確率を掛けて積分すればよく、

　 \langle\bar x\rangle &=\int dx_1\int dx_2\dots\int dx_n \Bigg(\frac{1}{n}\sum_{k=1}^n x_k \Bigg)\Bigg(\prod_{m=1}^n p(x_m)\Bigg)\\ &=\frac{1}{n}\sum_{k=1}^n\underbrace{\int p(x_1) dx_1}_{=\,1}\underbrace{\int p(x_2) dx_2}_{=\,1}\dots\underbrace{\int x_kp(x_k)dx_k}_{=\,x^*}\dots\underbrace{\int p(x_n) dx_n}_{=\,1} \\ &=\frac{1}{n}\sum_{k=1}^n x^*\\ &=x^*

したがって、平均値の期待値は真値と一致する。

一方、分散の期待値は、

　 \langle\sigma_x^2\rangle &=\int dx_1\int dx_2\dots\int dx_n \Bigg(\frac{1}{n}\sum_{k=1}^n (x_k-\bar x)^2 \Bigg)\Bigg(\prod_{m=1}^n p(x_m)\Bigg)\\

二乗のかかっているところを展開すると、

　 \sum_{k=1}^n( x_k-\bar x )^2 &=\sum_{k=1}^n x_k^2-2\bar x\sum_{k=1}^nx_k +\sum_{k=1}^n \bar x^2\\ &=\sum_{k=1}^n x_k^2-2n\bar x^2 +n \bar x^2\\ &=\sum_{k=1}^n x_k^2-n\bar x^2\\ &=\sum_{k=1}^n x_k^2-\frac{1}{n}\sum_{k=1}^n\sum_{l=1}^n x_kx_l\\ &=\frac{n-1}{n}\sum_{k=1}^n x_k^2-\frac{1}{n}\sum_{k=1}^n\sum_{l\ne k} x_kx_l\\

したがって、

　 \langle\sigma_x^2\rangle &=\int dx_1\int dx_2\dots\int dx_n \Bigg(\frac{n-1}{n^2}\sum_{k=1}^n x_k^2-\frac{1}{n^2}\sum_{k=1}^n\sum_{l\ne k} x_kx_l \Bigg)\Bigg(\prod_{m=1}^n p(x_m)\Bigg)\\ &= \frac{n-1}{n^2}\sum_{k=1}^n (x^2)^*-\frac{1}{n^2}\sum_{k=1}^n\sum_{l\ne k} (x^*)^2 \\ &= \frac{n-1}{n} (x^2)^*-\frac{n-1}{n} (x^*)^2 \\ &= \frac{n-1}{n}(\sigma_x^*)^2

すなわち、データの分散の期待値は真の分散とは等しくならず、 その \frac{n-1}{n} と一致する。

したがって、データの標準偏差から真の標準偏差を推定するには

　 \sigma_x^*=\sqrt\frac{n}{n-1}\sigma_x

とすべきなのである。

具体例（正規分布の時） †

特に、データが正規分布に従うと仮定して、平均値と標準偏差の推定を行うことで上記結果に対する理解を深めよう。

データの確率分布が

　 p(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^2/2\sigma^2}

であれば、あるデータセット x_1,x_2,\dots,x_n が得られる確率は、

　 \prod_{k=1}^n p(x_k)=\bigg(\frac{1}{\sqrt{2\pi}\sigma}\bigg)^n\exp\bigg[-\frac{1}{2\sigma^2}\sum_{k=1}^n(x_k-\mu)\bigg]

である。

この確率を最大化するように \mu,\sigma を選べば、それらは x^*,\sigma_x^* の最良推定値となるであろう。

[ ] 内は、

　 \sum_{k=1}^n(x_k-\mu)^2 &=\sum_{k=1}^nx_k^2-2\mu \sum_{k=1}^nx_k+\sum_{k=1}^n\mu^2\\ &=\sum_{k=1}^n(x_k^2-\bar x^2)+n\bar x^2-2n\mu\bar x+n\mu^2\\ &=n\sigma_x^2+n(\bar x-\mu)^2\\

となるから、上記の確率は平均値と分散で書き表せて、

　 P(\mu,\sigma) &=\bigg(\frac{1}{\sqrt{2\pi}\sigma}\bigg)^n\exp\bigg[-\frac{n}{2\sigma^2}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg]\\ &=\bigg(\frac{1}{\sqrt{2\pi}\sigma}\exp\bigg[-\frac{1}{2\sigma^2}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg]\bigg)^n\\

n=1,\bar x=0,\sigma_x=1 のときこの関数をプロットしてみると次のようになる。 横軸が \mu 縦軸が \sigma であり、明るい部分ほど確率が高い。 赤線については後述する。

mu_sigma.jpg

式の形から明らかなとおり、この関数は \mu に対して偶関数であり、 任意の \sigma に対して \mu=\bar x で最大値を取る。

一方、 \sigma に対しては、

　 \frac{\partial}{\partial \sigma}P(\mu,\sigma) &=n\bigg(\frac{1}{\sqrt{2\pi}\sigma}\exp\bigg[-\frac{1}{2\sigma^2}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg]\bigg)^{n-1}\\ &\hspace{2cm}\frac{1}{\sqrt{2\pi}}\bigg(-\frac{1}{\sigma^2}+\frac{1}{\sigma^4}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg)\exp\bigg[-\frac{1}{2\sigma}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg]\\ &=0

と置けば、

　 -\frac{1}{\sigma^2}+\frac{1}{\sigma^4}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}=0

　 \sigma^2=\sigma_x^2+(\bar x-\mu)^2

となって、 \mu=\bar x,\sigma=\sigma_x で最大値を取ることが分かる。 (上のグラフの赤線が最大値を与える \sigma を表す)

それにも関わらず分散の最良推定値が \sigma_x ではなく \frac{n}{n-1}\sigma_x で与えられるのはどうしてかというと、 \mu\ne\bar x では確率が最大になるのが \sigma>\sigma_x のときであるためだ。

尤度分布を \mu 方向に積分して \sigma のみの関数とすると、

　 P(\sigma)&=\int d\mu\,P(\mu,\sigma)\\ &=\int d\mu\bigg(\frac{1}{\sqrt{2\pi}\sigma}\exp\bigg[-\frac{1}{2\sigma^2}\big\{\sigma_x^2+(\bar x-\mu)^2\big\}\bigg]\bigg)^n\\ &=\bigg(\frac{1}{\sqrt{2\pi}\sigma}\exp\bigg[-\frac{\sigma_x^2}{2\sigma^2}\bigg]\bigg)^n \int d\mu\exp\bigg[-\frac{n}{2\sigma^2}(\bar x-\mu)^2\bigg]\\ &=\bigg(\frac{1}{\sqrt{2\pi}\sigma}\exp\bigg[-\frac{\sigma_x^2}{2\sigma^2}\bigg]\bigg)^n \frac{\sqrt{2\pi}\sigma}{\sqrt n}\\ &=\frac{1}{\sqrt n}\bigg(\frac{1}{\sqrt{2\pi}\sigma}\bigg)^{n-1}\exp\bigg[-\frac{n\sigma_x^2}{2\sigma^2}\bigg]\\

この最大値は、

　 \frac{\partial}{\partial\sigma}P(\sigma) &=\frac{1}{\sqrt n}\bigg(\frac{1}{\sqrt{2\pi}\sigma}\bigg)^{n-1}\bigg(-\frac{n-1}{\sigma}+\frac{n\sigma_x^2}{\sigma^3}\bigg)\exp\bigg[-\frac{n\sigma_x^2}{2\sigma^2}\bigg]\\ &=0

と置いて、

　 \sigma^2=\frac{n}{n-1}\sigma_x^2

の時に実現する。

これが分散の最良推定値の意味である。

Mathematica コード †

DensityPlot への GridLines の入れ方がなかなか分からなかった。
https://mathematica.stackexchange.com/questions/28025/how-to-draw-grid-lines-on-top-of-a-plot

LANG:mathematica
Method -> {"GridLinesInFront" -> True}

がキモになる。

LANG:mathematica
Show[
  DensityPlot[
    1/(Sqrt[2 Pi] s) Exp[-(1 + m^2)/(2 s^2)], {m, -2, 2}, {s, 0, 4}, 
    MeshFunctions -> {#3 &, #3 &}, Mesh -> 10, 
    GridLines -> Automatic, Method -> {"GridLinesInFront" -> True},
    PlotPoints -> 80, ImageSize -> Large
  ],
  Plot[ Sqrt[1 + m^2], {m, -2, 2}, PlotStyle -> {Thick, Red}]
]

コメント・質問 †




添付ファイル: filemu_sigma.jpg 373件 [詳細]
Counter: 6759 (from 2010/06/03), today: 1, yesterday: 4