2016-04-15

1 Empirical Bayes and the James-Stein Estimator:もう一度ぱらぱらめくる『Large-Scale Inference』

何かしらの観測をしたときに、観測値から真値を推定しようとしている。今、何の事前情報もなければ、観測値そのものが真値である尤度が最も高い。単一対象に関する、無情報を前提とした最尤推定値が得られる。
一方、真値について事前分布を置き、観測値と真値のずれについても事前分布を置けば、真値の事後分布が得られる。ベイズ推定値。
では、何かしらの分布を持っていると想定される多数の対象があるとする。それらを(同時に)観測すると、観測値には何かしらの分布が現れる。その分布は、その対象に想定していた真値の分布の様相についての情報を持っている。その真値分布の様相は事前分布として用いることができるから、それに基づいてベイズ推定することができる。これがEmpirical Bayes法。
本章では、真値からのずれを平均0分散1の正規分布とし、多数の対象が平均0分散Aの正規分布としたうえで、上記３つの推定値について検討する。この場合、Empirical Bayes法をJames-Stein推定と言う。
３つの推定方法のどれがよいかは、何かの指標で比較するべきであり、その指標としては、真値からのずれに関する最小二乗和を使っている。
James-Stein推定では、よさの指標が、対象数が２より大きいときには、事前無情報での最尤推定よりも、小さくなることが、数式的に示せる。
例
- 例１
  - 真値分布が $\mu \sim N(0,A)$ 、観測誤差分布が $z-\mu \sim N(0,1)$ であるとき
  - $\hat{\mu}^{(JS)} = 0 + (1-\frac{A}{A+1})z = (1-\frac{N-2}{S})z$ 、ただし $S=\sum z^2$
- 例２
  - 真値分布が $\mu \sim N(0,A)$ 、観測誤差分布が $z-\mu \sim N(0,\sigmaa_0^2)$ であるとき
  - $\hat{\mu}^{(JS)} = 0 + (1-\frac{A}{A+\sigma_0^2})z = (1-\frac{(N-2)\sigma_0^2}{S})z$
