回帰分析の理解まとめ〜ロジスティック回帰分析〜

これまで回帰分析のうち単回帰分析と重回帰分析についてまとめてきました。

takarotoooooo.hatenablog.com
takarotoooooo.hatenablog.com

今回はロジスティック回帰に関してまとめます。

ロジスティック回帰分析とは？

ロジスティック回帰は、事象が発生するか、発生しないかのように二値をとるような問題（コインを投げて表が出るか裏が出るかのような問題です）に関して、発生する（または発生しない）確率を予測するための分析手法です。

ロジスティック回帰ではシグモイド関数というものを利用します。

シグモイド関数

シグモイド関数は下記の式で表される関数で、下記のようなグラフを取ります。

$\begin{eqnarray} \sigma(x) = \frac{1}{1+e^{-x}} \end{eqnarray}$

ご覧の通り、シグモイド関数は、 $x$ の値が大きければ1に近づき、小さければ0に近づくような関数です。

このシグモイド関数に対して、重回帰分析の時と同様に $M$ 個の説明変数 $x_1, x_2, ... , x_M$ とパラメータ $w_1, w_2, ... , w_M$ を用いて、事象の発生確率の予測値 $\hat{y}$ を求めます。
$\begin{eqnarray} \hat{p}(x_1, x_2, ... , x_M) &=& \sigma(w_{0}x_{0} + w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M}) \\ &=& \frac{1}{1+e^{-(w_{0}x_{0} + w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M} )}} \end{eqnarray}$

なお、ここで $w_{0}x_{0}$ をバイアス項とし、
$\boldsymbol{x} = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_M \end{pmatrix}, \boldsymbol{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_M \end{pmatrix}$
とすると
$\begin{eqnarray} \hat{p}(\boldsymbol{x}) &=& \sigma( \boldsymbol{x}^{T}\boldsymbol{w} ) &=& \frac{1}{1+e^{-( \boldsymbol{x}^{T}\boldsymbol{w} )}} \end{eqnarray}$
と表すことができます。

評価関数

単回帰分析、重回帰分析同様に最適なパラメータ $\boldsymbol{w}$ を求めるために評価関数を用意します。

ロジスティック回帰分析の場合、観測される目的変数 $y$ の値は事象が発生した場合=1, 発生しなかった場合=0のように二値になります。
しかしながらロジスティック回帰分析では、発生する（ $y = 1$ となる）確率を予測するため、二乗誤差がそのままでは利用することができません。

そこで $N$ を観測データの個数として、下記のような式を考えます。
$\begin{eqnarray} \mathcal{L( \boldsymbol{w} )} = \sum_{i=1}^{N}{ -y_{i}\log{ (\hat{p}(\boldsymbol{x}_i)) } - (1 - y_{i})\log{ (1 - \hat{p}(\boldsymbol{x}_i)) } } \end{eqnarray}$

この式のうち前の項である下記は、事象が発生する（ $y = 1$ である）場合の誤差を表しています。
$\begin{eqnarray} {-}y_{i}\log{ (\hat{p}(\boldsymbol{x}_i)) } \end{eqnarray}$
下記のグラフからわかるように、確率の予測値 $\hat{p}(\boldsymbol{x}_i)$ が小さい（= 発生する確率が低いと予測した）場合は誤差が大きくなり、大きい（= 発生する確率が高いと予測した）場合は誤差が小さくなります

一方後の項である下記は、事象が発生しない（ $y = 0$ である）場合の誤差を表しています。
$\begin{eqnarray} {-}(1 - y_{i}) \log{ (1 - \hat{p}(\boldsymbol{x}_i)) } \end{eqnarray}$
こちらの項は逆に、確率の予測値が小さい場合は誤差も小さくなり、大きい場合は誤差も大きくなります。

つまり、事象が発生した目的変数に対する説明変数 $\boldsymbol{x}$ による予測確率 $\hat{p}(\boldsymbol{x})$ はより大きく、逆に事象が発生しなかった目的変数に対する説明変数による予測確率はより小さくなれば、評価関数 $\mathcal{L( \boldsymbol{w} )}$ が最小になるので、それを満たすパラメータ $\boldsymbol{w}$ を求めます。

パラメータを求めるための準備

最適なパラメータ $\boldsymbol{w}$ を求めるために上記の $\mathcal{L}$ の式をまずは変形します。

$\begin{eqnarray} \mathcal{L} &=& \sum_{i=1}^{N}{ -y_{i}\log{\hat{p}(\boldsymbol{x}_i)} - (1-y_{i})\log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}\log{\hat{p}(\boldsymbol{x}_i)} + \log{(1 - \hat{p}(\boldsymbol{x}_i))} - y_{i}\log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}\{ \log{\hat{p}(\boldsymbol{x}_i)} - \log{(1 - \hat{p}(\boldsymbol{x}_i))} \} + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}(\log{\frac{\hat{p}(\boldsymbol{x}_i)}{(1 - \hat{p}(\boldsymbol{x}_i))}}) + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \end{eqnarray}$

