共分散0の多変量正規分布の母平均ベクトルのJames-Stein推定量は、誤差二乗和をリスクとしたとき、標本平均に優越している。これはあっちこっちに書いてある有名な話だが、証明をいくら読んでもピンとこなかった。
このたび業を煮やし、真剣に考えた。
まず、次の補題を示しておく。これはSteinの補題と呼ばれている。
補題. \(Z \sim N(\mu, 1)\)とする。期待値が存在し、$$ \lim_{z \rightarrow \pm \infty} h(z) \exp \left( – \frac{1}{2}(z-\mu)^2 \right) = 0$$ ならば、\(h'(Z) = \partial h(Z)/\partial Z \) として$$ E[h(Z) (Z-\mu)] = E[h'(Z)]$$ である。
証明: 正規分布の密度関数より $$ E[h(Z)(Z-\mu)] = \frac{1}{\sqrt{2\pi}} \int^\infty_{-\infty} h(z) (z-\mu) \exp\left\{ -\frac{1}{2} (z-\mu)^2 \right\} dz$$ 指数関数の部分を\(z\)で微分してみよう。覚えてますか、合成関数の微分。\(x = u(z) = -\frac{1}{2}(z-\mu)^2\)と略記して、$$ \frac{d}{dz} \exp(u(z)) = \frac{d}{dx} \exp(x) \times \frac{d}{dz} u(z) $$ 私の青春時代は戦後の混乱期で微分公式も習わなかったのだが(適当)、なんでも\( (e^x)’ = e^x \)なのだそうでして、第1項は\(\exp(x)\)のままになる。第2項は \( \frac{d}{dz} u(z) = – \frac{1}{2} \times 2(z-u) = -(z-u) \) なので、$$ \frac{d}{dz} \exp(u(z)) = – (z-\mu) \exp \left\{-\frac{1}{2}(z-\mu)^2\right\} $$ これを元の式に代入して、$$ E[h(Z)(Z-\mu)] = – \frac{1}{\sqrt{2\pi}} \int^\infty_{-\infty} h(z) \frac{d}{dz} \exp(u(z)) dz$$ なんでこんなへんなことをしておるかというと、部分積分したいのである。覚えてますか、$$ \int_a^b f(z)g'(z)dz = \left. f(z)g(z) \right|_a^b – \int_a^b f'(z)g(z) dz $$ですよね。 \(f(z) = h(z), g'(z) = \frac{d}{dz} \exp(u(z))\)とすると $$ = – \frac{1}{\sqrt{2\pi}} \left. h(z) \exp(u(z)) \right|_{-\infty}^{\infty} + \frac{1}{\sqrt{2\pi}} \int^\infty_{-\infty} h'(z) \exp(u(z)) dz$$ 仮定より第1項はゼロ。第2項をよくみると、結局 $$ = E[h'(Z)] $$ である。
なお、この補題には \( \lim_{z \rightarrow \pm \infty} h(z) \exp \left( – \frac{1}{2}(z-\mu)^2 \right) = 0 \) というへんな条件が出てくるが、これは要するに、\(h(z)\)はフツウの関数だよ、なんか指数関数的に増減したりはしないよ、だからそれに正規密度を掛ければ右にいっても左に行っても0に近づくよ、ということであろうと思う。
いよいよ、定理の証明である。
定理. \(i=1, \ldots, m, m\geq 3\)について\(Z_i \sim_{ind} N(\mu_i, 1)\)とし、\(\mathbf{\mu} = (\mu_1, \ldots, \mu_m)^\top, \mathbf{Z} = (Z_1, \ldots, Z_m)^\top\)とする。$$ \hat{\mu}_{JS} (a) = \left( 1- \frac{a}{||\mathbf{Z}||^2} \right) \mathbf{Z} $$ ただし\(0 \lt a \lt 2(m-2)\)とする。
リスクを全MSEとしたとき、\(\hat{\mu}_{JS}(a)\)は\(Z\)に優越する。最適点は\(a=m-2\)である。
証明: $$ R(\mu, \hat{\mu}_{JS}(a)) = E \left[ || \hat{\mu}_{JS}(a) – \mu ||^2 \right] = E\left[ || \mathbf{Z} – \mu – \frac{a \mathbf{Z}}{||\mathbf{Z}||^2} ||^2 \right] $$ $$ = E \left[ || \mathbf{Z}-\mu ||^2 \right] + a^2 E \left[ \frac{1}{||\mathbf{Z}||^2} \right] – 2a E \left[ \frac{\mathbf{Z}^\top (\mathbf{Z} -\mu)}{||\mathbf{Z}||^2} \right]$$ 第1項は分散の和、つまり\(m\)である。第2項もこれでよいとしよう。問題は第3項である。期待値のところに注目すると、$$ E \left[ \frac{\mathbf{Z}^\top (\mathbf{Z} -\mu)}{||\mathbf{Z}||^2} \right] = \sum_{i=1}^m E \left[ \frac{Z_i (Z_i -\mu_i)}{||\mathbf{Z}||^2} \right] $$ \(h(Z_i) = \frac{Z_i}{|| \mathbf{Z} ||^2} \)と置くと $$ = \sum_{i=1}^m E \left[ h(Z_i)(Z_i-\mu_i) \right] $$ 補題より $$ = \sum_{i=1}^m E \left[ \frac{\partial}{\partial Z_i} \frac{Z_i}{||\mathbf{Z}||^2} \right] $$ 覚えてますか、分数の微分を。$$ \frac{d}{dx} \frac{f(x)}{g(x)} = \frac{f'(x)g(x) – f(x)g'(x)}{ (g(x))^2 } $$ですよね。$$ f(Z_i) = Z_i, \ f'(Z_i) = 1$$ $$ g(Z_i) = || \mathbf{Z} || ^2 = Z_1^2 + \ldots + Z_m^2, \ g'(Z_i) = 2Z_i $$ とおいて$$ = \sum_{i=1}^m E \left[ \frac{1 \times ||\mathbf{Z}||^2 – Z_i \times 2Z_i}{||\mathbf{Z}||^4} \right] = m E \left[ \frac{1}{||\mathbf{Z}||^2} \right] – \sum_{i=1}^m E \left[ \frac{2Z_i^2}{||\mathbf{Z}||^4} \right] = (m-2) E \left[ \frac{1}{||\mathbf{Z}||^2} \right] $$ というわけで、$$ R(\mu, \hat{\mu}_{JS}(a)) = m + [a^2 – 2a(m-2)] E \left[ \frac{1}{||\mathbf{Z}||^2} \right] $$ \(a^2 – 2a(m-2)\)は下に凸な放物線で、\(0 \lt a \lt 2(m-2)\)の範囲で負であり、\(a=m-2\)で最小とある。つまり、\(R(\mu, \hat{\mu}_{JS}(a))\)は常に\(R(\mu, Z) = m\)より小さく、\(a = m-2\)のときに最適となる。
やれやれ、やっとわかったような気がする…
以上、Rao & Molina (2015) “Small Area Estimation”, 2nd ed. の3.5節を参考にしたが、全ての誤りはもちろん私のせいである。