最尤推定値の計算法

数理統計　第11章 §11.3　★★ 尤度方程式を反復で解く。Newton 法は2次収束だが初期値に敏感、EM アルゴリズムは尤度を単調に増やすが1次収束。どちらも更新式の導出が小問になる

Newton–Raphson 法

対数尤度 $\ell(\theta)$（1母数）の最大点を、尤度方程式 $\ell'(\theta)=0$ の解として求める。現在の値 $\theta_k$ のまわりで $\ell'$ を1次まで Taylor 展開すると

$$\ell'(\theta)\approx\ell'(\theta_k)+\ell''(\theta_k)(\theta-\theta_k)$$

である。右辺を 0 とおいて解いた点を次の値にする。

$$\theta_{k+1}=\theta_k-\frac{\ell'(\theta_k)}{\ell''(\theta_k)} .\tag{11.12}$$

これを Newton–Raphson 法（Newton 法）と呼ぶ。$\ell$ を $\theta_k$ のまわりで2次式で近似し、その2次式の最大点へ移る、と言い換えることもできる。$\ell''(\theta_k)<0$ ならその2次式は上に凸で、移る先は最大点である。$\ell''(\theta_k)>0$ のところでは2次式が下に凸になり、Newton 法は極小点へ向かう。

母数がベクトル $\bm\theta\in\mathbb{R}^p$ のときは、次の2つを使う。

$$\bm s(\bm\theta)=\frac{\partial\ell}{\partial\bm\theta}\ \ (\text{スコアベクトル}),\qquad J(\bm\theta)=-\frac{\partial^2\ell}{\partial\bm\theta\partial\bm\theta^\top}\ \ (\text{観測情報行列})$$

更新式は

$$\bm\theta_{k+1}=\bm\theta_k+J(\bm\theta_k)^{-1}\bm s(\bm\theta_k)\tag{11.13}$$

とする。$J(\hat{\bm\theta})^{-1}$ は最尤推定量の漸近分散の推定にもなる（第6章）ので、収束したときに標準誤差も同時に得られる。

証明

$\psi$ を $\theta_k$ のまわりで、2次の剰余つきで展開する。$\hat\theta$ で評価すると次を得る。

$$0=\psi(\hat\theta)=\psi(\theta_k)+\psi'(\theta_k)(\hat\theta-\theta_k)+\frac12\psi''(\xi_k)(\hat\theta-\theta_k)^2$$

両辺を $\psi'(\theta_k)$ で割り、Newton の式 (11.12) の $\theta_{k+1}=\theta_k-\psi(\theta_k)/\psi'(\theta_k)$ を使うと

$$0=\frac{\psi(\theta_k)}{\psi'(\theta_k)}-e_k+\frac{\psi''(\xi_k)}{2\psi'(\theta_k)}e_k^2=(\theta_k-\theta_{k+1})-e_k+\frac{\psi''(\xi_k)}{2\psi'(\theta_k)}e_k^2 =-e_{k+1}+\frac{\psi''(\xi_k)}{2\psi'(\theta_k)}e_k^2$$

となり (11.14) を得る。

次に収束を示す。$\hat\theta$ の近くの区間 $I$ をとり、その上で次が成り立つとする。

$$\lvert \psi''\rvert \le A,\qquad\lvert \psi'\rvert \ge B>0,\qquad C=\frac A{2B}$$

$\theta_k\in I$ かつ $C\lvert e_k\rvert <1$ とする。このとき

$$\lvert e_{k+1}\rvert \le C\lvert e_k\rvert ^2<\lvert e_k\rvert$$

なので、$\theta_{k+1}$ も $I$ に入る。帰納的に

$$C\lvert e_k\rvert \le(C\lvert e_0\rvert )^{2^k}\to0$$

となる。極限の比は $\xi_k\to\hat\theta$、$\theta_k\to\hat\theta$ と連続性から従う。

定理は「十分近くから始めれば」という条件つきである。初期値が遠いと、$\ell''>0$ の領域に入って極小点へ向かったり、発散したりする（確認テスト 11.3-5）。実際には、モーメント推定量のような一致推定量を初期値に使う。$\sqrt n$ 一致推定量（$\sqrt n(\tilde\theta-\theta)$ が確率有界になる推定量）から Newton 法を1回だけ行った推定量（ワンステップ推定量）は、最尤推定量と同じ漸近分布を持つことが知られている。

Fisher のスコア法

Newton 法の $-\ell''(\theta_k)$（観測情報量）を、その期待値である Fisher 情報量 $I_n(\theta_k)=nI(\theta_k)$ に置き換えたもの

$$\theta_{k+1}=\theta_k+\frac{\ell'(\theta_k)}{I_n(\theta_k)}\tag{11.15}$$

を Fisher のスコア法と呼ぶ。$I_n(\theta)>0$ なので、$\ell''$ の符号で向きが反転することがなく、Newton 法より安定する。

収束の速さは一般に1次に落ちる。しかし $-\ell''(\hat\theta)/n$ と $I(\hat\theta)$ の差は、標本が大きいと小さい。そのため実際の反復回数はあまり変わらない。

自然母数で書いた指数型分布族 $f(x\mid\theta)=h(x)\exp\{\theta T(x)-A(\theta)\}$ では、$\ell''(\theta)=-nA''(\theta)$ がデータによらないので、観測情報量とFisher 情報量が一致し、スコア法と Newton 法は同じ反復になる（確認テスト 11.3-1）。第9章のロジスティック回帰の最尤推定でスコア法を使うと、各反復が重み付き最小2乗法の形になる。これを反復重み付き最小2乗法（IRLS）と呼ぶ。

EM アルゴリズム

観測データ $\bm y$ だけでは尤度が複雑でも、観測されなかったデータ $\bm z$ があれば尤度が簡単になる、という場合がある。典型的な例は次の3つである。

• 打ち切りのある寿命データ：打ち切られた個体の本当の寿命が $\bm z$。
• 混合分布：各観測がどの成分から来たかのラベルが $\bm z$。
• いくつかのセルがまとめて観測された多項分布：まとめられる前の内訳が $\bm z$。

$(\bm y,\bm z)$ を完全データと呼び、その密度を $f_c(\bm y,\bm z\mid\theta)$ と書く。完全データの対数尤度を $\ell_c(\theta)=\log f_c(\bm y,\bm z\mid\theta)$ と書く。観測データの尤度と対数尤度は

$$L(\theta)=\int f_c(\bm y,\bm z\mid\theta)\,\mathrm{d}\bm z,\qquad \ell(\theta)=\log L(\theta)$$

である。

名前は Expectation（期待値）と Maximization（最大化）の頭文字である（Dempster–Laird–Rubin が一般の形にまとめた）。完全データの分布が、指数型分布族

$$\exp\{\theta^\top T(\bm y,\bm z)-A(\theta)\}$$

の形だとする。このとき $\ell_c$ は十分統計量 $T$ について1次式なので、計算は次のように簡単になる。

• E ステップは、$T$ の条件付き期待値 $T_k=E_{\theta_k}[T\mid\bm y]$ を計算するだけでよい。
• M ステップは、完全データの最尤推定量の式の $T$ に $T_k$ を入れるだけになる。

証明

まず次の分解から始める。

$$f_c(\bm y,\bm z\mid\theta)=L(\theta)\,k(\bm z\mid\bm y,\theta)$$

この対数をとり、両辺を $k(\bm z\mid\bm y,\theta_k)$ で積分する。$\ell(\theta)$ は $\bm z$ によらないので

$$\ell(\theta)=Q(\theta\mid\theta_k)-H(\theta\mid\theta_k),\qquad H(\theta\mid\theta_k)=\int\log k(\bm z\mid\bm y,\theta)\;k(\bm z\mid\bm y,\theta_k)\,\mathrm{d}\bm z\tag{11.17}$$

を得る。$H$ について、Jensen の不等式（$\log$ は凹）から

$$H(\theta\mid\theta_k)-H(\theta_k\mid\theta_k)=\int\log\frac{k(\bm z\mid\bm y,\theta)}{k(\bm z\mid\bm y,\theta_k)}\,k(\bm z\mid\bm y,\theta_k)\,\mathrm{d}\bm z \le\log\int k(\bm z\mid\bm y,\theta)\,\mathrm{d}\bm z=\log1=0$$

である（左辺は Kullback–Leibler 情報量に負号を付けたもの）。したがって

$$\ell(\theta_{k+1})-\ell(\theta_k)=\bigl\{Q(\theta_{k+1}\mid\theta_k)-Q(\theta_k\mid\theta_k)\bigr\}-\bigl\{H(\theta_{k+1}\mid\theta_k)-H(\theta_k\mid\theta_k)\bigr\}\ge0 .$$

第1の括弧は M ステップの定義から 0 以上、第2の括弧は今示したとおり 0 以下である。

証明からわかるとおり、M ステップで $Q$ を最大化しきらなくても単調性は保たれる。$Q(\theta_{k+1}\mid\theta_k)\ge Q(\theta_k\mid\theta_k)$ となる $\theta_{k+1}$ を選べばよい（一般化 EM）。

また (11.17) を $\theta$ で微分すると、次が分かる。$H(\cdot\mid\theta_k)$ は $\theta=\theta_k$ で最大なので、

$$\ell'(\theta_k)=\frac{\partial Q(\theta\mid\theta_k)}{\partial\theta}\Bigr|_{\theta=\theta_k}$$

である。

EM の不動点 $\theta^\ast$ をとる。つまり $\theta_k=\theta^\ast$ なら $\theta_{k+1}=\theta^\ast$ となる点である。

これが $Q(\cdot\mid\theta^\ast)$ の内点の最大点なら、そこで $\partial Q/\partial\theta=0$ である。すなわち $\ell'(\theta^\ast)=0$ となり、不動点は尤度方程式の解である。

ただし単調性が保証するのは尤度が増え続けることだけで、大域的な最大点への収束は保証されない。尤度に複数の山がある混合分布では、初期値を変えて何回か走らせる。

収束の速さと欠測情報

例題 11.5 の収束率 $(n-r)/n$ は、一般の EM アルゴリズムの性質の特別な場合である。(11.17) を2回微分して $\hat\theta$ で評価する。次の2つの情報量を定める。

$$I_c=-\frac{\partial^2Q(\theta\mid\hat\theta)}{\partial\theta^2}\Bigr|_{\hat\theta}\ (\text{完全データの情報量}),\qquad I_m=-\frac{\partial^2H(\theta\mid\hat\theta)}{\partial\theta^2}\Bigr|_{\hat\theta}\ (\text{欠測データの情報量})$$

観測データの情報量は、その差 $-\ell''(\hat\theta)=I_c-I_m$ になる。EM の収束率（更新の写像の $\hat\theta$ での微分）は

$$g'(\hat\theta)=\frac{I_m}{I_c}=1-\frac{-\ell''(\hat\theta)}{I_c}\tag{11.18}$$

で、完全データの情報のうち欠測で失われた割合に等しい。例題では $I_c=n/\hat\lambda^2$、$-\ell''(\hat\lambda)=r/\hat\lambda^2$ である。よって $g'=1-r/n=(n-r)/n$ となり、(3) と一致する。

図 11.2 のとおり、EM は1次収束で Newton 法より遅いが、尤度を減らさないので初期値が遠くても破綻しない。Newton 法は2次収束だが、この例では初期値 $\lambda_0=1$ から始めると最初の1回で負の値に飛んで発散する。実務では、EM で解の近くまで進めてから Newton 法に切り替えることもある。また EM では $\ell''$ が計算されないので、最尤推定量の標準誤差は収束後に別途 (11.18) の関係などを使って求める必要がある。

有限混合分布

EM アルゴリズムが最もよく使われる場面は混合分布である。密度が

$$f(x\mid\bm\theta)=\pi\,\phi(x;\mu_1,\sigma_1^2)+(1-\pi)\,\phi(x;\mu_2,\sigma_2^2)$$

の2成分の正規混合を考える（$\phi(\cdot;\mu,\sigma^2)$ は正規密度）。対数尤度 $\sum_i\log f(x_i\mid\bm\theta)$ では対数の中に和があるので、尤度方程式は解けない。各観測が成分1から来たかどうかを表すラベル $z_i\in\{0,1\}$（成分1なら 1）を欠測データとすると、完全データの対数尤度は

$$\ell_c=\sum_{i=1}^n\Bigl[z_i\{\log\pi+\log\phi(x_i;\mu_1,\sigma_1^2)\}+(1-z_i)\{\log(1-\pi)+\log\phi(x_i;\mu_2,\sigma_2^2)\}\Bigr]$$

で、$z_i$ について1次式である。E ステップでは $z_i$ をその条件付き期待値（Bayes の定理による事後確率）

$$w_i^{(k)}=E_{\bm\theta_k}[z_i\mid x_i]=\frac{\pi_k\phi(x_i;\mu_{1,k},\sigma_{1,k}^2)}{\pi_k\phi(x_i;\mu_{1,k},\sigma_{1,k}^2)+(1-\pi_k)\phi(x_i;\mu_{2,k},\sigma_{2,k}^2)}\tag{11.19}$$

で置き換える。M ステップは、$w_i^{(k)}$ を重みとした重み付きの最尤推定で

$$\pi_{k+1}=\frac1n\sum_iw_i^{(k)},\qquad \mu_{1,k+1}=\frac{\sum_iw_i^{(k)}x_i}{\sum_iw_i^{(k)}},\qquad \sigma_{1,k+1}^2=\frac{\sum_iw_i^{(k)}(x_i-\mu_{1,k+1})^2}{\sum_iw_i^{(k)}}$$

となる（成分2は $w_i^{(k)}$ を $1-w_i^{(k)}$ に替える）。導出は確認テスト 11.3-4 で行う。

なお、2つの分散を別々に推定する正規混合では、大域的な最大点は存在しない。$\mu_1=x_i$ として $\sigma_1\to0$ とすると、尤度が無限大に発散するからである。EM で求めるのは尤度方程式の解のうち、分散が 0 に退化していない局所的な最大点である。

確認テスト

オリジナル

1. 11.3-1 自然母数の1母数指数型分布族

$$f(x\mid\theta)=h(x)\exp\{\theta T(x)-A(\theta)\}$$

から $n$ 個の観測を得る。観測情報量 $-\ell''(\theta)$ が Fisher 情報量 $I_n(\theta)$ に一致することを示せ。したがって Newton 法とスコア法が同じ反復になることを説明せよ。2. 11.3-2 例題 11.5 の対数尤度 $\ell(\lambda)=r\log\lambda-\lambda T$ に Newton 法を使う。更新式が $\lambda_{k+1}=2\lambda_k-\lambda_k^2T/r$ になることを示せ。誤差 $e_k=\lambda_k-\hat\lambda$ が $e_{k+1}=-(T/r)e_k^2$ を厳密に満たすことを示せ。$0<\lambda_0<2\hat\lambda$ のとき、かつそのときに限り $\hat\lambda$ に収束することを示せ。3. 11.3-3 4つのセルの多項分布を考える。セルの確率は

$$\Bigl(\frac12+\frac\theta4,\ \frac{1-\theta}4,\ \frac{1-\theta}4,\ \frac\theta4\Bigr)\qquad(0<\theta<1)$$

で、観測度数は $(y_1,y_2,y_3,y_4)=(80,25,21,24)$ とする。第1セルを確率 $\frac12$ の部分と $\frac\theta4$ の部分に分けた内訳を欠測データとして、EM の更新式を導け。$\theta_0=0.5$ から $\theta_1$ を求めよ。また、尤度方程式を2次方程式に直して最尤推定値を求めよ。4. 11.3-4 2成分の正規混合で、両成分の分散が 1 と既知、未知母数が $(\pi,\mu_1,\mu_2)$ であるとする。M ステップの更新式

$$\pi_{k+1}=\frac1n\sum w_i^{(k)},\qquad \mu_{1,k+1}=\frac{\sum w_i^{(k)}x_i}{\sum w_i^{(k)}}$$

を $Q$ の最大化から導け。5. 11.3-5 位置母数の Cauchy 分布 $f(x\mid\theta)=\{\pi(1+(x-\theta)^2)\}^{-1}$ から1個の観測 $x=0$ を得た。Newton 法の更新式が

$$\theta_{k+1}=-\frac{2\theta_k^3}{1-\theta_k^2}$$

になることを示せ。さらに次を示せ。

• $\lvert \theta_0\rvert <1/\sqrt3$ なら $\hat\theta=0$ に収束する。
• $\lvert \theta_0\rvert >1/\sqrt3$（$\lvert \theta_0\rvert \ne1$）なら収束しない。