NumPyでロジスティック回帰を実装する

ディープラーニング以前 (〜2010年) の機械学習について、はじパタを使って整理しています。

今回は「第6章線形識別関数」を参考にしながら、ロジスティック回帰をNumPyで実装してみます。

はじめてのパターン認識

作者: 平井有三
出版社/メーカー: 森北出版
発売日: 2012/07/31
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 7回
この商品を含むブログ (5件) を見る

ロジスティック回帰

実装に入る前に、ロジスティック回帰の理論を整理します。ここでは、2値 $\{C_1, C_2\}$ に分類するタスクを考えます。

シグモイド関数

まずクラス $C_1$ に属する確率を計算するための定式化を行います。

2値分類ですので、 $C_1$ の確率と $C_2$ の確率の合計は1です

入力 $\vec{x}$ (列ベクトル) が与えられたときにクラス $C_1$ に分類される事後確率は以下となります。

$P(C_1 | \vec{x}) = \frac{P(\vec{x} | C_1)P(C_1)}{P(\vec{x} | C_1)P(C_1) + P(\vec{x} | C_2)P(C_2)}$

このときaを導入すると、

$a := ln \frac{P(\vec{x} | C_1)P(C_1)}{P(\vec{x} | C_2)P(C_2)}$

事後確率 $P(C_1 | \vec{x})$ は以下のシンプルな式になります。

$P(C_1 | \vec{x}) = \frac{1}{1 + exp(-a)} = \sigma(a)$

$\sigma(a)$ をシグモイド関数といいます。シグモイド関数の値は [0, 1] となりますので、確率を表現するために好都合な関数です。

このシグモイド関数の逆関数をロジット関数といいます。

$a = ln \frac{\sigma(a)}{1 - \sigma(a)} = ln \frac{P(C_1 | \vec{x})}{P(C_2 | \vec{x})}$

ロジスティック回帰モデル

次にaの計算方法を定義します。

ここで重み $\vec{w}$ (列ベクトル) を導入し $a=\vec{w}^T\vec{x}$ で表現します。

$\sigma(a) = \frac{1}{1 + exp(-\vec{w}^T\vec{x})} = \frac{exp(\vec{w}^T\vec{x})}{1 + exp(\vec{w}^T\vec{x})}$

ロジスティック回帰モデルは、二項分布をリンク関数とした一般化線形モデルに分類されます。こちらの投稿も参考にしてみてください。緑本とstatsmodelsを使って、他の一般化線形モデルと比較しています。

Pythonで実装しながら緑本を学ぶ (第6章 GLMの応用範囲を広げる -ロジスティック回帰など-) - け日記

交差エントロピー誤差関数

最後に、重み $\vec{w}$ を教師データから獲得するために、最適化すべき尤度関数を定義します。

あるクラスに分類される確率で表現できる2値分類のため、尤度関数として交差エントロピー誤差関数を導入します。

$N$ は教師データ数
$t_i$ は教師データ
$y_i$ は入力 $\vec{x_i}$ に対して予測した確率

$L(\{y_1, ..., y_N\}) = \prod_{i=1}^{N} y_{i}^{t_i}(1 - y_{i})^{(1-t_{i})}$

この尤度関数を対数化して、負の符号をつけます。

$lL(\{y_1, ..., y_N\})=-\sum_{i=1}^{N}(t_{i} ln y_{i} + (1-t_i) ln (1-y_i))$

$y_i$ に代入します。

$lL(\vec{w})=-\sum_{i=1}^{N}(t_{i} \vec{w}^T \vec{x_i} - ln(1 + exp(\vec{w}^T\vec{x_i})))$

これを $\vec{w}$ で微分します。この $\frac{\partial lL}{\partial \vec{w}}$ が0となるように $\vec{w}$ を学習していきます。

