【自分用メモ】ロジスティック回帰モデルの係数の求め方

完全に自分用メモ。

# モチベーション

自分のバイタルデータ（体重、身長、性別、...）を用いて、来年病気にかかるかどうかを推定したい。


# ベルヌーイ分布
もっと一般化すると、
あるパラメータ$\boldsymbol{x}$を用いて、ある事象$Y$が起こる($Y=1$)か、起こらない($Y=0$)かを推定したい、ということになる。

2種類のみの結果しか得られない今回の場合のような試行はベルヌーイ試行と呼ばれ、$Y=1$のときの確率$P(Y=1)$を$p$とすると、$P(Y=0)$のときの確率は$1-p$となり、期待値は$p$、分散は$p(1-p)$となる。

また、確率質量関数は、
```math
P(Y=k) = p^k(1-p)^{1-k} \qquad for\qquad k ∈ {0,1}
```
となる。

# 尤度関数
ある患者のバイタルデータを$\boldsymbol{x}$とし、来年病気にかかるという結果を$Y=1$、病気にかからないという結果を$Y=0$とすると、この患者さんのデータが観測されたときの、病気の有無$Y$に対するベルヌーイの確率質量関数は、上で示した通り、
```math
P(Y=k) = p(\boldsymbol{x})^k(1-p(\boldsymbol{x}))^{1-k} \qquad for\qquad k ∈ {0,1}
```
となる。すなわち、$p$は$x$から計算される病気にかかる確率を表している。
もしこの患者さんが病気にかかる（$Y=1$）と観測されたとする。この「観測データ」が与えられたときに、モデルのパラメータである$p$がどれくらい尤もらしいかを示すのが尤度関数$L(p|Y=1)$である。
```math
L(p|Y=1)=p^1(1-p)^0=p
```
逆に病気にかからなかった$(Y=0)$とすると、
```math
L(p|Y=0)=p^0(1-p)^1=1-p
```
となる。
こういったデータが多く取得している場合を考える。
$n$人の患者さんのデータ$\boldsymbol{x}_1$~$\boldsymbol{x}_n$とそれぞれの患者さんの病気にかかったかどうかのデータがあるとする。それぞれの患者は独立だと仮定すると、全体の尤度関数は各患者の個別の尤度関数の積になる。
```math
L(\boldsymbol{p}|\boldsymbol{Y})=\prod^n_{i=1}p_i^{Y_i}(1-p_i)^{1-Y_i}
```
この尤度関数を最大にするようなモデルのパラメータを見つけるのが最尤推定である。

