前書き
記事の概要
今回は教師あり学習(分類問題)の有名な手法である ロジスティック回帰 について解説していきます。
この記事には
- 数式によるロジスティック回帰の解説
- プログラム実装例
が含まれています。
対象読者
大学1年生で勉強する「線形代数学」および「統計学」の基礎事項は既知であるとして話を進めていきます。
キーワードは 最尤推定 / ベルヌーイ分布 / 交差エントロピー です。
もし記事中でわからない言葉が出てきた場合には、線形代数学・統計学の教科書を読んだり他のサイトを見たりして調べてみてください。
導入
教師あり学習における2値分類問題について考えます。
2値分類問題では、各データ ${\bf x}$ に対してラベル $y \in \{ 0, 1 \}$ が付与されているとします。$y = 0$ のときには ${\bf x}$ はクラス $C_0$ に属しているとし、$y = 1$ のときには ${\bf x}$ はクラス $C_1$ に属していると解釈します。このとき、与えられた ${\bf x}$ に対してそのラベルを予測するのが2値分類問題です。
ロジスティック回帰 モデルではデータ ${\bf x}$ のラベル $y = 1$ である確率(クラス $C_1$ に属する確率)を予測として出力します。確率ではなく予測ラベルを出力する必要がある場合には、予測確率が閾値 $\theta$ 以上ならば $\hat{y} = 1$ 、そうでなければ $\hat{y} = 0$ とします。この予測確率はパラメータ ${\bf w}$ を用いて
$$
\hat{y} = \sigma({\bf w}^T {\bf x})
$$
とモデリングします。ここで $\sigma$ は (ロジスティック)シグモイド関数 と呼ばれる関数です。
モデルに含まれるパラメータ ${\bf w}$ は与えられた学習データ集合 $D = \{ ({\bf x_n}, y_n) | y_n \in \{ 0, 1 \} \}_{n=1}^N$ から推定します。ただし、一般線形回帰の場合とは異なり最尤解が解析的に求められないため、反復法を利用することになります。
理論
基礎的な話
識別的なアプローチ
観測データ集合 $D = \{ ({\bf x_n}, y_n) | y_n \in \{ 0, 1 \} \}_{n=1}^N$ が与えられたとき、$y = 1$ となる確率 $\hat{y}$ を
$$
\hat{y} = \sigma({\bf w}^T {\bf x})
$$
とします。 $\sigma$ は(ロジスティック)シグモイド関数で
$$
\sigma(x) \equiv \frac{1}{1 + \exp(- x)}
$$
という式で表されます。グラフを描いてみると
という感じのS字カーブになります。実数全体を定義域、 $[0, 1]$ を値域とする単調連続関数で、点 $(0.0, 0.5)$ に関して点対称な形をしています。 $\lim_{x \to \infty} \sigma(x) = 1, \lim_{x \to - \infty} \sigma(x) = 0$ の性質から、 $y = 1$ のデータに対しては ${\bf w}^T {\bf x}$ を大きく、 $y = 0$ のデータに対しては ${\bf w}^T {\bf x}$ を小さくすることが目標となります。
このモデリングにおいて観測データをうまく説明できるようにパラメータ ${\bf w}$ を調整することになりますが、その尤もらしさは 交差エントロピー誤差 で測ることにします。
交差エントロピー誤差
一般に( $K$ 値離散確率分布に関する) 交差エントロピー(クロスエントロピー) ${\mathcal H}({\bf p}, {\bf q})$ は
$$
{\mathcal H}({\bf p}, {\bf q}) = - \sum_{k=1}^K p_k \log{q_k}
$$
と定義されます。ここで ${\bf p}$ は真の確率分布、 ${\bf q}$ は推定された確率分布を表しています。
交差エントロピーはKLダイバージェンス $D_{KL}(\cdot \| \cdot)$ を用いて
$$
{\mathcal H}({\bf p}, {\bf q}) = {\mathcal H}({\bf p}) + D_{KL}({\bf p} \| {\bf q})
$$
と書くことができるので、 ${\bf p}$ が固定されている下では 推定された ${\bf q}$ が ${\bf p}$ とどれくらい離れているかを表す指標になっています。
機械学習の文脈では、観測されたクラスラベル ${\bf y}$ を真の確率 ${\bf p}$ 、予測したクラスラベルの確率分布 ${\bf \hat{y}}$ を ${\bf q}$ として評価します。このとき、 交差エントロピー誤差 は
$$
{\mathcal H}({\bf y}, {\bf \hat{y}}) = - \sum_{k=1}^K y_k \log{\hat{y_k}}
$$
となります。 ${\bf y}$ はone hot encodingで表現されており、1つの成分 $l$ に限って $y_l = 1$ でそれ以外の成分は $y_k = 0 \ (k \ne l)$ です。そのため、実質的には
$$
{\mathcal H}({\bf y}, {\bf \hat{y}}) = - \log{\hat{y_l}}
$$
という計算をしていることになります。この式から、正解のクラスラベルに大きな確率を割り当てると誤差は小さくなる(確率1を割り振ると誤差は0になる)ことがわかります。特に $K = 2$ の場合には観測されたクラスラベルはbinaryで表現すれば良いので
$$
{\mathcal H}(y, \hat{y}) = - y \log{\hat{y}} - (1 - y) \log{(1 - \hat{y})}
$$
となります。
誤差関数の最小化
では、(2値)交差エントロピー誤差関数を最小化するようなパラメータ ${\bf w}$ を求めていきます。
誤差関数 $L({\bf w})$ は
$$
L({\bf w}) = - \sum_{n=1}^N { y_n \log{\hat{y_n}} + (1 - y_n) \log{(1 - \hat{y_n})} }
$$
となります。 $\hat{y_n}$ は ${\bf x_n}$ に対する予測値で $\hat{y_n} = \sigma({\bf w}^T {\bf x_n})$ です。これを ${\bf w}$ で微分すると
\begin{eqnarray*}
\frac{\partial}{\partial {\bf w}} L
&=& - \frac{\partial}{\partial {\bf w}} \sum_{n=1}^N \{ y_n \log{ \hat{y_n} } + (1 - y_n) \log{ (1 - \hat{y_n}) }\} \\
&=& - \sum_{n=1}^N \left\{ y_n \frac{\partial}{\partial {\bf w}} \log{ \hat{y_n} } + (1 - y_n) \frac{\partial}{\partial {\bf w}} \log{ (1 - \hat{y_n}) } \right\} \\
&=& - \sum_{n=1}^N \left\{ \frac{y_n}{ \hat{y_n} } \frac{\partial}{\partial {\bf w}} \hat{y_n} - \frac{1 - y_n}{1 - \hat{y_n}} \frac{\partial}{\partial {\bf w}} \hat{y_n} \right\} \\
&=& - \sum_{n=1}^N \left\{ \left(\frac{y_n}{ \hat{y_n} } - \frac{1 - y_n}{1 - \hat{y_n}} \right) \frac{\partial}{\partial {\bf w}} \hat{y_n} \right\} \\
\end{eqnarray*}
が得られます。表記の都合上 $v_n = {\bf w}^T {\bf x_n}$ とすると
\begin{eqnarray*}
\frac{\partial}{\partial {\bf w}} \hat{y_n}
&=& \frac{\partial}{\partial {\bf w}} \sigma(v_n) \\
&=& \frac{\partial v_n}{\partial {\bf w}} \frac{\partial}{\partial v_n} \sigma(v_n) \\
&=& {\bf x_n} \frac{\partial}{\partial v_n} \sigma(v_n) \\
&=& {\bf x_n} \sigma(v_n) (1 - \sigma(v_n)) \\
&=& \hat{y_n} (1 - \hat{y_n}) {\bf x_n}
\end{eqnarray*}
が成立します。式変形の途中でシグモイド関数の微分に関する公式
$$
\frac{\partial}{\partial x} \sigma(x) = \sigma(x) \left\{ 1 - \sigma(x) \right\}
$$
を利用しました。導出過程は折りたたんであります。
シグモイド関数の微分公式の導出
\begin{eqnarray*}
\frac{\partial}{\partial x} \sigma(x)
&=& \frac{\partial}{\partial x} \frac{1}{1 + \exp{(-x)}} \\
&=& \frac{\exp{(-x)}}{(1 + \exp{(-x)})^2} \\
&=& \frac{1}{1 + \exp{(-x)}} \cdot \frac{\exp{(-x)}}{1 + \exp{(-x)}} \\
&=& \frac{1}{1 + \exp{(-x)}} \cdot \left( 1 - \frac{1}{1 + \exp{(-x)}} \right) \\
&=& \sigma(x) (1 - \sigma(x))
\end{eqnarray*}
以上のようにしてシグモイド関数の微分の公式を導出することができます。シグモイド関数の微分の計算には入力 $x$ は必要なく、出力 $y = \sigma(x)$ さえあれば良いです。
したがって
\begin{eqnarray*}
\frac{\partial}{\partial {\bf w}} {\mathcal L}
&=& \sum_{n=1}^N \left(\frac{y_n}{ \hat{y_n} } - \frac{1 - y_n}{1 - \hat{y_n}} \right) \hat{y_n} (1 - \hat{y_n}) {\bf x_n} \\
&=& - \sum_{n=1}^N \{ y_n (1 - \hat{y_n}) - (1 - y_n) \hat{y_n} \} {\bf x_n} \\
&=& - \sum_{n=1}^N (y_n - \hat{y_n}) {\bf x_n} \\
&=& \sum_{n=1}^N (\hat{y_n} - y_n) {\bf x_n} \\
\end{eqnarray*}
となっていることがわかります。
この勾配を ${\bf 0}$ にするパラメータ ${\bf w^{\ast}}$ を求めたいですが、方程式 $\frac{\partial}{\partial {\bf w}} {\mathcal L} = {\bf 0}$ を解析的に解くことはできません。シグモイド関数が複雑すぎるためです。
そのため、最適解 ${\bf w^{\ast}}$ を反復法(勾配法)によって求めます。具体的な更新式は
$$
{\bf w} := {\bf w} - \eta \frac{\partial}{\partial {\bf w}} {\mathcal L}
= {\bf w} - \eta \sum_{n=1}^N (\hat{y_n} - y_n) {\bf x_n} \
$$
です。( $\eta$ は学習率を表しています。)
発展的な話
生成的なアプローチ
この項では、ロジスティック回帰を確率推論の枠組みから考えることで交差エントロピー誤差が自然に現れることを見ていきたいと思います。
観測データ集合 $D = \{ ({\bf x_n}, y_n) | y_n \in \{ 0, 1 \} \}_{n=1}^N$ が与えられたとき、それを表現するモデルとして
$$
y \sim \text{Ber}(p) \text{ where } p \equiv \hat{y} = \sigma({\bf w}^T {\bf x})
$$
を採用します。以下ではモデルに含まれるパラメータ ${\bf w}$ を最尤推定することを考えます。
各データ組 $({\bf x_n}, y_n)$ が独立に得られたと仮定すると、観測データ集合 $D$ に対する負の対数尤度 ${\mathcal L({\bf w})}$ は
\begin{eqnarray*}
{\mathcal L({\bf w})}
&=& - \log{ \prod_{n=1}^N p(y_n | {\bf x_n}, {\bf w}) } \\
&=& - \log{ \prod_{n=1}^N p^{y_n} (1 - p)^{1 - y_n} } \\
&=& - \sum_{n=1}^N \{ y_n \log{ p } + (1 - y_n) \log{ (1 - p) } \} \\
&=& - \sum_{n=1}^N \{ y_n \log{ \hat{y_n} } + (1 - y_n) \log{ (1 - \hat{y_n}) } \}
\end{eqnarray*}
と求めることができます。よく見ると、この式は 交差エントロピー誤差 に等しくなっています。
「一般線形回帰モデルにおける負の対数尤度の最小化」が「二乗和誤差の最小化」と等価であったように、「ロジスティック回帰モデルにおける負の対数尤度の最小化」が「交差エントロピー誤差の最小化」と等価になっています。(一般線形回帰モデルについては 過去の記事 を参照してください。)
シグモイド関数の導入について
ロジスティック回帰モデルでは確率値の出力のためにシグモイド関数を導入しました。しかし、これ以外にも値を区間 $[0, 1]$ に収めるような関数の候補はたくさんあるはずです。この項では、シグモイド関数を導入する背景について見ていきます。
生成的なアプローチによるロジスティック回帰の定式化をした際、ラベル $y$ はベルヌーイ分布に従うと仮定していましたが、データ ${\bf x}$ に関しては何も前提を置いていませんでした。そこで、データとラベルの組 $({\bf x}, y)$ に関して
\begin{eqnarray*}
y &=& 0 \Rightarrow {\bf x} \sim {\mathcal N}({\bf \mu_0}, \Sigma) \\
y &=& 1 \Rightarrow {\bf x} \sim {\mathcal N}({\bf \mu_1}, \Sigma)
\end{eqnarray*}
という仮定を置くことにします。つまり、クラス $C_0$ に属するデータ ${\bf x}$ はガウス分布 ${\mathcal N}({\bf \mu_0}, \Sigma)$ に従うとし、クラス $C_1$ に属するデータ${\bf x}$ はガウス分布 ${\mathcal N}({\bf \mu_1}, \Sigma)$ に従うと考えます。 共分散行列 $\Sigma$ が共通であるというのがとても重要です。
求めるべきは(データ ${\bf x}$ が与えられたときに) $y = 1$ となる確率でした。この条件付き確率を $\Pr(y=1 | {\bf x})$ と書くことにします。するとベイズの定理より
\begin{eqnarray*}
\Pr(y=1 | {\bf x})
&=& \frac{ \Pr(y=1) \Pr({\bf x} | y=1) }{ \Pr(y=1) \Pr({\bf x} | y=1) + \Pr(y=0) \Pr({\bf x} | y=0)} \\
&=& \frac{1}{1 + \frac{\Pr(y=0) \Pr({\bf x} | y=0)}{\Pr(y=1) \Pr({\bf x} | y=1)}} \\
&=& \frac{1}{1 + \exp{ \left( - \log{ \left( \frac{\Pr(y=1) \Pr({\bf x} | y=1)}{\Pr(y=0) \Pr({\bf x} | y=0)} \right) } \right) } }
\end{eqnarray*}
が成立します。分母第2項についてもう少し見てみましょう。仮に $v = \log{ \left( \frac{\Pr(y=1) \Pr({\bf x} | y=1)}{\Pr(y=0) \Pr({\bf x} | y=0)} \right) }$ とおけば
\begin{eqnarray*}
v
&=& \log{ \left( \frac{\Pr(y=1) \Pr({\bf x} | y=1)}{\Pr(y=0) \Pr({\bf x} | y=0)} \right) } \\
&=& \log{\Pr(y=1)} + \log{\Pr({\bf x} | y=1)} - \log{\Pr(y=0)} - \log{\Pr({\bf x} | y=0)}
\end{eqnarray*}
となります。 ${\bf x}$ に関係する項だけに注目すると
\begin{eqnarray*}
v
&=& \log{\Pr({\bf x} | y=1)} - \log{\Pr({\bf x} | y=0)} + const \\
&=& - \frac{1}{2} ({\bf x} - {\bf \mu_1})^T \Sigma^{-1} ({\bf x} - {\bf \mu_1}) + \frac{1}{2} ({\bf x} - {\bf \mu_0})^T \Sigma^{-1} ({\bf x} - {\bf \mu_0}) + const \\
&=& - \frac{1}{2} \left( {\bf x}^T \Sigma^{-1} {\bf x} - 2 {\bf x}^T \Sigma^{-1} {\bf \mu_1} + {\bf \mu_1}^T \Sigma^{-1} {\bf \mu_1} \right) + \frac{1}{2} \left( {\bf x}^T \Sigma^{-1} {\bf x} - 2 {\bf x}^T \Sigma^{-1} {\bf \mu_0} + {\bf \mu_0}^T \Sigma^{-1} {\bf \mu_0} \right) + const \\
&=& {\bf x}^T \Sigma^{-1} ({\bf \mu_1} - {\bf \mu_0}) + const
\end{eqnarray*}
というように簡約化することができます。クラス $C_0$ に属するデータの共分散行列とクラス $C_1$ に属するデータの共分散行列が共通だったことで ${\bf x}$ に関する2次の項が消去されました。これで $v$ は ${\bf x}$ に関する線形結合で表せることがわかります。 ${\bf w} = \Sigma^{-1} ({\bf \mu_1} - {\bf \mu_0})$ とすれば
$$
v = {\bf w}^T {\bf x}
$$
となり、予測確率 $\Pr(y=1 | {\bf x})$ は
\begin{eqnarray*}
\Pr(y=1 | {\bf x})
&=& \frac{1}{1 + \exp{(- v)} } \\
&=& \sigma(v)
\end{eqnarray*}
であることが従います。
こうした背景を踏まえると、線形結合 ${\bf w}^T {\bf x}$ にシグモイド関数を適用したものを予測確率とする理由が少しは納得できるかと思います。
実装
ロジスティック回帰の実装はGitHubにあります。
fit メソッドでパラメータを最適化しますが、ここでは単純な最急降下法を利用しています。(パラメータ max_iter で更新を行う回数を設定することができます。)
predict メソッドで予測ラベルを、 predict_proba メソッドでラベルが1になる確率 $Pr(y = 1 | {\bf x})$ を出力するようになっています。(scikit-learnの predict_proba メソッドとは少しインターフェースが異なります。)