【JDLA E資格】ロジスティック回帰

#はじめに
JDLA E資格試験で出題されるロジスティック回帰について解説した記事です。
E資格試験の機械学習パートでは、ロジスティック回帰のモデルやオッズ、学習方法(順伝播と逆伝播)などが出題されます。
また、E資格試験では、ロジスティック回帰は実装問題含めて頻出ですので、本稿の内容を理解しておくのは必須です。

なお、他パートの具体的な解説については、下記をご覧ください。
[E資格試験に関する私の投稿記事リスト][link-1]
[link-1]:https://qiita.com/fridericusgauss/items/5a97f2645cdcefe15ce0

###目次

1. 機械学習の分類
2. ロジスティック回帰
3. ロジスティック回帰の学習
4. おわりに

#機械学習の分類
E資格試験で登場する機械学習モデル・アルゴリズムは下記表のように分類されます。

モデル	学習タイプ	タスク
線形回帰	教師あり	回帰
ロジスティック回帰	教師あり	分類
SVM	教師あり	分類
k-近傍法	教師あり	分類
k-means	教師なし	クラスタリング
k-means++	教師なし	クラスタリング
主成分分析	教師なし	次元削減
動的計画法	強化学習	ー
モンテカルロ法	強化学習	ー
Q-Learning	強化学習	ー
SARSA	強化学習	ー

#ロジスティック回帰
###概要
ロジスティック回帰は、教師あり学習で2値分類タスクのための(説明変数を用いて2値の目的変数を表現する)の数理モデルです。
説明変数と目的変数の組のデータを学習させ、目的変数の値(クラス)を予測するモデルを構築できます。
__ロジスティック回帰は、一般化線形モデルの一種__です。

###ロジスティック回帰のモデル
目的変数$y$、説明変数$\boldsymbol{x}=(x_1,\cdots,x_M)^{\mathrm{T}}$としたとき、ロジスティック回帰モデルは式(1)で表されます。

y=\sigma(\boldsymbol{w}^{\mathrm{T}}\boldsymbol{x}+w_{0})
\tag{1}

$\boldsymbol{w}=(w_1,\cdots,w_M)^{\mathrm{T}},w_{0}$はモデルパラメータで、データを学習することで決定します。

式(1)は、ダミー変数を用いて式(2)に変形できます。

y=\sigma(\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})
\tag{2}

\begin{align}
\boldsymbol{w}^{\prime}=&(w_1 ,\cdots,w_M ,w_0 )^{\mathrm{T}}\\
\boldsymbol{x}^{\prime}=&(x_1 ,\cdots,x_M ,1)^{\mathrm{T}}
\end{align}

活性化関数$\sigma:\mathbb{R}\rightarrow [0,1]$はシグモイド関数で、式(3)で表されます。

\sigma(z)=\frac{1}{1+\exp{(-z)}}
\tag{3}

シグモイド関数によって、線形結合した出力を非線形変換し、$[0,1]$に収めます。
シグモイド関数は、$z=0$が変曲点で、$z<0$の領域で値が0に近くなり、$z>0$の領域で値が1に近くなるので、2値(2クラス)の出力に対応しています。

###確率的解釈とオッズ
天下り的な説明ですが、過去の記事を理解していれば覚えるのは苦労はないと思います。

ロジスティック回帰モデル（式(1)）による目的変数の出力$\hat{y}$は、条件付き確率$P
(y=1|\boldsymbol{x})$に対応しています。
よって、式(4)のような関係で表せます。

\begin{align}
P(y=1|\boldsymbol{x})=&\hat{y}\\
P(y=0|\boldsymbol{x})=&1-\hat{y}
\end{align}
\tag{4}

また、$\boldsymbol{x}$が与えられた条件で試行した結果、確率変数が$y\in \{0,1\}$、その条件付き確率が式(4)として与えられるとき、オッズ$o$は式(5)で表されます。

o=\frac{P(y=1|\boldsymbol{x})}{P(y=0|\boldsymbol{x})}
\tag{5}

ロジスティック回帰のオッズ$o$は式(6)で表されます。

\begin{align}
o
&=\frac{\hat{y}}{1-\hat{y}}\\
&=\frac{\sigma(\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})}{1-\sigma(\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})}\\
&=\frac{(1+\exp{(-\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime}}))^{-1}}
{1-(1+\exp{(-\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime}}))^{-1}}\\
&=\frac{1}
{(1+\exp{(-\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})})-1}\\
&=\frac{1}
{\exp{(-\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})}}\\
&=\exp{(\boldsymbol{w}^{\prime\mathrm{T}}\boldsymbol{x}^{\prime})}
\end{align}
\tag{6}

