ロジスティック回帰分析！

説明変数K個・目的変数1個のロジスティック回帰分析のアルゴリズムについて、プログラムとして実装できるようにするために自分なりに理解してまとめたものです。

ロジスティック回帰分析の実装について調べると、ほとんどが基本的な内容とライブラリを使った実装方法の説明であり、ライブラリに頼らずに自分で実装する方法について知りたいのに、それを詳しく説明しているサイトがなかったので。

0. 前提条件

• 回帰分析（線型回帰、曲線回帰、重回帰等）についての基本的な知識があること。

1. アルゴリズム

ロジスティック回帰分析では、求める回帰式として次の式を想定する。
　（目的変数\(1\)個、説明変数\(K\)個、データ数\(N\)個を想定）

\[p = \frac{1}{1 + e^{-(\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_Kx_K)}} \tag{1}\]

　(ここの \(p\) は確率（目的変数; \(y\) でもよいが確率なので \(p\) とした）、 \(e\) はネイピア数)

また、「シグモイド関数」を次の式で定義すると、

\[\sigma(x) = \frac{1}{1 + e^{(-x)}} \tag{2}\]

式(1)は次のように表される。

\[p = \sigma(\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_Kx_K) \tag{3}\]

さて、式(1)を変形すると、

\[\frac{p}{1 - p} = e^{\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_Kx_K} \tag{4}\]

　（ここの \(\frac{p}{1 - p}\) は「オッズ(odds)」と呼ばれるもの）

対数をとると、

\[\log\frac{p}{1 - p} = \beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_Kx_K \tag{5}\]

　（ここの \(\log\frac{p}{1 - p}\) は「ロジット(logit)」と呼ばれるもの）

また、式(5)（もしくは元の式(1)）を変形し、 \(i\) 番目のデータについて表現すると、
　（また、 \(e\) で表現すると指数部分が小さくなるので、以降 \(\exp\) で表現する）

\[p_i = \frac{\exp(\beta_0+\beta_1 x_{i1}+\beta_2 x_{i2}+\cdots+\beta_Kx_{iK})}{1 + \exp(\beta_0+\beta_1 x_{i1}+\beta_2 x_{i2}+\cdots+\beta_Kx_{iK})} \tag{6}\]

以下、「最尤法」を用いて係数を推定する。

まず、\(i\)番目のデータについて、予測値との当てはまりを考えてみる。

\[P_i = y_i^{t_i}(1 - y_i)^{1 - t_i} \tag{7}\]

（\(y_i\)：\(i\)番目のデータにおける、ロジスティック回帰式から算出される確率（式(1)の\(p\)）、
　\(t_i\)：\(i\)番目のデータにおける実際の値（\(0\) or \(1\)））

\(P_i\)が大きい方が予測式の当てはまりが良いとすると、
　　・\(t_i = 1\)の時：\(P_i = y_i\)なので、\(y_i\)は大きい方がいい
　　・\(t_i = 0\)の時：\(P_i = 1 - y_i\)なので、\(y_i\)は小さい方がいい
ということになる。

この式(7)を全テータに対して適用すると、

\[L(\beta) = \prod_{i=1}^{N}y_i^{t_i}(1 - y_i)^{1 - t_i} \tag{8}\]

この\(L(\beta)\)を「尤度関数」と呼び、この式が最大になるように\(\beta\)の値を定めれば、当てはまりが最も良くなると言える。

式(8)の最大化が目的だが、乗算は計算が複雑な上に、確率の総乗を計算しようとすると、値が限りなく\(0\)に近付いてしまう。（アンダーフロー）　よって、この式(8)の対数をとることによって対処する。

\[-\log L(\beta) = -\sum_{i=1}^{N}\{t_i\log y_i + (1 - t_i)\log(1 - y_i)\} \tag{9}\]

対数をとったこの負の式(9)の最小化で元の式(8)の最大化となる。

式(9)の最小値は「勾配降下法」により求める。

まず、次の式を「交差エントロピー誤差関数（\(\subset{損失関数}\)）」とよぶ。（式(9)を\(E(\beta)\)と置いただけ）

\[E(\beta) = -\log L(\beta) = -\sum_{i=1}^{N}\{t_i\log y_i + (1 - t_i)\log(1 - y_i)\} \tag{10}\]

この関数をパラメータ\(\beta_i\)について偏微分する。（計算過程省略）

（\(\beta_0\)はバイアス（定数）なので、\(k=0\)と\(k\in\{1,2,\cdots,K\}\)（今回説明変数はK個を想定）で分ける）

\[\frac{\partial E(\beta)}{\partial\beta_0} = \sum_{i=1}^{N}(y_i - t_i)\ \ (k=0の場合) \tag{11}\] \[\frac{\partial E(\beta)}{\partial\beta_k} = \sum_{i=1}^{N}x_{ik}(y_i - t_i)\ \ (k\in\{1,2,\cdots,K\}の場合) \tag{12}\]

この式(11),(12)が\(0\)になるようにパラメータ\(\beta\)を算出すればよいが、解析的には算出できないので、勾配降下法を用いる。

以下のようにして、パラメータ\(\beta_k\)を更新していく。\(\alpha\)は（機会学習で言うところの）学習率で、\(0.01\)等の数値を使用する。要は、\(\beta_k\)で偏微分した値に学習率を乗じて元の\(\beta_k\)から減算していき、値がほとんど変化しなくなったら、値の更新を終了させる、ということ。

\[\beta_k \leftarrow \beta_k - \alpha \frac{\partial E(\beta)}{\partial\beta_k} \tag{13}\]

ちなみに、データ数が多いと式(11),(12)の値が大きくなってしまうので、式(11),(12)を次のようにすることもあるようだ。

\[\frac{\partial E(\beta)}{\partial\beta_0} = \frac{1}{N}\sum_{i=1}^{N}(y_i - t_i)\ \ (k=0の場合) \tag{14}\] \[\frac{\partial E(\beta)}{\partial\beta_k} = \frac{1}{N}\sum_{i=1}^{N}x_{ik}(y_i - t_i)\ \ (k\in\{1,2,\cdots,K\}の場合) \tag{15}\]

---

近いうちに、プログラムに実装してみたいです。

以上。