一般化線形モデル（特にロジスティック回帰）

 このスクラップは？一般化線形モデルのお勉強ノート。
主に、以下の書籍で学んだ内容をまとめて行く予定：
A. J. Dobson（田中 豊、森川 敏彦、山中 竹春、冨田 誠 翻訳）、一般化線形モデル入門 原著第2版（共立出版、2008）

https://www.kyoritsu-pub.co.jp/book/b10010684.html

https://amzn.asia/d/ahWa303

畳屋民也

 一般化線形モデルの最尤推定
 問題設定観測値 \left\{(\boldsymbol{x}_i, y_i)\right\}_{i=1}^n が与えられているものとする。
各サンプルは独立同時な過程により生成され、説明変数 \boldsymbol{x} が与えられた時の目的変数 y の（条件付き）確率密度分布関数はパラメータ \boldsymbol{\beta} を用いて f(y\vert \boldsymbol{x}; \boldsymbol{\beta}) と表せるものとする。
ここで、説明変数 \boldsymbol{x} で条件づけられた y の期待値を \mu = E_Y[y\vert \boldsymbol{x}; \boldsymbol{\beta}] と置いたとき、\mu と \boldsymbol{x}, \, \boldsymbol{\beta} との間に適当な関数 g(\mu) を用いて以下のような関係があるものとする：

g(\mu) = \boldsymbol{x}^\top \boldsymbol{\beta}.
Logistic 回帰の場合Logistic 回帰の場合、\mu と \boldsymbol{x}, \, \boldsymbol{\beta} の間に以下のような関係を仮定する：

\mu = \frac{1}{1 + e^{-\boldsymbol{x}^\top \boldsymbol{\beta}}}.
これは

\log\left(\frac{\mu}{1 + \mu}\right) = \boldsymbol{x}^\top \boldsymbol{\beta}
と書き換えられ、

g(\mu) = \log\left(\frac{\mu}{1 + \mu}\right)
に対応する。
このとき、尤度を最大化するパラメータ \boldsymbol{\beta} = \hat{\boldsymbol{\beta}} を求めたい。
ただし、対数尤度は以下のように定義する：

\begin{aligned}
l(\boldsymbol{\beta}) &= \sum_{i=1}^n l_i(\boldsymbol{\beta}),
\end{aligned}

\begin{aligned}
    l_i(\boldsymbol{\beta}) &= \log f(y_i \vert \boldsymbol{x}_i; \boldsymbol{\beta}).
\end{aligned}

 パラメータ最尤推定量の数値的な求め方対数尤度 l(\boldsymbol{\beta}) を最大化するパラメータ \hat{\boldsymbol{\beta}} を数値的に求めるには、以下のような漸化式を計算し \boldsymbol{b}^{(k)} の収束する先を求めれば良い：

\begin{aligned}
\boldsymbol{b}^{(k+1)} &= \boldsymbol{b}^{(k)} + E_Y\left[\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} \middle\vert \boldsymbol{X}; \boldsymbol{b}^{(k)}\right]^{-1} \left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta}=\boldsymbol{b}^{(k)}}. \tag{1}
\end{aligned}
ただし、 \boldsymbol{X} = (\boldsymbol{x}_1,\boldsymbol{x}_2,..., \boldsymbol{x}_n)^\top である。

 式(1)の導出\hat{\boldsymbol{\beta}} は l(\boldsymbol{\beta}) が最大となるパラメータであることから、

\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \boldsymbol{0}
となる \boldsymbol{\beta} を探したい。
\boldsymbol{\hat{\beta}} に十分近い値を持つ \boldsymbol{b} について、以下のような近似式を考える：

\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \hat{\boldsymbol{\beta}}}
\approx \left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \boldsymbol{b}} + \left.\frac{\partial^2 l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}} (\hat{\boldsymbol{\beta}} - \boldsymbol{b}).
ここで

\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \hat{\boldsymbol{\beta}}} = \boldsymbol{0}
を代入すると、上記の式は以下のように書き改めることができる：