なお、ロジスティック回帰のオッズの導出は、E資格試験では頻出です。

#ロジスティック回帰の学習
###尤度関数
学習データを$D=\{({\boldsymbol{x}}_{1},y_{1}),\cdots,({\boldsymbol{x}}_{N},y_{N})\}$と表します。
なお、$y_n\in \{0,1\}$とし、2つのクラスに割り当てます。
学習では、学習データ$D$を与え、最尤推定法によってモデルパラメータを求めます。

まず、尤度関数を定義します。
$\boldsymbol{x}$が与えられた条件で試行した結果、確率変数が$y\in \{0,1\}$、その条件付き確率が式(4)として与えられるとき、その確率変数はベルヌーイ分布に従います。
このため、尤度$L$は式(7)のベルヌーイ分布の確率質量関数で与えられます。

\begin{align}
L(\boldsymbol{w}^{\prime};D)
&=P(y=1|\boldsymbol{x})^y P(y=0|\boldsymbol{x})^{1-y}\\
&=\hat{y}^y (1-\hat{y})^{1-y}
\end{align}
\tag{7}

[ベルヌーイ分布][link-2]
[link-2]:https://qiita.com/fridericusgauss/items/14bbcdc0423fc39067c2

学習データ$D$の尤度関数$L(\boldsymbol{w}^{\prime};D)$は、各データの尤度の積で与えられるので、式(8)で表されます。

\begin{align}
L(\boldsymbol{w}^{\prime};D)
&=\prod_{n=1}^N \hat{y}_n^{y_n} (1-\hat{y}_n)^{1-y_n}
\end{align}
\tag{8}

負の対数尤度関数$L_{\log}(\boldsymbol{w}^{\prime}; D)$は式(9)で表されます。

\begin{array}{ll}
L_{\log}(\boldsymbol{w}^{\prime}; D)
 &= -\ln L(\boldsymbol{w}^{\prime}; D) \\
&= -\ln \prod_{n=1}^{N} \hat{y}_n^{y_n} (1-\hat{y}_n)^{1-y_n}\\
&= - \sum_{i=1}^{N} \ln (\hat{y}_n^{y_n} (1-\hat{y}_n)^{1-y_n})\\
&= - \sum_{i=1}^{N} (\ln \hat{y}_n^{y_n} +\ln (1-\hat{y}_n)^{1-y_n})\\
&= - \sum_{i=1}^{N} (y_n\ln \hat{y}_n +(1-y_n)\ln (1-\hat{y}_n))
\end{array}
\tag{9}

ただし、$\ln x$は自然対数$\log_{\mathrm{e}} x$です。

最尤推定法では、負の対数尤度関数$L_{\log}(\boldsymbol{w}^{\prime}; D)$が最小となるようなモデルパラメータを求めます。
なお、ロジスティック回帰における負の対数尤度関数$L_{\log}(\boldsymbol{w}^{\prime}; D)$は、2値クロスエントロピー誤差関数と等価になります。
これは、クラス分類における誤差関数としては、クロスエントロピー誤差関数が使用されることを表しています。

[クロスエントロピー][link-3]
[link-3]:https://qiita.com/fridericusgauss/items/0716ffd297c4f7b2dca7

###最尤推定法によるパラメータ推定
式(9)の負の対数尤度関数$L_{\log}(\boldsymbol{w}^{\prime}; D)$の最小化問題は、解析的に解けないため、最適化アルゴリズムを使用して数値的に解きます。
なお、この最適化問題においては、学習データがパラメータ(与えられて固定される)で、モデルパラメータが最適化変数(変化させる変数)ですので、ご注意ください。
機械学習では最適化アルゴリズムとして、勾配法がよく使用されます。
勾配法は、目的関数の勾配ベクトルを利用して、逐次的に解を更新する最適化アルゴリズムの枠組みです。

特に深層学習では、様々な勾配法が使用されますが、本記事ではロジスティック回帰の解説ですので、深層学習における最適化アルゴリズムについては別の記事でまとめる予定です。

