真値と標準偏差の推定

更新

はじめての誤差論

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

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

この系に対して n 回の測定を行った時の平均値や標準偏差がどのような確率分布を持つか、これらの期待値がいくつになるか、を考えよう。(母集団から n 個の標本を取り出したと考えてもよい)

確率分布 p(x) を使うと、あるデータセット 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 1032件 [詳細]
Counter: 15463 (from 2010/06/03), today: 2, yesterday: 6