はじめに
どのサイトも小手先だけのロジスティック回帰分析を行っているが、理論的な数式はほとんどない。
どこよりも丁寧な"数式"でロジスティック回帰分析のパラメータ取得アルゴリズムを解説します。
※またここでは、最尤法により導出される対数尤度関数を最大にするパラメータを Newton 法を用いて求めるアルゴリズムを紹介します
ロジスティック回帰の歴史は生存確率の推定
ロジスティック回帰分析（ロジットモデル）は、その歴史から薬の薬効を調べる生存確率の推定に用いられてきました。どの薬を用いると生存率が上がり、どの薬を用いると生存率が下がるのか、そのパラメータを推定してきたのです。
そこから応用し、現在では企業の倒産確率導出でも利用されることがしばしばあります。
例えば、ある財務諸表の値ひとつひとつをリスクファクターと考えて、どの財務諸表が企業を倒産に導き、どのリスクファクターが企業を存続させるために有用なのか。そういった分析に用いられてくるようになりました。
そこで！！ここでは、専門用語が並んでもわかりづらいので企業の倒産確率を想定した数式構築を前提に進めていきます。
ロジスティック回帰分析（ロジットモデル）の構築
$K$番目のリスクファクター$X$の$i$番目の企業に関して、信用リスク度を
\tilde{Z}_i=\beta_0+\beta_1X_{i,1}+\beta_2X_{i,2}+...+\beta_kX_{i,k}+\epsilon_i
とし、倒産確率を
PD_i=F(\tilde{Z_i})=F(\beta_0+\beta_1X_{i,1}+\beta_2X_{i,2}+...+\beta_kX_{i,k}) ・・・(*)
と表す。このとき
PD_i=\frac{1}{1+exp(-\bar{Z_i})},\bar{Z_i}\equiv E[\bar{Z_i}] = \beta_0+\beta_1X_{i,1}+\beta_2X_{i,2}+...+\beta_kX_{i,k}
となる。ここで、$\bar{Z_i}$は信用リスク指標の期待値を示している。
また、リスク指標の期待値は、ロジット回帰分析の場合、対数オッズと同値である。すなわち、(＊)の逆関数より
\bar{Z_i} = \frac{PD_i}{1-PD_i}
が得られる。この式の右辺は、倒産の発生確率と倒産が発生しない確率の比、つまりオッズを対数変換したものを表している。したがって、対数オッズ比は、信用リスク度の期待値になる。また
\frac{\partial \bar{Z_i}}{\partial X_{i,k}}=\frac{\partial log(PD_i/1-PD_i)}{\partial X_{i,k}}=\beta_k
であることから、企業$i$の$k$番目のリスクファクターの対数オッズに対する感応度は，$\beta_k$ で表すことができる。対数オッズは倒産確率の 1 つの表現であり，その値が大きいほど信用リスクが高いと解釈できる。
パラメータ推定方法
最尤法：尤度関数の導出
ここでは、
PD_i=\frac{1}{1+exp(-(\beta_0+\sum_{p=1}^{K}\beta_pX_{i,p}))}
と定義する。また、このとき$\beta_0(p=1,...,K)$を$p$番目のパラメータとし、企業$i(=1,...,N)$の$p$番目のリスクファクターを$X_{i,p}$とする。また、企業$i$が倒産企業であれば1、そうでなければ0とする
Y_i = \left\{
\begin{array}{ll}
1 & (,if企業iが倒産企業) \\
0 & (,if企業iが非倒産企業)
\end{array}
\right.
を定義する。また、各企業が独立に倒産、非倒産のいずれかの状態をとるものとする。
以上のように定義するとき、パラメータ$\beta_p$に対する対数尤度関数$L$は、
L=\log \left(\prod_{i=1}^NPD_i^{Y_i} \left(1-PD_i\right)^{1-Y_i}\right)
=\sum_{i=1}^N\left(Y_i・\log\left(PD_i\right)+\left(1-Y_i\right)・\log\left( 1-PD_i\right)\right).
ここで、
\begin{align}
\log\left(1-PD_i\right)
&= \log\left( 1-\frac{1}{1+exp\left(-(\beta_0+\sum_{p=1}^K\beta_pX_{i,p})\right)} \right)
\\
&= \log\left( \frac{exp\left(-(\beta_0+\sum_{p=1}^K\beta_pX_{i,p})\right)}{1+exp\left(-(\beta_0+\sum_{p=1}^K\beta_pX_{i,p})\right)} \right)
\\
&= \log\left( exp\left(-(\beta_0+\sum_{p=1}^K\beta_pX_{i,p})\right)PD_i \right)
\\
&= \log\left( exp\left(-(\beta_0+\sum_{p=1}^K\beta_pX_{i,p})\right)\right) + \log\left(PD_i \right)
\\
&= -\beta_0-\sum_{p=1}^K \beta_pX_{i,p}+ \log(PD_i)
\end{align}
より、
\begin{align}
L
&=\sum_{i=1}^N\left(
Y_i・\log\left(PD_i\right)+\left(1-Y_i\right)・\left(-\beta_0-\sum_{p=1}^K \beta_pX_{i,p}+ \log(PD_i) \right)
\right)
\\
&=\sum_{i=1}^N\left(
\log\left(PD_i\right)+\left(1-Y_i\right)・\left(-\beta_0-\sum_{p=1}^K \beta_pX_{i,p} \right)
\right)
\\
&=\sum_{i=1}^N\left(
-\log\left(1+\exp\left(-\beta_0-\sum_{p=1}^K\beta_pX_{i,p} \right)\right)+\left(1-Y_i\right)・\left(-\beta_0-\sum_{p=1}^K \beta_pX_{i,p} \right)
\right)
\end{align}
となる。最尤法では，対数尤度が最大になるようにパラメータを決定することになる。したがって、ここでは尤度関数$L$を数値計算 (ここではニュートン法) により最適化問題を解く。
ここで、Newton法の際に利用する勾配ベクトル，Hessian 行列を以下のように定義する．
定義1：勾配ベクトル
\Delta f(\beta) = \left[
\frac{\partial f}{\partial \beta_0}
,\frac{\partial f}{\partial \beta_1}
,...
,\frac{\partial f}{\partial \beta_k}
\right].
ただし、$\beta=(\beta_0,\beta1,...,\beta_k)$とする。
定義２：Hessian行列
Hessian行列$\Delta^2f(\beta_k)$を以下のように定義する
\Delta^2f(\beta_k)=
\left(
\begin{array}{cccc}
\frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
& \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1}
& \cdots
& \frac{\partial^2 f}{\partial \beta_0 \partial \beta_K}
\\
\frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
& \frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
& \cdots & \frac{\partial^2 f}{\partial \beta_0 \partial \beta_K}
\\
\vdots & \vdots & \ddots & \vdots
\\
\frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
& \frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
& \cdots
& \frac{\partial^2 f}{\partial \beta_0 \partial \beta_0}
\end{array}
\right).
対数尤度関数の一階微分・二階微分
このとき勾配ベクトル及びHessian行列を求めるには、対数尤度関数の一階微分と２階微分が必要であることから、
\begin{align}
L
&=\sum_{i=1}^N\left(
-\log\left(1+\exp\left(-\beta_0-\sum_{p=1}^K\beta_pX_{i,p} \right)\right)+\left(1-Y_i\right)・\left(-\beta_0-\sum_{p=1}^K \beta_pX_{i,p} \right)
\right)
\end{align}
より、$j=1,...K$に対して
\begin{align}
\frac{ \partial L(\beta)}{\partial \beta_j}
&=\sum_{i=1}^N
\left[
\frac{
X_{i,j}・
\left\{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
-1
\right\}
}{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}
+
\left(
1-Y_i
\right)
・
\left(
-0-X_{i,j}
\right)
\right]
\\
&=\sum_{i=1}^N
\left[
(X_{i,j}-X_{i,j}・PD_i)-(X_{i,j}-Y・X_{i,j})
\right]
\\
&=-\sum_{i=1}^N(PD_i-Y_i)X_{i,j}
\end{align}
となる。よって、対数尤度関数の一階微分は、
\frac{ \partial L(\beta)}{\partial \beta_0}
=-\sum_{i=1}^N(PD_i-Y_i)X_{i,j},
\\
\frac{ \partial L(\beta)}{\partial \beta_j}
=-\sum_{i=1}^N(PD_i-Y_i)
となる。また、$jj=1,...,K$とし、一階微分の結果を利用することで
\begin{align}
\frac{\partial^2L(\beta)}{\partial\beta_j \partial\beta_{jj}}
&=\sum_{i=1}^N
\left\{
\frac{
0+1・X_{i,jj}・
\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}{
\left\{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
\right\}^2
}
X_{i,j}
-
0・X_{i,j}
\right\}
\\
&=\sum_{i=1}^N
\left\{
\frac{
1
}
{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}
\frac{
\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}
X_{i,j}・X_{i,jj}
\right\}
\\
&=\sum_{i=1}^N
\left\{
\frac{
1
}
{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}
\left\{
1-
\frac{
1
}{
1+\exp
\left(
-\beta_0-\sum_{p=1}^K\beta_pX_{i,p}
\right)
}
\right\}
X_{i,j}・X_{i,jj}
\right\}
\\
&=\sum_{i=1}^NPD_i・(1-PD_i)・X_{i,j}・X_{i,jj}
\end{align}
となる。よって、対数尤度関数の２階微分は、
\frac{
\partial^2L(\beta)
}{
\partial \beta_j \partial \beta_{jj}
}
= \sum_{i=1}^NPD_i・(1-PD_i)・X_
{i,j}・X_{i,jj}
となる。
Newton法
$\beta$を$N$次のベクトルとし、
y=f(\beta)
とする。ここでは、Newton法を用いて、これらが最大になる$\beta$を求める。ここで、$\beta(k)$の周りでテイラー展開より、
f(\beta)\simeq f(\beta(k))+\Delta f(\beta(k))(\beta-\beta(k))
+ \frac{1}{2}(\beta-\beta(k))^\tau\Delta^2 f(\beta(k))(\beta-\beta(k))
であり、$f$が２階連続微分可能であれば$\Delta^2f$は対象行列になる。ここで、上記の近似の精度が十分に良ければ、右辺が最大になる$\beta$を求めることで、良い近似が得られる。ここで簡易化のために
H=\Delta^2 f(\beta(k)),
\\
p=(\Delta(\beta(k)))^\tau,
\\
\xi=\beta-\beta(k)
とおき、右辺を書き直すと、
f(\beta(k))-p^\tau\xi + \frac{1}{2}\xi^\tau H\xi
=
f(\beta(k))-\frac{1}{2} p^\tau H^{-1}p
+\frac{1}{2}(\xi + H^{-1}p)^\tau H(\xi + H^{-1}p)
であり、
このHessian行列は最尤法の対数尤度関数に対するものであり、このときは、つねに判不定値であることに注意すると、
\tau = -H^{-1}p
のときに、最大値をとることがわかる。
したがって、改めて$\tau + \beta(k)$を$\beta(k+1)$とおき、もとの記号をつかって書き表すと、以下の漸化式が得られる。
\beta(k+1) = \beta(k)
-
\left(
\Delta^2 f(\beta(k))
\right)^{-1}
(\Delta f(\beta(k)))^\tau
ロジスティック回帰分析のアルゴリズム
以上の最尤法とNewton法より、以下のようなアルゴリズムにより、パラメータ推定を行う。
１. $\bar{Z_i}$を適当に定め、$k=0$とする
２. $PD_i$を計算して求める。
３. 最尤法により導出した
\frac{ \partial L(\beta)}{\partial \beta_j}
=-\sum_{i=1}^N(PD_i-Y_i)
,
\\
\frac{
\partial^2L(\beta)
}{
\partial \beta_j \partial \beta_{jj}
}
= \sum_{i=1}^NPD_i・(1-PD_i)・X_
{i,j}・X_{i,jj}
,
\\
for:j,jj=0,...,K and X_{i,0}=1
の右辺をそれぞれ求め、$\Delta f(\beta))$、$\Delta^2 f(\beta(k))$を作成する。
４. $\left( \Delta^2 f(\beta(k)) \right)^{-1}$を計算し、$\left( \Delta^2 f(\beta(k)) \right)^{-1}(\Delta f(\beta(k)))^\tau$を求め、
\beta(k+1) = \beta(k)
-
\left(
\Delta^2 f(\beta(k))
\right)^{-1}
(\Delta f(\beta(k)))^\tau
により、解$\beta(k)$を更新する。
５. $k=k+1$とし、更新された$\beta(k+1)$を用いて$\bar{Z_i}$を更新し、対数尤度関数の更新を行う。
６. 対数尤度関数の更新にはほとんど変化がないのであれば、ループを終了させる。変化があれば2にもどる。