ここでは勾配法の中でも、最も基本的なアルゴリズムである、最急降下法について解説します。
下記では、便宜上、目的関数である負の対数尤度関数を$L_{\log}(\boldsymbol{w})$と表記します。
最急降下法は、式(10)のように、最適化変数を逐次更新する方法です。

\begin{align}
\boldsymbol{w}:&=\boldsymbol{w}-\eta \boldsymbol{\nabla}(L_{\log})\\
\boldsymbol{\nabla}(L_{\log})
&=\frac{\partial L_{\log}}{\partial\boldsymbol{w}}
=\left(\frac{\partial L_{\log}}{\partial w_1},\cdots,\frac{\partial L_{\log}}{\partial w_M}\right)^{\mathrm{T}}
\end{align}
\tag{10}

ただし、$\eta>0$はステップ幅(学習率)、$\boldsymbol{\nabla}(L_{\log})$は勾配ベクトルです。
最急降下法は非常に有名なアルゴリズムなので、具体的な説明は割愛しますが、勾配ベクトルを逐次計算する必要がある点がポイントです。

###順伝播と逆伝播
ロジスティック回帰における、勾配ベクトルを導出してみます。
下記では、式(2)から、$\hat{y}$が$\boldsymbol{w}$の関数であることに着目し、$\hat{y}(\boldsymbol{w})$と明記します。
勾配ベクトルは式(11)のように変形できます。

\begin{align}
\boldsymbol{\nabla}(L_{\log})
&=\frac{\partial L_{\log}}{\partial\boldsymbol{w}}\\
&=- \frac{\partial}{\partial\boldsymbol{w}}
\sum_{i=1}^{N} (y_n\ln \hat{y}_n(\boldsymbol{w}) +(1-y_n)\ln (1-\hat{y}_n(\boldsymbol{w})))\\
&=- \sum_{i=1}^{N} 
\left(y_n \frac{\partial}{\partial\boldsymbol{w}} \ln \hat{y}_n(\boldsymbol{w}) 
+(1-y_n) \frac{\partial}{\partial\boldsymbol{w}}
\ln (1-\hat{y}_n(\boldsymbol{w}))\right)\\
&=- \sum_{i=1}^{N} 
\left(\frac{y_n}{\hat{y}_n(\boldsymbol{w})} \frac{\partial}{\partial\boldsymbol{w}} \hat{y}_n(\boldsymbol{w}) 
+ \frac{1-y_n}{1-\hat{y}_n(\boldsymbol{w})} \frac{\partial}{\partial\boldsymbol{w}}
 (1-\hat{y}_n(\boldsymbol{w}))\right)\\
&=- \sum_{i=1}^{N} 
\left(\frac{y_n}{\hat{y}_n(\boldsymbol{w})} \frac{\partial \hat{y}_n(\boldsymbol{w}) }{\partial\boldsymbol{w}} 
- \frac{1-y_n}{1-\hat{y}_n(\boldsymbol{w})} \frac{\partial \hat{y}_n(\boldsymbol{w})}{\partial\boldsymbol{w}}\right)\\
&=- \sum_{i=1}^{N} 
\left(\left(\frac{y_n}{\hat{y}_n(\boldsymbol{w})} 
- \frac{1-y_n}{1-\hat{y}_n(\boldsymbol{w})}\right) \frac{\partial \hat{y}_n(\boldsymbol{w})}{\partial\boldsymbol{w}}\right)
\end{align}
\tag{11}

次に、式(12)のように書き直すと、$\frac{\partial \hat{y}_n(\boldsymbol{w})}{\partial\boldsymbol{w}}$は式(13)で表せます。

\begin{align}
z_n=\boldsymbol{w}^{\mathrm{T}}\boldsymbol{x}_n\\
\hat{y}_n=\sigma(z_n)
\end{align}
\tag{12}

\begin{align}
\frac{\partial \hat{y}_n(\boldsymbol{w})}{\partial\boldsymbol{w}}
&=\frac{\partial \hat{y}_n}{\partial z_n}\frac{\partial z_n}{\partial\boldsymbol{w}}\\
&=\frac{\partial }{\partial z_n}(\sigma(z_n))\frac{\partial}{\partial\boldsymbol{w}}(\boldsymbol{w}^{\mathrm{T}}\boldsymbol{x}_n)\\
&=\sigma(z_n)(1-\sigma(z_n))\boldsymbol{x}_n\\
&=\hat{y}_n(\boldsymbol{w})(1-\hat{y}_n(\boldsymbol{w}))\boldsymbol{x}_n
\end{align}
\tag{13}