※最尤推定は[この動画](https://www.youtube.com/watch?v=iSCKAQCLiUg&list=PLUB2-8oQULNgVWS1C4PfnJwXnfwUthXsz)がわかりやすい。

# ロジスティック回帰
$p$が
```math
\ln(\frac{p}{1-p})=\beta_0+\Sigma\beta_k\boldsymbol{x}_k
```
というモデルで表せると仮定したとき、このモデルをロジスティック回帰モデルという。
これを変形すると、
```math
p=\frac{1}{1+\exp(-(\beta_0+\Sigma\beta_k\boldsymbol{x}_k))}
```
となる。
私達の目標は、パラメータ$\beta_0$、$\beta_1$、...が尤もらしいものを見つけることである。

## 損失関数
ここで、尤度関数$L(\boldsymbol{p}|\boldsymbol{Y})$を思い出す。この対数$\ln(L(\boldsymbol{p}|\boldsymbol{Y}))$ともとの尤度関数の最大値をとるパラメータは同じである。計算を簡単にするため、対数を取り、$p=\frac{1}{1+\exp(-(z_i))}$を代入する。$z_i=\beta^i_0+\Sigma_k\beta^i_k\boldsymbol{x}^i_k$とした。

```math
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_i\ln(p_i)+(1-Y_i)\ln(1-p_i)]
```

```math
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_i\ln(\frac{1}{1+\exp(-z_i)})+(1-Y_i)\ln(1-\frac{1}{1+\exp(-z_i)})]
```
式変形すると、
```math
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_iz_i-\ln(1+\exp(z_i)]
```
となる。
$x_0=1$として、$z_i$をまとめると、
```math
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_i\boldsymbol{x}^T_i\boldsymbol\beta-\ln(1+\exp(\boldsymbol{x}^T\boldsymbol\beta)]
```
となり、これは損失関数と呼ばれる。
損失関数を$\boldsymbol\beta$で偏微分したときにゼロとなる点が、尤もらしいパラメータになる。この偏微分=0の解は解析的に解くことが難しいので、ニュートン・ラプソン法を用いる。

## ニュートン・ラプソン法
一般に、関数$f(\beta)=0$の解を求めるニュートン・ラプソン法の更新式は以下の通りである。
```math
\beta_{k+1}=\beta_k-\frac{f(\beta_k)}{f'(\beta_k)}
```
多次元展開し、ベクトル表記すると
```math
\boldsymbol\beta_{k+1}=\boldsymbol\beta_k-\frac{f(\boldsymbol\beta_k)}{f'(\boldsymbol\beta_k)}
```
$f(\beta)$は損失関数を偏微分したものであり、
$f(\beta_j)=\frac{\partial}{\partial\beta_j}\ln L$
となる。なので、更新式を作成するにあたって、損失関数の1回偏微分と2回偏微分した式が必要となる。
## 1回偏微分
```math
\displaylines{
\frac{\partial}{\partial\beta_j}\ln L = \Sigma^n_{i=1}[Y_ix_{ij}-\frac{\exp(x^T_i \beta_j)}{1+\exp(x^T_i \beta_j)}x_{ij}] \\
 = \Sigma^n_{i=1}x_{ij}[Y_i-p_i]
}
```
※$p_i=\frac{\exp(x^T_i \beta_j)}{1+\exp(x^T_i \beta_j)}$を使用している。

ベクトル形式で表現すると、勾配ベクトル$\nabla L$は以下のようになる
```math
\nabla L(\boldsymbol\beta)=\boldsymbol X^T(\boldsymbol Y - \boldsymbol P)
```
## 2回偏微分
```math
\frac{\partial^2L}{\partial \beta_j \partial \beta_k} = \Sigma^n_{i=1} -x_{ij} x_{ik} p_i(1-p_i)
```
$\frac{\partial p_i}{\partial \beta_k}=x_{ik} p_i(1-p_i)$となることを用いている。
ヘッセ行列$\boldsymbol H$の各要素$\boldsymbol H_{jk}$はこの2改変微分に対応する。行列形式で表現すると、
```math
\boldsymbol H = - \boldsymbol X^T \boldsymbol W \boldsymbol X
```
ここで、$\boldsymbol W$はn×nの対角行列で、対角要素は$w_{ii}=p_i(1-p_i)$である。

これで、ニュートン・ラプソン法の計算ができるようになる。

## 反復過程
1. 初期値の設定
まず$\boldsymbol\beta^{(0)}$を設定します。例えば、すべて0。
1. 収束するまで以下のステップを繰り返します。
    1. 現在の$\boldsymbol\beta^{(k)}$を用いて、各患者の$p^{(k)}_i$を計算する。
    1. $p^{(k)}_i$を用いて、勾配ベクトル$\nabla L$を計算する
    1. $p^{(k)}_i$を用いて、重み行列$\boldsymbol W^{(k)}$を計算する。
    1. $\boldsymbol W^{(k)}$を用いて、ヘッセ行列$\boldsymbol H^{(k)}$を計算する。
    1. パラメータを更新する。
    ```math
    \boldsymbol\beta^{(k+1)} = \boldsymbol\beta^{(k)} - [\boldsymbol H^{(k)}]^{-1} \nabla L
    ```
    ||$\boldsymbol\beta^{(k+1)}-\boldsymbol\beta^{(k)}||$が非常に小さくなった場合に計算を終了する。
    

これで、ロジスティック回帰モデルの係数が求まる。