🟠

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

2023/09/30に公開

はじめに

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

導出の流れ

$x \in \mathbb R^p, \beta_0 \in \mathbb R , \beta \in \mathbb R^p$ に対して

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

p=(1+\exp\{-(\beta_0+x^T\beta)\})^{-1}

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

p^y(1-p)^{1-y}

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

尤度は

L = \Pi_{i=1}^N p^{y_i}(1-p)^{1-y_i}

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

つまり

\begin{align*} \hat{\beta}_0,\hat{\beta} &=\argmax_{\beta_0,\beta}L \end{align*}

マイナス対数尤度は

\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*}

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

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

\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*}

$\beta_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を $\beta_0,\beta$ で微分すると

\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)

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

\frac{\partial^2 l}{\partial \beta_j\partial \beta_k}=\sum_{i=1}^Nx_{ij}x_{ik}p_i(1-p_i)

となる。

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

\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}

\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\utilde{\beta}))-\utilde{y}^TX\utilde{\beta}

1回微分は

\dot{l} = X^T(\utilde{p}-\utilde{y})=0

2回微分は

\ddot{l} = X^TVX, V=diag(v_1,\cdots,v_N),v_i=p_i(1-p_i)

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

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

すると、

\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*}

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

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