- 例３
  - 真値分布が $\mu \sim N(M,A)$ 、観測誤差分布が $z-\mu \sim N(0,\sigma_0^2)$ であるとき
  - $\hat{\mu}^{(JS)} = \bar{z} + (1-\frac{(N-3)\sigma_0^2}{S})(z-\bar{z}_i$ ，ここで $\bar{z}=\frac{\sum z}{N}$
  - 真値分布の中央が不明なときは、その値を観測データから推定する必要があり、そのために、N-2がN-3に変じている。
講義では、ここまでを一気に示し、例を検討することにするのがよく、最後に式変形で納得したい人に、それを課題として出すことにする
- 例４　１次線形回帰への拡張
  - 真値分布が $\mu_i \sim N(M_0 + M_1\times age_i,A)$ 、観測誤差分布が $z_i-\mu_i \sim N(0,\sigma_0^2)$ であるとき
  - $\hat{\mu}_i^{(JS)} = \hat{\mu}_i^{regression} + (1-\frac{(N-4)\sigma_0^2}{S})(z_i-\hat{\mu}_i^{regression}$ ，ここで $\bar{z}=\frac{\sum z}{N}$
1.1 ベイズルールと多変量正規分布推定
- 複数の対象パラメタがある(多変量)
- その複数の対象パラメタの値はある分布に従っている( $\mu \sim g(.)$ )
- その複数の対象パラメタの値を観測するとする。それは $\mu$ に依存したものとなる( $z|\mu$ )
- その観察がある分布に従うとする( $z|\mu \sim f_{\mu}(z)$ )
- ベイズの定理では、 $z$ という観察をしたときに、 $\mu$ の事前確率 $g(\mu)$ が、 $\mu$ の下での $z$ の生起確率 $f_{\mu}(z)$ を使って、 $g(\mu|z) \sim g(\mu) f_{\mu}(z)$ であるとする
- 例
  - $\mu$ はn個の変数であり、それが平均0、分散Aの正規分布を事前分布と考えているとする( $\mu \sim N(0,A)$ )
  - 観察は、真値を平均とした、分散1の正規分布に従うとする( $z|\mu \sim N(\mu,1)$ )
  - このとき、zは平均0、分散A+1の正規分布に従う
  - Rでやってみる

n <- 10^6
A <- runif(1)*5
mu <- rnorm(n,0,sqrt(A))
z <- mu + rnorm(n,0,1)
mean(z)
var(z)
A
var(z) - A

- - 今、興味があるのは、この分散がA+1になっている分布に従っている観測値を見て、それのもともとのパラメタ値 $\mu$ の尤度を知ること
  - ある $\mu_i$ がある $z_j$ を観測値とする確率は $\frac{1}{\sqrt{2\pi}}e^{-\frac{(\mu_i-z_j)^2}{2}}$
  - パラメタ値が $\mu_i$ である事前確率は $\frac{1}{\sqrt{2\pi A}}e^{-\frac{(\mu_i)^2}{2A}}$ であるから
  - $z_j$ をもたらす $_mu_i$ の尤度は、次の値に比例する
  - $\frac{1}{\sqrt{2\pi}}e^{-\frac{(\mu_i-z_j)^2}{2}}\frac{1}{\sqrt{2\pi A}}e^{-\frac{(\mu_i)^2}{2A}}$
  - これを整理すると
  - $\frac{1}{2\pi\sqrt{A}} e^{-\frac{A(\mu_i-z_j)^2 + \mu_i^2}{2A}}$
  - さらに、 $e^{x}$ 以外の成分は、比例することだけを考慮すれば気にしなくて良いことに留意して変形すると
  - $Ce^{-\frac{(\mu_i-\frac{A}{A+1}z_j)^2}{2\frac{A}{A+1}}}$ となる
  - これの意味するところは、 $z_j$ の観測があったとき、それをもたらす $\mu$ の値の事後分布は平均 $\frac{A}{A+1}z_j$ 、分散 $\frac{A}{A+1}$ の正規分布であるということである
  - 今、観察 $z_j$ に対して、 $\mu_j$ の推定量を $z_j$ とするという方法と $z_j\frac{A}{A+1}$ とする方法があるが、すべてのパラメタについて、その真値からのずれを最小二乗和として評価すると $z_j$ そのものとするよりも、 $\frac{A}{A+1}$ で補正した値にする方が、小さくなる(全体としてのフィッティングがよくなる)
  - ずれの総和は、 $z_j$ そのものの場合には、乱雑項を $N(0,1)$ としているから、パラメタ数そのものとなるし、補正すると、その分、小さくなる

m.mle <- z
m.jbz <- z*A/(A+1)

sum((mu-m.mle)^2)
n
sum((mu-m.jbz)^2)
n*A/(A+1)

- Aの値を事前想定してやれば、そのAの値を使って推定できる
1.2 Empirical Bayes Estimation
- この $\mu \frac{A}{A+1}=\mu(1-\frac{1}{A+1})$ であるが、 $\frac{1}{A+1}$ の値を推定してやれば、Aの事前想定値を決めなくても使える
- では、 $\frac{1}{A+1}$ はどうなるか、というと、次のように考える
- 今、 $\mu \sim N(0,A+1)$ なるN次元正規分布なので、 $S=||z||^2$ は、自由度Nの $\chi^2$ 分布と関係する。元の正規分布の分散がA+1なので、 $S \sim (A+1) \chi_N^2$ となる
- 今、Sについて、 $\frac{N-2}{S}$ の期待値が $\frac{1}{A+1}$ であることから、観測データから $\frac{N-2}{S}$ を算出し、 $\frac{1}{A+1}$ の代用とすることができる
- ちなみに、 $\frac{N-2}{S}$ の期待値が $\frac{1}{A+1}$ であることを示しておく
- $\frac{N-2}{S}$ の期待値は $\int \frac{N-2}{x(A+1)} \frac{1}{\frac{N}{2} \Gamma{\frac{N}{2}}} x^{\frac{N}{2}-1}e^{-\frac{x}{2}} dx$ は $\chi^2$ 分布に従う変数の期待値の定義式
- これを $M=N-2,N=M+2$ を用いて変形する。そのとき、 $\Gamma(a+1) = \Gamma(a)\times a$ であることと、 $\chi^2$ 分布の確率密度の積分が１になることを使う
- $=\int \frac{M}{x(A+1)} \frac{1}{\frac{M+2}{2} \Gamma{\frac{M+2}{2}}} x^{\frac{M+2}{2}-1}e^{-\frac{x}{2}} dx$
- $=\frac{M}{2\frac{M}{2}(A+1)}\int \frac{1}{\frac{1}{M}\Gamma{\frac{M}{2}}} x^{\frac{M}{2}-1}e^{-\frac{x}{2}} dx$
- $=\frac{1}{A+1}$
- Rでもやっておこう

n <- 10^5
S <- (A+1) * rchisq(n,n)
mean((n-2)/S)
1/(A+1)

- 今、n個のパラメタについて観測値が得られたときにとしたとき、
  - $\hat{\mu}^{JS} = (1-\frac{n-2}{S})z$ をJames-Stein estimatorと呼ぶ
３つの推定値
- 観測値を見て、他の要素とは全く独立に、観測誤差が正規乱数になっていると考えれば、z自体を推定値とするし、そのようなとき、それは最尤推定値になっている
- たくさんの観測値がえられるときに、もとのパラメタが依存する分散と、観測誤差の分散とを事前に想定してやれば、それを用いて、補正することができる。これは、事前想定を入れるから、ベイズ法
- そして、観測値からJames-Stein式に $\frac{1}{A+1}$ を「経験的に」求め、それを持って補正するのがEmpirical Bayes法
- 真値からのずれを、全てのパラメタの最小二乗で測れば、Aが正しければベイズ法が最小。James-Steinは最尤推定よりは小さくなる