真値と標準偏差の推定 のバックアップ差分(No.3)

更新


  • 追加された行はこの色です。
  • 削除された行はこの色です。
[[はじめての誤差論]]
&mathjax();

* 平均値および標準偏差の期待値 [#l5b8c341]

&math(x); の確率分布を &math(p(x)); とすると、
あるデータセット &math(x_1,x_2,\dots,x_n); が得られる確率は、

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

である。

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

 &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\int p(x_1) dx_1\int p(x_2) dx_2\dots\int x_kp(x_k)dx_k\dots\int p(x_n) dx_n \\
&=\frac{1}{n}\sum_{k=1}^n\underbrace{\int x_k p(x_k)dx_k}_{x^*}\\
&=\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(\frac{n-1}{n}); と一致する。

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

 &math(\sigma_x^*=\sqrt\frac{n}{n-1}\sigma_x);

とすべきなのである。

* 具体例(正規分布の時) [#p3014645]

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

データの確率分布が

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

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

 &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(\mu,\sigma); を選べば、それらは &math(x^*,\sigma_x^*); 
の最良推定値となるであろう。

[ ] 内は、

 &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(n=1,\bar x=0,\sigma_x=1); のときこの関数をプロットしてみると次のようになる。

&attachref(mu_sigma.jpg);

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

一方、&math(\sigma); に対しては、

 &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(\mu=\bar x,\sigma=\sigma_x); で最大値を取ることが分かる。

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

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

 &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 コード [#qe24f101]

 LANG:mathematica
 DensityPlot[
   1/(Sqrt[2 Pi] s) Exp[-(1 + m^2)/(2 s^2)], {m, -2, 2}, {s, 0, 4}, 
   MeshFunctions -> {#3 &, #3 &}, Mesh -> 10, 
   PlotPoints -> 80, ImageSize -> Large
 ]
* コメント・質問 [#b447598b]

#article_kcaptcha


Counter: 6936 (from 2010/06/03), today: 1, yesterday: 4