ここで下記二つの式を利用することで
$\begin{eqnarray} \hat{p}(\boldsymbol{x}_i) &=& \frac{1}{1+e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )}} \\ (1 - \hat{p}(\boldsymbol{x}_i)) &=& \frac{ e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )} }{1+e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )}} \end{eqnarray}$

下記が成り立ちます。
$\begin{eqnarray} \frac{ \hat{p}(\boldsymbol{x}_i) }{ 1 - \hat{p}(\boldsymbol{x}_i) } &=& \frac{ 1 }{ e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )} } &=& e^{ \boldsymbol{x}_{i}^{T}\boldsymbol{w} } \end{eqnarray}$
この値をオッズと呼び、両辺対数をとった下記をロジットと呼びます。

$\begin{eqnarray} \log{ \frac{ \hat{p}(\boldsymbol{x}_i) }{ 1 - \hat{p}(\boldsymbol{x}_i) } } &=& \boldsymbol{x}_{i}^{T}\boldsymbol{w} \end{eqnarray}$

このロジットを利用してさらに式を変形します。
$\begin{eqnarray} \mathcal{L} &=& -\sum_{i=1}^{N}{ y_{i}(\log{\frac{\hat{p}(\boldsymbol{x}_i)}{(1 - \hat{p}(\boldsymbol{x}_i))}}) + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) + \log{(1 - \frac{1}{ 1 + e^{(- \boldsymbol{x}_{i}^{T} \boldsymbol{w} )}})}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) + \log{(\frac{e^{( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} )}}{1 + e^{( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} )}})}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) } + \log{(e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) } + ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} \\ &=& -\sum_{i=1}^{N}{ (y_{i} - 1) ( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })}} \end{eqnarray}$

次に下記の式を考えていきます
$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}}{\mathcal{L(\boldsymbol{w})}} \end{eqnarray}$

勾配を求める

$j$ 番目のパラメータ $w_j$ における偏微分を考えます。

$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}_{j}}{\mathcal{L(\boldsymbol{w})}} &=& \frac{\partial}{\partial\boldsymbol{w}_{j}}{(-\sum_{i=1}^{N}{ (y_{i} - 1) ( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })}})} \\ &=& -\sum_{i=1}^{N}({ (y_i - 1)\boldsymbol{x}_{ij} - \frac{e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}{1 + e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}\boldsymbol{x}_{ij} }) \\ &=& -\sum_{i=1}^{N}({ (y_i - \frac{1}{1 + e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}){\boldsymbol{x}_{ij}} }) \\ &=& -\sum_{i=1}^{N}({ (y_i - \sigma({\boldsymbol{x}_{i}^{T} \boldsymbol{w}})){\boldsymbol{x}_{ij}} }) \\ &=& \sum_{i=1}^{N}({ (\sigma({\boldsymbol{x}_{i}^{T} \boldsymbol{w}}) - y_i){\boldsymbol{x}_{ij}} }) \\ &=& \begin{pmatrix} {\boldsymbol{x}_{1j}} & {\boldsymbol{x}_{2j}} & \dots & {\boldsymbol{x}_{Nj}} \end{pmatrix} \begin{pmatrix} \sigma({\boldsymbol{x}_{1}^{T} \boldsymbol{w}}) - y_1 \\ \sigma({\boldsymbol{x}_{2}^{T} \boldsymbol{w}}) - y_2 \\ \vdots \\ \sigma({\boldsymbol{x}_{N}^{T} \boldsymbol{w}}) - y_N \end{pmatrix} \\ &=& \begin{pmatrix} {\boldsymbol{x}_{1j}} & {\boldsymbol{x}_{2j}} & \dots & {\boldsymbol{x}_{Nj}} \end{pmatrix} \begin{pmatrix} \sigma({\boldsymbol{x}_{1}^{T} \boldsymbol{w}}) \\ \sigma({\boldsymbol{x}_{2}^{T} \boldsymbol{w}}) \\ \vdots \\ \sigma({\boldsymbol{x}_{N}^{T} \boldsymbol{w}}) \end{pmatrix} - \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \end{eqnarray}$