\hat{\boldsymbol{\beta}} \approx \boldsymbol{b} - \left(\left.\frac{\partial^2 l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}}\right)^{-1}\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta}=\boldsymbol{b}}
このとき、

\left.\frac{\partial^2 l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}}
をその期待値

E_Y\left[\left.\frac{\partial^2 l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}}\middle\vert \boldsymbol{x}; \boldsymbol{b}\right]
で置き換えられるものとし、さらに

E_Y\left[\frac{\partial^2 l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top} \middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right]
= -E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right]
という関係を用いると、

\hat{\boldsymbol{\beta}}
\approx \boldsymbol{b} + E_Y\left[\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}}\middle\vert \boldsymbol{x}; \boldsymbol{b}\right]^{-1}\left.
    \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta}=\boldsymbol{b}}
のように書き換えることができる。
これをもとに、\boldsymbol{b} \to \boldsymbol{b}^{(k)}, \, \hat{\boldsymbol{\beta}} \to \boldsymbol{b}^{(k+1)} と置き換えると、式(1)が導かれる。

畳屋民也

 一般化線形モデルの最尤推定その２: 重み付き最小二乗法による別表現一般化線形モデルにおいて、対数尤度 l(\boldsymbol{\beta}) を最大化するパラメータ \boldsymbol{\beta} = \hat{\boldsymbol{\beta}} を数値的に求めるには、以下のような漸化式を計算すればよかった：

\begin{aligned}
\boldsymbol{b}^{(k+1)} &= \boldsymbol{b}^{(k)} + E_Y\left[\left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} \middle\vert \boldsymbol{X}; \boldsymbol{b}^{(k)}\right]^{-1} \left.\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta}=\boldsymbol{b}^{(k)}} \tag{1}
\end{aligned}
https://zenn.dev/link/comments/0fa9737c2e71b3
ここで、説明変数 \boldsymbol{x} とパラメータ \boldsymbol{\beta} が与えられている時の目的変数 y の（条件付き）確率密度分布関数 f(y\vert \boldsymbol{x}; \boldsymbol{\beta}) が以下のような形で表せるものとする：

\begin{aligned}
    f(y\vert \boldsymbol{x}; \boldsymbol{\beta}) = \exp\left\{ y b(\boldsymbol{x}; \boldsymbol{\beta}) + c(\boldsymbol{x}; \boldsymbol{\beta}) + d(y)\right\}.
\end{aligned}
なお、この時、l_i(\boldsymbol{\beta}) = y_i b(\boldsymbol{x}_i; \boldsymbol{\beta}) + c(\boldsymbol{x}_i; \boldsymbol{\beta}) + d(y_i) と表すことができる。
すると、式(1)は以下のように表すことができる：

\begin{aligned}
\boldsymbol{b}^{(k+1)} &=  \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)} \boldsymbol{z}^{(k)} \tag{2}
\end{aligned}
ただし、\boldsymbol{W}^{(k)}, \, \boldsymbol{z}^{(k)} は以下のように定義する：

z_i^{(k)} =  \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i^{(k)}} \right)^{-1} (y_i - \mu_i^{(k)})+ \boldsymbol{x}_i^\top \boldsymbol{b}^{(k)}

W_{ij}^{(k)} =
\begin{cases}
    \frac{1}{V_Y[y_i\vert \boldsymbol{x}_i; \boldsymbol{b}^{(k)}]} \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i^{(k)}} \right)^{-2} &\quad (i=j)\\
    0 & \quad (i \ne j)
\end{cases}
なお、

\mu_i^{(k)} = E_{Y}[y_i \vert \boldsymbol{x}_i; \boldsymbol{b}^{(k)}].
以下、これを示す。

 導出式(1)のうち、

\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}, \quad
E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top} \middle\vert \boldsymbol{X}; \boldsymbol{\beta}\right]
について式変形を行う。
なお、以降では、

