1. Introduction : 平均2乗誤差
母集団のパラメータ$\theta\in\mathbb{R}^d$に対して、標本$X=(X_1,\cdots,X_n)$を用いた推定量$\hat{\theta}=\hat{\theta}(X_1,\cdots,X_n)$を考えたとき、その推定量$\hat{\theta}$の良さを
\begin{eqnarray*}
\mathrm{MSE}(\hat{\theta};\theta)&:=&\mathbb{E}\left[\|\hat{\theta}-\theta\|^2\right]
\end{eqnarray*}
の大きさで測るアイディアは自然に見えるのではないでしょうか。これを推定量$\hat{\theta}$の平均2乗誤差 (mean squared error, MSE) といいます。
2. James-Stein推定量
1956年、Steinが発見したJames-Stein推定量は、平均2乗誤差の小さい推定量が必ずしも直感的なものとは限らないことを教えてくれます。この記事で一緒に考えてみましょう。
問題 : 標本$X=(X_1,\cdots,X_d)$は独立に期待値の異なる正規分布$X_i\sim\mathrm{N}(\mu_i,1^2)$に従って得られるものとします。このとき、期待値ベクトル$\mu=(\mu_1,\cdots,\mu_d)$に対して平均2乗誤差がなるべく小さくなるような推定量を設計してください。
この問題に対する単純な答案の一つは$\hat{\mu}=X$だと思います。おそらく多くの人は$\mu_i$の推定に$X_i$以外の標本点が寄与するとは思わないだろうからです。実際、$\hat{\mu}=X$は多くの良い性質をみたします（証明はAppendix参照）。
- minimax推定量
- 最尤推定量
- 不偏推定量のなかでvarianceが最小
しかし、Steinらは以下の定理が成り立つことに気づきました。
定理 : 以下の定理が成り立つ。
(1) 以下の推定量は、$\hat{\mu}=X$よりも平均2乗誤差が小さい。この推定量$\hat{\mu}^{JS}$をJames-Stein推定量という。
\begin{eqnarray*}
\hat{\mu}^{JS} &:=& \left(1-\frac{d-2}{\|X\|^2}\right)X
\end{eqnarray*}
(2) James-Stein推定量は事前分布$\mu\sim\mathrm{N}(0,\sigma^2E_d)$のもとでの経験Bayes推定量と一致する。
今回の目標は、この意外な定理を証明することです。
2. 定理の証明
定理(1)を証明する際に、以下の恒等式が便利です。事前に証明しておきましょう。
補題 (Steinの恒等式) : $X$を期待値が$\mu$, 分散が$1^2$の正規分布$\mathrm{N}(\mu,1^2)$に従う確率変数、$f(x)$を$\mathrm{N}(\mu,1^2)$の確率密度関数とします。関数$g(x)$が
\begin{eqnarray*}
\lim_{x\rightarrow-\infty}g(x)f(x)=0 &,& \lim_{x\rightarrow\infty}g(x)f(x)=0
\end{eqnarray*}
みたすならば、$\mathbb{E}[g'(X)]=\mathbb{E}[(X-\mu)g(X)]$が成り立つ。
補題の証明 : 部分積分を用いましょう。
\begin{eqnarray*}
\mathbb{E}[g'(X)] &=& \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}g'(x)\exp\left(-\frac{(x-\mu)^2}{2}\right)dx\\
&=& \left[g(x)f(x)\right]_{-\infty}^{\infty}+\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}(x-\mu)g(x)\exp\left(-\frac{(x-\mu)^2}{2}\right)dx
\end{eqnarray*}
ここで、2段目の第1項が仮定から0になるので、$\mathbb{E}[g'(X)]=\mathbb{E}[(X-\mu)g(X)]$が従います。■
2.1 定理(1)の証明
James-Stein推定量の平均2乗誤差を具体的に計算しましょう。
\begin{eqnarray*}
\mathrm{MSE}(\hat{\mu}^{JS},\mu) &=& \mathbb{E}\left[\left\|\left(1-\frac{d-2}{\|X\|^2}\right)X-\mu\right\|^2\right]\\
&=& \mathbb{E}\left[\|x-\mu\|^2\right]-2\mathbb{E}\left[\frac{d-2}{\|X\|^2}X^T(X-\mu)\right]+\mathbb{E}\left[\left(\frac{d-2}{\|X\|^2}\right)^2\|X\|^2\right]\\
&=& \mathbb{E}\left[\|x-\mu\|^2\right]-2(d-2)\mathbb{E}\left[\frac{1}{\|X\|^2}X^T(X-\mu)\right]+(d-2)^2\mathbb{E}\left[\frac{1}{\|X\|^2}\right]\\
\end{eqnarray*}
ここで、3段目の第2項をSteinの恒等式を用いて展開します。
\begin{eqnarray*}
\mathbb{E}\left[\frac{1}{\|X\|^2}X_1(X_1-\mu_1)\right] &=& \mathbb{E}\left[\mathbb{E}\left[\frac{X_1}{X_1^2+\sum_{i=2}^{d}x_i^2}(X_1-\mu_1)\mid X_2=x_2,\cdots,X_d=x_d\right]\right]\\
&=& \mathbb{E}\left[\mathbb{E}\left[\frac{\left(X_1^2+\sum_{i=2}^{d}x_i^2\right)-2X_1^2}{\left(X_1^2+\sum_{i=2}^{d}x_i^2\right)^2}\mid X_2=x_2,\cdots,X_d=x_d\right]\right]\\
&=& \mathbb{E}\left[\frac{\|X\|^2-2X_1^2}{\|X\|^4}\right]
\end{eqnarray*}
が成り立つことから、3段目の第2項は次のようになることが分かるわけです。
\begin{eqnarray*}
\mathbb{E}\left[\frac{1}{\|X\|^2}X^T(X-\mu)\right] &=& \sum_{i=1}^{d}\mathbb{E}\left[\frac{\|X\|^2-2X_i^2}{\|X\|^4}\right]\\
&=& \mathbb{E}\left[\frac{d\|X\|^2-2\|X\|^2}{\|X\|^4}\right]\\
&=& (d-2)\mathbb{E}\left[\frac{1}{\|X\|^2}\right]
\end{eqnarray*}
以上により、
\begin{eqnarray*}
\mathrm{MSE}(\hat{\mu}^{JS},\mu)
&=& \mathbb{E}\left[\|X-\mu\|^2\right]-(d-2)^2\mathbb{E}\left[\frac{1}{\|X\|^2}\right]\\
\end{eqnarray*}
が得られます。ここで、第1項が$\hat{\mu}=X$の平均2乗誤差であることに注意すれば、James-Stein推定量が$\hat{\mu}=X$より小さな平均2乗誤差をもつことが従います。■
2.2 定理(2)の証明
事前分布$\mu\sim\mathrm{N}(0,\sigma^2E_d)$を標本の実現値$X=(x_1,\cdots,x_n)$でBayes更新します。$X\sim\mathrm{N}(\mu,E_d)$なので、事後分布の確率密度関数$f(\mu\mid x)$は
\begin{eqnarray*}
f(\mu\mid x)&\propto&\exp\left(-\frac{\|x-\mu\|^2}{2}\right)\exp\left(-\frac{\|\mu\|^2}{2\sigma^2}\right)\\
&\propto&\exp\left(-\frac{(\sigma^2+1)\|\mu-\frac{\sigma^2}{\sigma^2+1}x\|^2}{2\sigma^2}\right)\\
&=& \exp\left(-\frac{\|\mu-\frac{\sigma^2}{\sigma^2+1}x\|^2}{2\frac{\sigma^2}{\sigma^2+1}}\right)
\end{eqnarray*}
となり、
\mu\mid x\sim\mathrm{N}\left(\frac{\sigma^2}{\sigma^2+1}x, \frac{\sigma^2}{\sigma^2+1}E_d\right)
が成り立ちます。ゆえに、$\mu$のMAP推定量は事前分布の分散パラメータ$\sigma^2$を用いて
\begin{eqnarray*}
\hat{\mu}&=&\frac{\sigma^2}{\sigma^2+1}X\\
&=& \left(1-\frac{1}{\sigma^2+1}\right)X
\end{eqnarray*}
と求まります。ところで、$a=\frac{1}{\sigma^2+1}$は標本から次のように不偏推定できます。事前分布で周辺化した$X$の分布の確率密度関数を$g(x)$と書くと、
\begin{eqnarray*}
g(x)&\propto&\int_{\mathbb{R}^d}\exp\left(-\frac{\|x-\mu\|^2}{2}\right)\exp\left(-\frac{\|\mu\|^2}{2\sigma^2}\right)d\mu\\
&\propto& \int_{\mathbb{R}^d}\exp\left(-\frac{\|x-\mu\|^2}{2}\right)\exp\left(-\frac{\|\mu\|^2}{2\sigma^2}\right)d\mu\\
&\propto&\exp\left(-\frac{\|x\|^2}{2(\sigma^2+1)}\right)
\end{eqnarray*}
ゆえに、$X$の$\mu$の事前分布による周辺分布は$\mathrm{N}\left(0,aE_d\right)$となることがわかります。あとは、逆$\chi^2$-分布の定義とその期待値から
\begin{eqnarray*}
\mathbb{E}\left[\frac{a}{\|X\|^2}\right] &=& \frac{1}{d-2}
\end{eqnarray*}
が従うので、$a$の不偏推定量として$\hat{a}=\frac{d-2}{||x||^2}$を得ます。要するに、$\mu$の経験ベイズ推定量$\hat{\mu}^{EB}$として
\begin{eqnarray*}
\hat{\mu}^{EB} &=& \left(1-\frac{d-2}{\|X\|^2}\right)X
\end{eqnarray*}
を得ることができ、これはJames-Stein推定量$\hat{\mu}^{JS}$と一致しているというわけです。■
Appendix.
$\hat{\mu}=X$がminimax推定量（subsection A）, 最尤推定量（subsection B）, 一様最小分散不偏推定量（subsection C）であることを示しておきましょう。
A. minimax推定量であることの証明
命題1 : $\mu$のminimax推定量を$\hat{\mu}^{MM}$と書くものとする。$\hat{\mu}^{MM}=X$が成り立つ。
命題1の証明 : 推定量を$\hat{\mu}=X$としたとき、平均2乗誤差の上限は$\mathrm{MSE}(X;\mu)=\mathbb{E}\left[||X-\mu||^2\right]=d$なので$\sup_{\mu}\mathrm{MSE}(X;\mu)=d$です。この推定量が$\sup_{\mu}\mathrm{MSE}(\hat{\mu};\mu)$の下限を与えることを示します。
ヒントになるのは、Bayes riskです。どんな$\mu$の事前分布$\pi$に対しても、Bayes riskは
\begin{eqnarray*}
B_{\pi}(\hat{\mu}):=\int_{-\infty}^{\infty}\mathrm{MSE}(\hat{\mu};\mu)\pi(\mu)d\mu &\leq& \sup_{\mu}\mathrm{MSE}(\hat{\mu};\mu)
\end{eqnarray*}
を満たします。さらに左辺の最小値を求めてみましょう。左辺は、
\begin{eqnarray*}
\int_{-\infty}^{\infty}\mathrm{MSE}(\hat{\mu};\mu)\pi(\mu)d\mu &=& \int_{-\infty}^{\infty}\left(\int_{\mathbb{R}^d}\|\hat{\mu}(x)-\mu\|^2f(x;\mu)dx\right)\pi(\mu)d\mu\\
&=& \int_{-\infty}^{\infty}\int_{\mathbb{R}^d}\|\hat{\mu}(x)-\mu\|^2p(x,\mu)d\mu dx\\
&=& \int_{\mathbb{R}^d}\left(\int_{-\infty}^{\infty}\|\hat{\mu}(x)-\mu\|^2\pi(\mu\mid x))d\mu\right)m(x)dx\\
\end{eqnarray*}
と等式変形できます。ここで、$f(x;\mu)$は$X$の確率密度関数、$\pi(\mu)$は$\mu$の事前分布の確率密度関数、$p(x,\mu)$は同時密度関数、$\pi(\mu\mid x)$は$\mu$の事後分布の確率密度関数、$m(x)$は事前分布により周辺化された$X$の確率密度関数です。この式の3段目は、左辺の最小値が$\mu$のEAP推定量によって与えられることを示しています。
例えば、事前分布を$N(0,\sigma^2E_d)$とすれば事後分布は$N\left(\frac{\sigma^2}{\sigma^2+1}x, \frac{\sigma^2}{\sigma^2+1}E_d\right)$となり、特に、
\begin{eqnarray*}
B_{\pi}(\hat{\mu}^{EAP}) &=& \int_{\mathbb{R}^d}\left(\int_{-\infty}^{\infty}\|\hat{\mu}^{EAP}-\mu\|^2\pi(\mu\mid x)d\mu\right)m(x)dx\\ &=& \frac{d\sigma^2}{\sigma^2+1}
\end{eqnarray*}
が成り立ちます。どんな$\sigma^2>0$に対しても、
\begin{eqnarray*}
\sup_{\mu}\mathrm{MSE}(\hat{\mu};\mu) &\geq& \frac{d\sigma^2}{\sigma^2+1}
\end{eqnarray*}
が従うわけなので、特に$\sup_{\mu}\mathrm{MSE}(\hat{\mu};\mu) \geq d$です。以上から、$\hat{\mu}=X$がminimax推定量であることが従います。■
B. 最尤推定量であることの証明
命題2 : $\mu$の最尤推定量を$\hat{\mu}^{ML}$と書くものとする。$\hat{\mu}^{ML}=X$が成り立つ。
命題2の証明 : $\mu$の尤度関数$L(\mu;x_1,\cdots,x_n)$は以下のようになります。
\begin{eqnarray*}
L(\mu;x_1,\cdots,x_n) &=& \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x_i-\mu_i)^2}{2}\right)\times\cdots\times\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x_d-\mu_d)^2}{2}\right)\\
&=& \left(\frac{1}{\sqrt{2\pi}}\right)^d\exp\left(-\frac{\sum_{i=1}^{d}(x_i-\mu_i)^2}{2}\right)
\end{eqnarray*}
特に、$\mu$の対数尤度関数は
\begin{eqnarray*}
l(\mu;x_1,\cdots,x_n) &=& -\frac{1}{2}\sum_{i=1}^{d}(x_i-\mu_i)^2 + \mathrm{const.}
\end{eqnarray*}
となります。最尤推定量は対数尤度関数を最大にするような$\mu$を求めることで得られますが、これは$\hat{\mu}_i^{ML}=x_i$, $(i=1,\cdots,d)$、要するに$\hat{\mu}^{ML}=X$です。■
C. 不偏推定量のなかでvarianceが最小であることの証明
平均2乗誤差はbias-variance分解
\begin{eqnarray*}
\mathrm{MSE}(\hat{\theta};\theta) &=& \|\mathbb{E}[\hat{\theta}]-\theta\|^2+\mathrm{tr}(\mathbb{V}[\hat{\theta}])
\end{eqnarray*}
を持ちます。$\mathrm{bias}(\hat{\theta})=\mathbb{E}[\hat{\theta}]-\theta$を推定量$\hat{\theta}$のbias, $\mathrm{tr}(\mathbb{V}[\hat{\theta}])$を推定量$\hat{\theta}$のvarianceといいます。また、biasが0の推定量を**不偏推定量 (unbiased estimator)**といいます。
命題3 : $\hat{\mu}=X$は$\mu$の不偏推定量のなかで最もvarianceが小さい。
命題3の証明 : $\mu$のFisher情報量行列は$I_{X}(\mu)=E_d$です。多次元の場合のCramer-Raoの定理から、$\hat{\mu}$が不偏推定量のとき$\mathbb{V}[\hat{\mu_i}]\geq1$が成り立ちますが、$\hat{\mu}=X$は不偏推定量でありながらこの等号が成立します。これで、$\hat{\mu}=X$は不偏推定量のなかで最もvarianceが小さいことが示されました。■
Acknowledgement
今週は気分が優れず、あまり記事を書く気が起こらなかったんですが、周囲の応援のお蔭で無事に書くことが出来ました。応援してくれた方々に感謝申し上げます。