$\begin{eqnarray} \boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}_{1}^{T} \\ \boldsymbol{x}_{2}^{T} \\ \vdots \\ \boldsymbol{x}_{N}^{T} \end{pmatrix}, \quad \boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \end{eqnarray}, \quad \boldsymbol{x}^{(j)} = \begin{pmatrix} \boldsymbol{x}_{1j} \\ \boldsymbol{x}_{2j} \\ \vdots \\ \boldsymbol{x}_{Nj} \\ \end{pmatrix}$
とすることで、下記のように表すことができます。
$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}_{j}}{\mathcal{L(\boldsymbol{w})}} &=& {\boldsymbol{x}^{(j)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \end{eqnarray}$

これが勾配の $j$ 番目の要素となるので、勾配をまとめると

$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}}{\mathcal{L(\boldsymbol{w})}} &=& \begin{pmatrix} {\boldsymbol{x}^{(0)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ {\boldsymbol{x}^{(1)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ {\boldsymbol{x}^{(2)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ \vdots \\ {\boldsymbol{x}^{(M)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ \end{pmatrix} \\ &=& \begin{pmatrix} {\boldsymbol{x}^{(0)}} & {\boldsymbol{x}^{(1)}} & {\boldsymbol{x}^{(2)}} & \dots & {\boldsymbol{x}^{(M)}} \end{pmatrix}^T (\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ &=& \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \end{eqnarray}$

パラメータを求める

評価関数が最小値をとるパラメータ $\boldsymbol{w}$ は上記の勾配が $\boldsymbol{0}$ となるところになります。
勾配は傾きが最も大きくなる方向を示すベクトルなので、勾配が示す方向に $\boldsymbol{w}$ を変化させていけば、勾配が $\boldsymbol{0}$ の極小値にたどり着くことができます。
勾配を利用して最適なパラメータを求める方式の一つに最急降下法があるので、ここではこの方法で説明します。

最急降下法

パラメータの初期値 $\boldsymbol{w}_0$ を決定した上で、次の手順を繰り返すことで最適なパラメータを求めるのが最急降下法です。

現在のパラメータ $\boldsymbol{w}_{k}$ での勾配を求める
求められた勾配の反対方向にあらかじめ決めた学習率 $\eta$ 分、パラメータを移動させ新たなパラメータ $\boldsymbol{w}_{k+1}$ とする

つまり、下記の式を繰り返し計算することで、最適なパラメータが求められます。
$\begin{eqnarray} \boldsymbol{w}_{k+1} = \boldsymbol{w}_{k} - \eta \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w_k}) - \boldsymbol{y}) \end{eqnarray}$

上記の手順を繰り返してパラメータを求めたイメージが下記になります。

まとめ

ロジスティック回帰分析はある事象が発生する確率を予測する分析手法
シグモイド関数を利用する
最適なパラメータは $\boldsymbol{w}_{k+1} = \boldsymbol{w}_{k} - \eta \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w_k}) - \boldsymbol{y})$ を繰り返し解くことで求めることができる

今回はロジスティック回帰分析に関して概要からパラメータの求め方までをまとめてみました。
具体的なデータを使っての求め方は長くなってしまったので、別でまとめてようと思います。