\begin{aligned}
\mu_i &= \mu(\boldsymbol{x}_i; \boldsymbol{\beta})\\
&= E_{Y}[y_i \vert \boldsymbol{x}_i; \boldsymbol{\beta}]
\end{aligned}
と表すことにする。
まず、

\begin{aligned}
\frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= \left. \frac{\partial \mu}{\partial \boldsymbol{\beta}}\right|_{\mu=\mu_i} \frac{\partial l_i(\boldsymbol{\beta})}{\partial \mu}\\
&= \frac{1}{V_Y[y_i \vert \boldsymbol{x}_i; \boldsymbol{\beta}]} \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i} \right)^{-1} \boldsymbol{x}_i (y_i - \mu_i) \tag{3}
\end{aligned}
と表すことができる。
なぜなら、

\frac{\partial \mu}{\partial \boldsymbol{\beta}} = \left( \frac{dg(\mu)}{d\mu} \right)^{-1} \boldsymbol{x} \tag{4}

\begin{aligned}
\frac{\partial l_i(\boldsymbol{\beta})}{\partial \mu} &= \left.\frac{\partial b}{\partial \mu}\right|_{\mu = \mu_i} (y_i - \mu_i)\\
&= \frac{1}{V_Y[y_i \vert \boldsymbol{x}_i; \boldsymbol{\beta}]} (y_i - \mu_i)\tag{5}
\end{aligned}
が成り立つからである。
式(4)の導出g(\mu) = \boldsymbol{x}^\top \boldsymbol{\beta} を \boldsymbol{\beta} で偏微分して、

\begin{aligned}
    \frac{\partial g(\mu)}{\partial \boldsymbol{\beta}}
    &= \frac{\partial \mu}{\partial \boldsymbol{\beta}} \frac{dg(\mu)}{\partial \mu}\\
    &= \frac{\partial (\boldsymbol{\boldsymbol{x}^\top \boldsymbol{\beta}})}{\partial \boldsymbol{\beta}} = \boldsymbol{x}
\end{aligned}
したがって、

\frac{\partial \mu}{\partial \boldsymbol{\beta}} = \left( \frac{dg(\mu)}{d\mu} \right)^{-1} \boldsymbol{x}
が成り立つ。
式(5) の導出以下のような l_o(\boldsymbol{\beta}) を考える：

l_o(\boldsymbol{\beta}) = y b(\boldsymbol{x}; \boldsymbol{\beta}) + c(\boldsymbol{x}; \boldsymbol{\beta}) + d(y).
この両辺を \mu で微分すると、以下が成り立つ：

\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu} = y \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu} + \frac{\partial c(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}. \tag{6}
ここで、両辺の Y \sim f(Y\vert \boldsymbol{x}; \boldsymbol{\beta}) についての期待値をとると、以下のようになる：

E_Y\left[\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right] = \mu \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu} + \frac{\partial c(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}.
一方で

E_Y\left[\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right] = 0
より

\frac{\partial c(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu} = -\mu \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}
が成り立つので、これを式(6)に代入することで

\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu} = \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu} (y - \mu) \tag{7}
が得られる。
さらに、式(6)の分散については、以下のように表すことができる：

\begin{aligned}
    V_Y\left[\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right]
    &= \left\{\frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}\right\}^2 V_Y[Y \vert \boldsymbol{x}; \boldsymbol{\beta}]\\
    &= -E_Y\left[\frac{\partial^2 l_o(\boldsymbol{\beta})}{\partial \mu^2}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right]. \tag{8}
\end{aligned}
ここで、式(7)の両辺を \mu で偏微分することにより

\frac{\partial^2 l_o(\boldsymbol{\beta})}{\partial \mu^2} =
    \frac{\partial^2 b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu^2} (y - \mu) - \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}
となるので

\begin{aligned}
E_Y\left[\frac{\partial^2 l_o(\boldsymbol{\beta})}{\partial \mu^2}\middle\vert \boldsymbol{x}; \boldsymbol{\beta}\right]
&= - \frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}
\end{aligned}
が成り立つ。
これを式(8)に代入することで、

