真値と標準偏差の推定 の履歴(No.5)

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

の確率分布を とし、その期待値を 、分散を とする。

解の測定で、あるデータセット が得られる確率は、

　

である。

このデータセットに対して平均値は だから、 「平均値の期待値」を求めるにはこれに確率を掛けて積分すればよく、

　&math( \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^* );

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

一方、分散の期待値は、

　&math( \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)\\ );

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

　&math( \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\\ );

したがって、

　&math( \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 );

すなわち、データの分散の期待値は真の分散とは等しくならず、 その と一致する。

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

　

とすべきなのである。

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

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

データの確率分布が

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

であれば、あるデータセット が得られる確率は、

　&math( \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] );

である。

この確率を最大化するように を選べば、それらは の最良推定値となるであろう。

[ ] 内は、

　&math( \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\\ );

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

　&math( 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\\ );

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

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

一方、 に対しては、

　&math( \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 );

と置けば、

　&math(

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

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

となって、 で最大値を取ることが分かる。 (上のグラフの赤線が最大値を与える を表す)

それにも関わらず分散の最良推定値が ではなく で与えられるのはどうしてかというと、 では確率が最大になるのが のときであるためだ。

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

　&math( 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]\\ );

この最大値は、

　&math( \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 );

と置いて、

　&math( \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
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}]
]

コメント・質問†