このスクラップは？

一般化線形モデルのお勉強ノート。

主に、以下の書籍で学んだ内容をまとめて行く予定：

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}$ でサンプルに重みをつけた重み付き最小二乗法による係数推定として捉えることもできる。