\left\{\frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu}\right\}^2 V_Y[y \vert \boldsymbol{x}; \boldsymbol{\beta}] = \frac{\partial b(\boldsymbol{x}, \boldsymbol{\beta})}{\partial \mu}
が成り立つことから

\frac{\partial b(\boldsymbol{x}; \boldsymbol{\beta})}{\partial \mu} = \frac{1}{V_Y[y\vert \boldsymbol{x}; \boldsymbol{\beta}]} \tag{9}
が成り立つ。
以上から、式(7)、(9)合わせて、

\frac{\partial l_o(\boldsymbol{\beta})}{\partial \mu}
= \frac{1}{V_Y[y \vert \boldsymbol{x}; \boldsymbol{\beta}]} (y - \mu)
が成り立ち、これに y=y_i, \, \boldsymbol{x}=\boldsymbol{x}_i を代入することで式(5)が示される。
したがって、

\begin{aligned}
\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}
&= \sum_{i=1}^n \frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\\
&= \sum_{i=1}^n \frac{1}{V_Y[y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}]} \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i} \right)^{-1} \boldsymbol{x}_i (y_i - \mu_i)
\end{aligned}
と表すことができる。
さらにこれを用いることで、

\begin{aligned}
    E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top} \middle\vert \boldsymbol{X}; \boldsymbol{\beta}\right]
    &= \sum_{i=1}^n E_Y\left[\frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\middle\vert \boldsymbol{x}_i; \boldsymbol{\beta}\right]\\
    &=\sum_{i=1}^n \frac{1}{V_Y[y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}]} \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i} \right)^{-2} \boldsymbol{x}_i \boldsymbol{x}_i^\top
\end{aligned}
と表すことができる。

ただし、 i \ne j のとき \partial l_i (\boldsymbol{\beta}) / \partial \boldsymbol{\beta}, \, \partial l_j (\boldsymbol{\beta}) / \partial \boldsymbol{\beta} が互いに独立であることを用いた。
ここで \boldsymbol{W}(\boldsymbol{\beta}), \, \tilde{\boldsymbol{z}} を以下のように定義する：

W_{ij}(\boldsymbol{\beta}) =
\begin{cases}
    \frac{1}{V_Y[y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}]} \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i} \right)^{-2} &\quad (i=j)\\
    0 & \quad (i \ne j)
\end{cases}

\tilde{z}_i = \left( \left. \frac{dg(\mu)}{d\mu} \right|_{\mu = \mu_i} \right)^{-1} (y_i - \mu_i)
これらを用いると、

\begin{aligned}
\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}
&= \boldsymbol{X}^\top \boldsymbol{W}(\boldsymbol{\beta}) \tilde{\boldsymbol{z}}\\
    E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top} \middle\vert \boldsymbol{X}; \boldsymbol{\beta}\right]
    &= \boldsymbol{X}^\top \boldsymbol{W}(\boldsymbol{\beta}) \boldsymbol{X}
\end{aligned}
と書き直すことができる。

したがって

\boldsymbol{z} = \tilde{\boldsymbol{z}} + \boldsymbol{X}\boldsymbol{\beta}
と置くことで

\begin{aligned}
\boldsymbol{\beta} + E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top} \middle\vert \boldsymbol{X}; \boldsymbol{\beta}\right]^{-1} \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\\
 = \left(\boldsymbol{X}^\top \boldsymbol{W}(\boldsymbol{\beta}) \boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}(\boldsymbol{\beta}) \boldsymbol{z}
\end{aligned}
が言えるので、\boldsymbol{\beta} = \boldsymbol{b}^{(k)} を代入することで式(2)が導かれる。

畳屋民也

 Logistic 回帰と最尤推定これまでの内容の例として、Logistic 回帰に当てはめて考えてみる。