$\frac{\partial lL}{\partial \vec{w}} = -\sum_{i=1}^N (t_i \vec{x_i} - \frac{\vec{x_i}exp(\vec{w}^T \vec{x_i})}{1 + exp(\vec{w}^T\vec{x_i})}) = \sum_{i=1}^N \vec{x_i}(y_i - t_i)$

最急降下法

学習方法として、今回は最急降下法で実装します。 $\alpha$ は更新量をコントロールするためのハイパパラメータです。

初期化: ランダムな重み $\vec{w_0}$ を設定
$w$ で勾配 $\frac{\partial dL}{\partial \vec{w}}$ を計算
$w$ <- $\vec{w} - \alpha \frac{\partial dL}{\partial \vec{w}}$ で重みを更新
収束するまで2.3.を繰り返す

NumPyで実装

scikit-learn付属のirisデータセットを使って、上を実装してみます。

入力は3次元 (萼片の長さ・幅とバイアス)
- 重み $\vec{w}$ も3次元
setosaとversicolorのデータのみを使う
- versicolorの場合に $t_i = 1$ とする

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# irisデータセットのロード
iris = load_iris()

# setosaとversicolorの萼片の長さ・幅のみを使う
X = iris.data[iris.target != 2][:, :2].T
X = np.concatenate([X, np.ones((1, X.shape[1]))])  # バイアス項
t = iris.target[iris.target != 2]

# 分布をチェック
plt.scatter(X[0][t == 0], X[1][t == 0], label=iris.target_names[0]) 
plt.scatter(X[0][t == 1], X[1][t == 1], label=iris.target_names[1])
plt.legend()
plt.show()

f:id:ohke:20190626111425p:plain

学習用のコードは以下です。上の式をナイーブに実装しています。

# wを乱数で初期化
def initialize_w():
    return np.random.normal(0, size=3).reshape((3, 1))

# versicolorの確率を予測
def predict_probability(w, X):
    return (np.exp(np.dot(w.T, X)) / (1.0 + np.exp(np.dot(w.T, X)))).reshape(-1)

# 損失関数
def loss_function(w, X, t):
    loss = 0

    for i in range(t.shape[0]):
        loss += t[i]*np.dot(w.T, X[:, i]) - np.log(1 + np.exp(np.dot(w.T, X[:, i])))

    return -loss

# 損失関数の勾配
def gradient_loss(w, X, t):
    y = predict_probability(w, X)
    n = t.shape[0]
    loss = np.zeros_like(w)
  
    for i in range(n):
        loss += ((y[i] - t[i]) * X[:, i])[:, np.newaxis]

    return loss

# 最急降下法で学習
def gradient_descent(X, t, alpha=0.001):
    w_t = initialize_w()
    loss_t = loss_function(w_t, X, t)

    w_history, loss_history = [w_t], [loss_t]

    for i in range(1000000):  # lossに変化がなくなるまでループする
        w_t -= gradient_loss(w_t, X, t) * alpha
        loss_t = loss_function(w_t, X, t)

        # 収束判定
        if loss_t >= loss_history[-1]:
            break

        w_history.append(w_t)
        loss_history.append(loss_t)
        
    return w_history, loss_history


# 学習
w_history, loss_history = gradient_descent(X, t)

損失関数の値の推移はこちらで、最初にぐっと下がって、あとはダラダラと少しずつ下げていってることがわかります。

f:id:ohke:20190626115100p:plain

得られた重みを使って分類させてみましたところ、精度100%で完全に分離できていることを確認できました。

# 重みの確認
w_hat = w_history[-1]
print(w_hat)
# [[ 12.30194381]
#  [-12.92596168]
#  [-26.27142723]]

# 得られた重みで確率を計算
y_probability = predict_probability(w_hat, X)

# 精度
y = np.where(y_probability > 0.5, 1, 0)
sum(y == t) / t.shape[0]  # 1.0

まとめ

今回はNumPyでナイーブなロジスティック回帰を実装しました。