ロジスティック回帰の最尤推定量の導出

はじめに

今回は前回の続き、ロジスティック回帰の最尤推定の証明をしたいと思います。
数式ばかりになってしまいますのでご了承ください。

導出の流れ

1. 尤度を求める
2. マイナス対数尤度に変換
3. $β\beta$で微分
4. 見やすくするために表記を全てベクトルに
5. ニュートンラフソン法により推定値を求める。

前回の記事
https://zenn.dev/totopironote/articles/48e39f9b7fef2b

本題

$x∈Rp,β0∈R,β∈Rpx \in \mathbb R^p, \beta_0 \in \mathbb R , \beta \in \mathbb R^p$ に対して

y=1となる確率pは以下で定義される

$$p=(1+exp⁡{−(β0+xTβ)})−1 p=(1+\exp\{-(\beta_0+x^T\beta)\})^{-1}$$

ベルヌーイ分布の密度関数は

$$py(1−p)1−y p^y(1-p)^{1-y}$$

よってデータ$(x1,y1),⋯ ,(xN,yN)∈Rp×{0,1}(x_1,y_1),\cdots,(x_N,y_N) \in \mathbb R^p \times \{0,1\}$が得られたとき。

尤度は

$$L=Πi=1Npyi(1−p)1−yi L = \Pi_{i=1}^N p^{y_i}(1-p)^{1-y_i}$$

最優推定量はこの尤度を最大化する$(β0^,β^)(\hat{\beta_0},\hat{\beta})$を求めたい。

つまり

$$β^0,β^=arg max⁡β0,βL \begin{align*} \hat{\beta}_0,\hat{\beta} &=\argmax_{\beta_0,\beta}L \end{align*}$$

マイナス対数尤度は

$$l=−log⁡L=−∑i=1N{ylog⁡p+(1−y)log⁡(1−p)}=−∑i=1N{ylog⁡p1−p+log⁡(1−p)}=∑i=1N{log⁡(1+exp⁡(β0+xTβ))−y(β0+xTβ)} \begin{align*} l = -\log L &= -\sum_{i=1}^N\{y\log p + (1-y)\log(1-p)\}\\ &=-\sum_{i=1}^N\{y\log\frac{p}{1-p}+\log(1-p)\}\\ &=\sum_{i=1}^N\{\log(1+\exp(\beta_0+x^T\beta))-y(\beta_0+x^T\beta)\} \end{align*}$$

ここでは尤度を最大化ではなくマイナス対数尤度の最小化を考える。

$βj,j=1,⋯ ,p\beta_j ,j=1,\cdots,p$ に関する微分は

$$∂l∂βj=∑i=1N{xijexp⁡(β0+xiTβ)1+exp⁡(β0+xiTβ)−xijyi}=∑i=1N{xij(pi−yi)}=0 \begin{align*} \frac{\partial l}{\partial \beta_j} &=\sum_{i=1}^N \{\frac{x_{ij}\exp(\beta_0+x_i^T\beta)}{1+\exp(\beta_0+x_i^T\beta)}-x_{ij}y_i\}\\ &=\sum_{i=1}^N\{x_{ij}(p_i-y_i)\}=0 \end{align*}$$

$β0\beta_0$に関する微分は

$$∂l∂β0=∑i=1N{exp⁡(β0+xiTβ)1+exp⁡(β0+xiTβ)−yi}=∑i=1N(pi−yi)=0 \begin{align*} \frac{\partial l}{\partial \beta_0} &=\sum_{i=1}^N \{\frac{\exp(\beta_0+x_i^T\beta)}{1+\exp(\beta_0+x_i^T\beta)}-y_i\}\\ &=\sum_{i=1}^N(p_i-y_i)=0 \end{align*}$$

pを$β0,β\beta_0,\beta$で微分すると

$$∂pi∂β0=pi(1−pi),∂pi∂βj=xijpi(1−pi) \frac{\partial p_i}{\partial \beta_0} = p_i(1-p_i) , \frac{\partial p_i}{\partial \beta_j} = x_{ij}p_i(1-p_i)$$

よって$xi0=1x_{i0}=1$として、$j,k=0,⋯ ,pj,k=0,\cdots,p$に対して

$$∂2l∂βj∂βk=∑i=1Nxijxikpi(1−pi) \frac{\partial^2 l}{\partial \beta_j\partial \beta_k}=\sum_{i=1}^Nx_{ij}x_{ik}p_i(1-p_i)$$

となる。

ここで全てベクトル表記になおす。

$$β~=(β0β),X==(1x1T1x2T⋮⋮1xNT)=(1x11⋯x1p1x21⋯x2p⋮1xN1⋯xNp),y~=(y1⋮yN) \utilde{\beta}=\begin{pmatrix} \beta_0\\ \beta \end{pmatrix}, X==\begin{pmatrix} 1 &x_{1}^T \\ 1 &x_{2}^T \\ \vdots & \vdots \\ 1 &x_{N}^T \end{pmatrix}=\begin{pmatrix} 1 &x_{11} & \cdots &x_{1p}\\ 1 &x_{21} & \cdots & x_{2p}\\ & & \vdots\\ 1 &x_{N1} & \cdots & x_{Np} \end{pmatrix},\utilde{y}=\begin{pmatrix} y_1\\ \vdots\\ y_N \end{pmatrix}$$

$$p~=(p1⋮pN),log⁡p~1N−p~=(log⁡p11−p1⋮log⁡pN1−pN)=Xβ~ \utilde{p} = \begin{pmatrix} p_1\\ \vdots\\ p_N \end{pmatrix}, \log\frac{\utilde{p}}{1_N-\utilde{p}}=\begin{pmatrix} \log\frac{p_1}{1-p_1}\\ \vdots\\ \log\frac{p_N}{1-p_N} \end{pmatrix}=X\utilde{\beta}$$

これらを用いると、

マイナス対数尤度は

$$l=log⁡(1+exp⁡(Xβ~))−y~TXβ~ l = \log(1+\exp(X\utilde{\beta}))-\utilde{y}^TX\utilde{\beta}$$

1回微分は

$$l˙=XT(p~−y~)=0 \dot{l} = X^T(\utilde{p}-\utilde{y})=0$$

2回微分は

$$l¨=XTVX,V=diag(v1,⋯ ,vN),vi=pi(1−pi) \ddot{l} = X^TVX, V=diag(v_1,\cdots,v_N),v_i=p_i(1-p_i)$$

ここでdiagは対角に成分を並べたもので非対角は0の行列。

解析的には解けないのでニュートンラフソン法を使う。

すると、

$$β~new=β~old−l¨−1l˙=β~old−(XTVX)−1XT(p~−y~) \begin{align*} \utilde{\beta}^{new}&=\utilde{\beta}^{old}-\ddot{l}^{-1}\dot{l}\\ &=\utilde{\beta}^{old}-(X^TVX)^{-1}X^T(\utilde{p}-\utilde{y}) \end{align*}$$

この更新式で収束したものが最尤推定量である。

まとめ

ロジスティック回帰の場合は最尤推定量が解析的に求められないので、ニュートンラフソンを使いました。最後の式はかなりスッキリしていて気持ちい。

次回はロジスティック回帰を回帰の文脈で用いる方法を説明します。
https://zenn.dev/totopironote/articles/5af3ece45d144c

参考文献

• 鈴木 譲 「統計的機械学習の数理100問 with Python」