なお、計算過程で、式(14)のシグモイド関数の微分の性質を用いました。

\begin{align}
\frac{d }{d x}\sigma(x)=\sigma(x)(1-\sigma(x))
\end{align}
\tag{14}

【折り畳み】式(14)の導出

> >```math >\begin{align} >\frac{d }{d x}\sigma(x) =&\frac{d }{d x}\frac{1}{1+\exp{(-x)}}\\ >=&-\frac{1}{(1+\exp{(-x)})^2}\frac{d }{d x}(1+\exp{(-x)})\\ >=&\frac{\exp{(-x)}}{(1+\exp{(-x)})^2}\\ >=&\frac{1+\exp{(-x)}}{(1+\exp{(-x)})^2}-\frac{1}{(1+\exp{(-x)})^2}\\ >=&\frac{1}{1+\exp{(-x)}}-\left(\frac{1}{1+\exp{(-x)}}\right)^2\\ >=&\sigma(x)-\sigma^2(x)\\ >=&\sigma(x)(1-\sigma(x)) >\end{align} >``` >

よって、勾配ベクトルは式(15)のように変形できます。
なお、$\hat{y}(\boldsymbol{w})$を$\hat{y}$に書き直します。

\begin{align}
\boldsymbol{\nabla}(L_{\log})
&=- \sum_{i=1}^{N} 
\left(\left(\frac{y_n}{\hat{y}_n} 
- \frac{1-y_n}{1-\hat{y}_n}\right) \hat{y}_n(1-\hat{y}_n)\boldsymbol{x}_n\right)\\
&=- \sum_{i=1}^{N} 
\left((y_n(1-\hat{y}_n)- (1-y_n) \hat{y}_n)\boldsymbol{x}_n\right)\\
&=- \sum_{i=1}^{N} 
\left((y_n- \hat{y}_n)\boldsymbol{x}_n\right)\\
&=-\boldsymbol{X}^{\mathrm{T}}(\boldsymbol{y}-\hat{\boldsymbol{y}})\\
&=-\boldsymbol{X}^{\mathrm{T}}\boldsymbol{\delta}\\
\boldsymbol{X}&=(\boldsymbol{x}_1^{\prime} ,\cdots,\boldsymbol{x}_N^{\prime})^{\mathrm{T}}
\end{align}
\tag{15}

勾配ベクトルは目的変数の誤差$\boldsymbol{\delta}=\boldsymbol{y}-\hat{\boldsymbol{y}}$と説明変数で表せることがわかります。
つまり、最適化問題を解くためには、目的関数(誤差関数)の勾配ベクトルの計算が必要ですが、そのためにはモデルの最終的な出力(誤差)を計算する必要があります。

説明変数のデータ$\boldsymbol{X}$を入力したとき、ロジスティック回帰モデルの出力$\hat{\boldsymbol{y}}$は、式(16)のように、線形結合の出力$\boldsymbol{z}$を経由し、それを入力とした
シグモイド関数の出力$\boldsymbol{\sigma}(\boldsymbol{z})$です。

\begin{align}
\boldsymbol{z}=\boldsymbol{X}\boldsymbol{w}\\
\hat{\boldsymbol{y}}=\boldsymbol{\sigma}(\boldsymbol{z})
\end{align}
\tag{16}

このように、データを入力した後、中間変数を経由して、モデルの最終的な出力(あるいは誤差)を計算することを、__順伝播(forward)__といいます。

一方、勾配ベクトルは、誤差を計算した後、中間変数を経由して、連鎖微分によって間接的に得られます。
このように、順伝播とは逆方向に、誤差を入力した後、中間変数を経由して、誤差関数の勾配ベクトルを計算することを、__逆伝播(backward)__といいます。

ロジスティック回帰のような、活性化関数などを経由するモデルの場合、順伝播と逆伝播を逐次行うことで、学習を進めます。
深層学習を含むニューラルネットワークの学習も同様なので、順伝播と逆伝播の計算はE資格試験において頻出です。
さらに、E資格試験では、順伝播と逆伝播の実装問題も頻出ですので、結果を覚えるのではなく、導出も覚えておきましょう。

#おわりに
E資格向けのロジスティック回帰について解説しました。
なお、上記は、2021年2月時点における内容であることにご注意ください。

[E資格試験に関する私の投稿記事リスト][link-1]