完全に自分用メモ。
モチベーション
自分のバイタルデータ（体重、身長、性別、...）を用いて、来年病気にかかるかどうかを推定したい。
ベルヌーイ分布
もっと一般化すると、
あるパラメータ$\boldsymbol{x}$を用いて、ある事象$Y$が起こる($Y=1$)か、起こらない($Y=0$)かを推定したい、ということになる。
2種類のみの結果しか得られない今回の場合のような試行はベルヌーイ試行と呼ばれ、$Y=1$のときの確率$P(Y=1)$を$p$とすると、$P(Y=0)$のときの確率は$1-p$となり、期待値は$p$、分散は$p(1-p)$となる。
また、確率質量関数は、
P(Y=k) = p^k(1-p)^{1-k} \qquad for\qquad k ∈ {0,1}
となる。
尤度関数
ある患者のバイタルデータを$\boldsymbol{x}$とし、来年病気にかかるという結果を$Y=1$、病気にかからないという結果を$Y=0$とすると、この患者さんのデータが観測されたときの、病気の有無$Y$に対するベルヌーイの確率質量関数は、上で示した通り、
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)$である。
L(p|Y=1)=p^1(1-p)^0=p
逆に病気にかからなかった$(Y=0)$とすると、
L(p|Y=0)=p^0(1-p)^1=1-p
となる。
こういったデータが多く取得している場合を考える。
$n$人の患者さんのデータ$\boldsymbol{x}_1$~$\boldsymbol{x}_n$とそれぞれの患者さんの病気にかかったかどうかのデータがあるとする。それぞれの患者は独立だと仮定すると、全体の尤度関数は各患者の個別の尤度関数の積になる。
L(\boldsymbol{p}|\boldsymbol{Y})=\prod^n_{i=1}p_i^{Y_i}(1-p_i)^{1-Y_i}
この尤度関数を最大にするようなモデルのパラメータを見つけるのが最尤推定である。
※最尤推定はこの動画がわかりやすい。
ロジスティック回帰
$p$が
\ln(\frac{p}{1-p})=\beta_0+\Sigma\beta_k\boldsymbol{x}_k
というモデルで表せると仮定したとき、このモデルをロジスティック回帰モデルという。
これを変形すると、
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$とした。
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_i\ln(p_i)+(1-Y_i)\ln(1-p_i)]
\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)})]
式変形すると、
\ln(L(\boldsymbol{p}|\boldsymbol{Y}))=\Sigma^n_{i=1}[Y_iz_i-\ln(1+\exp(z_i)]
となる。
$x_0=1$として、$z_i$をまとめると、
\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$の解を求めるニュートン・ラプソン法の更新式は以下の通りである。
\beta_{k+1}=\beta_k-\frac{f(\beta_k)}{f'(\beta_k)}
多次元展開し、ベクトル表記すると
\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回偏微分
\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$は以下のようになる
\nabla L(\boldsymbol\beta)=\boldsymbol X^T(\boldsymbol Y - \boldsymbol P)
2回偏微分
\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改変微分に対応する。行列形式で表現すると、
\boldsymbol H = - \boldsymbol X^T \boldsymbol W \boldsymbol X
ここで、$\boldsymbol W$はn×nの対角行列で、対角要素は$w_{ii}=p_i(1-p_i)$である。
これで、ニュートン・ラプソン法の計算ができるようになる。
反復過程
- 初期値の設定
まず$\boldsymbol\beta^{(0)}$を設定します。例えば、すべて0。 - 収束するまで以下のステップを繰り返します。
- 現在の$\boldsymbol\beta^{(k)}$を用いて、各患者の$p^{(k)}_i$を計算する。
- $p^{(k)}_i$を用いて、勾配ベクトル$\nabla L$を計算する
- $p^{(k)}_i$を用いて、重み行列$\boldsymbol W^{(k)}$を計算する。
- $\boldsymbol W^{(k)}$を用いて、ヘッセ行列$\boldsymbol H^{(k)}$を計算する。
- パラメータを更新する。
||$\boldsymbol\beta^{(k+1)}-\boldsymbol\beta^{(k)}||$が非常に小さくなった場合に計算を終了する。\boldsymbol\beta^{(k+1)} = \boldsymbol\beta^{(k)} - [\boldsymbol H^{(k)}]^{-1} \nabla L
これで、ロジスティック回帰モデルの係数が求まる。