https://zenn.dev/link/comments/0fa9737c2e71b3

 モデル設定Logistic 回帰の場合、y_i の確率密度分布は以下のように表すことができる：

f(y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}) = \pi(\boldsymbol{x}_i; \boldsymbol{\beta})^{y_i}\left(1 - \pi(\boldsymbol{x}_i; \boldsymbol{\beta})\right)^{1 - y_i}.
ただし、\pi(\boldsymbol{x}_i; \boldsymbol{\beta}) は y_i の（条件付き）期待値

E_Y[y_i\vert \boldsymbol{x}_i, \boldsymbol{\beta}] = \pi(\boldsymbol{x}_i; \boldsymbol{\beta})
であり、以下のように表される：

\pi(\boldsymbol{x}_i; \boldsymbol{\beta}) = \frac{1}{1 + e^{-\boldsymbol{x}_i^\top \boldsymbol{\beta}}}.
なお、上記は式変形すると

\log \left( \frac{\pi(\boldsymbol{x}_i; \boldsymbol{\beta})}{1 - \pi(\boldsymbol{x}_i; \boldsymbol{\beta})}\right) = \boldsymbol{x}_i^\top \boldsymbol{\beta}
と表せるので、

g(\mu) = \log\left(\frac{\mu}{1 - \mu}\right)
に相当する。

 最尤推定
 対数尤度とその微分ロジスティック回帰において、対数尤度 l(\boldsymbol{\beta}) は以下のように表される：

\begin{aligned}
l(\boldsymbol{\beta}) &= \sum_{i=1}^n \log f(y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}).
\end{aligned}
ここで、l(\boldsymbol{\beta}) = \sum_{i=1}^n l_i(\boldsymbol{\beta}) とすると、l_i(\boldsymbol{\beta}) = \log f(y_i; \boldsymbol{x}_i ,\boldsymbol{\beta}) は以下のように表すことができる：

l_i(\boldsymbol{\beta}) = y_i \log \left( \frac{\pi(\boldsymbol{x}_i; \boldsymbol{\beta})}{1 - \pi(\boldsymbol{x}_i; \boldsymbol{\beta})}\right) + \log \left( 1 - \pi(\boldsymbol{x}_i; \boldsymbol{\beta})\right)
l_i(\boldsymbol{\beta}) を \boldsymbol{\beta} で偏微分すると、以下のようになる：

\frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = (y_i -  \pi(\boldsymbol{x}_i; \boldsymbol{\beta}) ) \boldsymbol{x}_i.
従って、以下が成り立つ：

\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \sum_{i=1}^n(y_i -  \pi(\boldsymbol{x}_i; \boldsymbol{\beta}) ) \boldsymbol{x}_i.

\begin{aligned}
E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\middle\vert \boldsymbol{X}; \boldsymbol{\beta}\right]
&= \sum_{i=1}^n V_Y[y_i\vert \boldsymbol{x}_i; \boldsymbol{\beta}] \boldsymbol{x}_i \boldsymbol{x}_i^\top\\
&= \sum_{i=1}^n\pi(\boldsymbol{x}_i; \boldsymbol{\beta}) (1 - \pi(\boldsymbol{x}_i; \boldsymbol{\beta})) \boldsymbol{x}_i \boldsymbol{x}_i^\top\\
\end{aligned}

 最尤推定量を数値的に求める方法対数尤度 l(\boldsymbol{\beta}) を最大にする \boldsymbol{\beta} = \hat{\boldsymbol{\beta}} を数値的に求めるには、以下の更新式に従って \boldsymbol{b}^{(k+1)} を計算していく：

\begin{aligned}
    \boldsymbol{b}^{(k+1)} = \boldsymbol{b}^{(k)}+ \boldsymbol{M}(\boldsymbol{b}^{(k)})^{-1}\left. \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} \tag{1}
\end{aligned}
ただし、

\begin{aligned}
\left. \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} &= \sum_{i=1}^n(y_i -  \pi(\boldsymbol{x}_i; \boldsymbol{b}^{(k)}) ) \boldsymbol{x}_i,\\
\boldsymbol{M}(\boldsymbol{b}^{(k)}) &= \left.E_Y\left[\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^\top}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} \middle\vert \boldsymbol{X}; \boldsymbol{b}^{(k)}\right]\\
&=\sum_{i=1}^n\pi(\boldsymbol{x}_i; \boldsymbol{b}^{(k)}) (1 - \pi(\boldsymbol{x}_i; \boldsymbol{b}^{(k)})) \boldsymbol{x}_i \boldsymbol{x}_i^\top.
\end{aligned}

 式(1)の別表現https://zenn.dev/link/comments/9e95d9b4d45073
\pi_i^{(k)} = \pi(\boldsymbol{x}_i; \boldsymbol{b}^{(k)}) と表して以下のような \boldsymbol{W}^{(k)}, \, \boldsymbol{\Pi}^{(k)}, \boldsymbol{z}^{(k)} を定義する：

W_{ij}^{(k)} =
\begin{cases}
\pi_i^{(k)} \left(1 - \pi_i^{(k)} \right) &\quad (i = j)\\
0 &\quad (i \ne j)
\end{cases}

\boldsymbol{\Pi}^{(k)} = \left( \pi_1^{(k)}, \pi_2^{(k)}, ..., \pi_n^{(k)}\right)^\top

\boldsymbol{z}^{(k)} = \left(\boldsymbol{W}^{(k)}\right)^{-1} (\boldsymbol{Y} - \boldsymbol{\Pi}^{(k)}) + \boldsymbol{X}\boldsymbol{b}^{(k)}
すると、式(\ref{1}) は以下のように書き表すことができる：

\boldsymbol{b}^{(k+1)} = \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{z}^{(k)} \tag{2}
導出
\begin{aligned}
\left. \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}} &= \sum_{i=1}^n (y_i - \pi_i^{(k)}) \boldsymbol{x}_i\\
&= \boldsymbol{X}^\top \left(\boldsymbol{Y} - \boldsymbol{\Pi}^{(k)}\right)
\end{aligned}

\begin{aligned}
\boldsymbol{M}(\boldsymbol{b}^{(k)})
&=\sum_{i=1}^n\pi_i^{(k)} (1 - \pi_i^{(k)}) \boldsymbol{x}_i \boldsymbol{x}_i^\top\\
&= \boldsymbol{X}^\top \boldsymbol{W}^{(k)} \boldsymbol{X}
\end{aligned}
よって、式(1) は以下のように書き換えられる：

\begin{aligned}
    \boldsymbol{b}^{(k+1)} &= \boldsymbol{b}^{(k)}+ \boldsymbol{M}(\boldsymbol{b}^{(k)})^{-1}\left. \frac{\partial l(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right|_{\boldsymbol{\beta} = \boldsymbol{b}^{(k)}}\\
    &= \boldsymbol{b}^{(k)} + \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \left(\boldsymbol{Y} - \boldsymbol{\Pi}^{(k)}\right)\\
    &= \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)} \boldsymbol{X} \boldsymbol{b}^{(k)} + \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)}\left(\boldsymbol{W}^{(k)}\right)^{-1} \left(\boldsymbol{Y} - \boldsymbol{\Pi}^{(k)}\right)\\
    &= \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)} \left\{ \left(\boldsymbol{W}^{(k)}\right)^{-1} \left(\boldsymbol{Y} - \boldsymbol{\Pi}^{(k)}\right) + \boldsymbol{X}\boldsymbol{b}^{(k)}\right\}\\
    &= \left(\boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(k)}\boldsymbol{z}^{(k)}
\end{aligned}
なお、式(2)は \boldsymbol{z}^{(k)} を目的変数、\boldsymbol{X} を説明変数として \boldsymbol{W} でサンプルに重みをつけた重み付き最小二乗法による係数推定として捉えることもできる。