ロジスティック回帰による2値分類 - 薬剤師のプログラミング学習日記

今回はロジスティック回帰を使ってデータを分類してみようと思います。「回帰」とありますが、分類タスクで使用されるアルゴリズムです。なお、この記事はCourseraのMachine Learning（機械学習）講座『Logistic Regression』の内容を基に書いています。また、コードは講座の課題をPythonで書き直したものです。

2値の分類問題とは
ロジスティック回帰のモデル表現
ロジスティック回帰の目的関数
最急降下法を適用する
プログラムの概要とPythonによる実装
オーバーフィッティングを正則化で回避する
ロジスティック回帰と正則化まとめ
参考

2値の分類問題とは

Eメールがスパムかスパムでないか？
商用サイトのオンライン取引が不正か不正でないか？
腫瘍が悪性か良性か？

のように1または0、陽性または陰性など2つの値だけを取るような分類問題です。

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

仮説関数(hypothesis function) $h_\theta (x)$ がどのようになるかを考えます。
線形回帰の場合は $h_\theta (x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \cdot \cdot \cdot +\theta_n x_n = \theta^T x$ でした。

ロジスティック回帰の場合は、シグモイド関数(ロジスティック関数)と言われる以下のような形をした関数 $\displaystyle{g(z) = \frac{1}{1 + e ^{-z}}}$ を使って次のように定義します。
f:id:enokisaute:20200613113209p:plain
$h_\theta (x) = g( \theta^T x)$

このシグモイド関数は $\theta^T x$ の値を0～1の範囲に押し込めることができ、これによって $h_\theta (x)$ が出力した値を"1"(または陽性)になる確率として解釈することができます。
例えば、腫瘍が悪性か良性かを分類するタスクで、腫瘍の特徴X(サイズ)を与えたとき、 $h_\theta (x)$ の返す出力が0.7だったとしたら、y=1である確率は0.7というように解釈します。
$0 \le h_{\theta} (x) \le 1$ なので、0.5を閾値として

$h_\theta (x) \ge 0.5$ ならy=1(悪性)
$h_\theta (x) \lt 0.5$ ならy=0(良性)

と分類します。

ロジスティック回帰の目的関数

データと $h_\theta (x)$ の返す出力のずれが大きければ返す値も大きく、両者のずれが小さければ小さい値を返すことのできる以下のような関数を考えます。

f:id:enokisaute:20200613183507p:plain

$Cost(h_\theta (x), y) = \left\{ \begin{array}{ll} y=1のとき ~~ -\log(h_\theta (x)) \\ y=0のとき ~~ -\log(1 - h_\theta (x)) \end{array} \right.$

どちらも $h_\theta (x)$ 、すなわち予測と $y$ が近ければ $Cost$ (コスト)が0に近づき、 $h_\theta (x)$ の出力と $y$ の乖離が大きいほどコストは大きく( $∞$ )なります。
これらをy=0、y=1の両方の場合を含んだ1つの式にまとめると、最終的に目的関数*1は次のようになります。

$\displaystyle{J(\theta) = \frac{1}{m} \sum_{i=1}^m ~[ -y^{(i)} \log(h_\theta (x^{(i)})) - (1 - y^{(i)}) \log(1 - h_\theta (x^{(i)})) ] }$

この目的関数が返すものは、訓練データの一つ一つ1からmまでについてのデータ全体におけるコストです。

最急降下法を適用する

あとはこの目的関数を最小にするパラメータ $\theta$ を見つけるだけです。
最急降下法でのパラメータ $\theta$ の更新式は以下のように書けます。

$\theta_j := \displaystyle{ \theta_j - \alpha\frac{\partial}{\partial \theta_j}J(\theta)}$

目的関数の偏微分はこちら

$\displaystyle{ \frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m} \sum_{i=1}^m (h_\theta (x^{(i)}) - y^{(i)}) x_j ^{(i)}}$

最急降下法のアルゴリズムについては、線形回帰のときに書いたので、詳しくはこちらをご参照ください。

線形回帰でデータに直線を当てはめる

ロジスティック回帰の勾配降下法の更新式は線形回帰の更新式と同じに見えますが、 $h_\theta (x)$ の定義が異なるため、実際の計算式は異なります。

プログラムの概要とPythonによる実装

学生が入学試験をパスできるかどうかを予測(合格 or 不合格)するという内容です。過去の志願者の点数の履歴データを訓練データとして用いてロジスティック回帰モデルを構築し、新たな志願者について2つの試験の点数から判断します。
データ(ex2data1.txt)の内容は次のようになっており*2、プロットすると下図のようになります。

[[34.62365962 78.02469282]
 [30.28671077 43.89499752]
 [35.84740877 72.90219803]
 [60.18259939 86.3085521 ]
・・・

f:id:enokisaute:20200613224128p:plain
軸はそれぞれ2つの試験の点数で合格と不合格が異なるマーカーで表示されています。
以下のコードで、ロジスティック回帰モデルを構築しこれらのデータに適合するよう訓練します。

import numpy as np
import matplotlib.pyplot as plt
 
# ロジスティック回帰モデル
 
# xが行列の場合, すべての要素でシグモイド関数が実行される必要がある。
def sigmoid(x):
    return 1 / (1 + np.exp(-x))
 
def compute_cost(theta, x, y):
    y_ = y.reshape(m, 1)    # ベクトルとして 1 - y, y * log(h)がそのまま計算できる
    h = sigmoid(np.dot(x, theta))   # 予測
    j = 1 / m * np.sum(-y_ * np.log(h) - (1 - y_) * np.log(1 - h))
    grad = 1 / m * np.dot(x.T, h - y_)
    return j, grad
 
# 勾配はcompute_costで求めているのでそれを使う
def gradient_descent(x, y, theta, alpha, num_iters):
    j_history = np.zeros(num_iters)
    for i in range(0, num_iters):
        j_history[i], grad = compute_cost(theta, x, y)
        theta = theta - alpha * grad
    return j_history, theta
 
 
if __name__ == '__main__':
    data = np.loadtxt('./data/ex2data1.txt', delimiter=',', dtype=float)
    X = data[:, 0:2]
    y = data[:, 2]
 
    # 訓練データのプロット
    pos = X[data[:, 2] == 1]  # positive
    neg = X[data[:, 2] == 0]  # negative
    plt.scatter(pos[:, 0], pos[:, 1], color='red', marker='x', label='admitted')
    plt.scatter(neg[:, 0], neg[:, 1], color='blue', marker='o', label='Not admitted')
    plt.xlabel('Exam 1 score')
    plt.ylabel('Exam 2 score')
    plt.title('Scatter plot of training data')
    plt.legend()
    # plt.show()
 
    m, n = X.shape  # m:トレーニングデータ数, n:特徴量の数
    intercept_term = np.ones((m, 1)) 
    X = np.hstack((intercept_term, X.reshape(m, -1)))
    initial_theta = np.zeros((n + 1, 1))
 
    # initial_thetaのコスト
    cost, grad = compute_cost(initial_theta, X, y)
    print('initial_cost: ', cost)
 
    alpha = 0.001        # 学習率
    num_iters = 300000
    # 最急降下法でthetaを更新していく
    cost, theta = gradient_descent(X, y, initial_theta, alpha, num_iters)
    print('theta: \n', theta)
    print('cost: ', cost[-1])
 
    # 決定境界の描画
    exam1 = np.arange(0, 100, 1)
    plt.plot(exam1, -theta[1] / theta[2] * exam1 - theta[0] / theta[2])
    plt.xlim(30, 100)
    plt.ylim(30, 100)
 
    # イテレーション毎のコストをプロット
    plt.figure(2)
    plt.plot(cost)
    plt.xlim(0,)
    plt.xlabel('iteration')
    plt.ylabel('J(theta)')
    plt.show()
 
    # トレーニングデータにおける正解率
    predict = np.where(sigmoid(np.dot(X, theta)) >= 0.5, 1, 0)
    print('Train Accuracy: ', np.mean(y.reshape(m, 1) == predict))

課題ではパラメータベクトル $\theta$ を求める部分はOctaveの組み込み最適化ソルバーを使いましたが、ここでは実際に最急降下法で $\theta$ を求めています。
gradient_descent()の学習率αと繰り返し回数num_itersをいろいろ試しながらやってみましたが、なかなか収束しなかったり、発散したりしてうまく値が求まらず、最初は自分の書いたコードが間違っているのかと思いました。そこで、参考文献の方のコードを見ると、繰り返し回数を30万回に設定しており、やっと解決。収束するまでにこれほどの更新が必要とは思いませんでした。匙加減が難しいとのことです。

以下は、先の訓練データにおいて決定境界( $\theta^T x = 0$ )を描画したものと目的関数 $J(\theta)$ を $\theta$ を更新する毎にプロットしたものです。
f:id:enokisaute:20200614004245p:plain

オーバーフィッティングを正則化で回避する

機械学習モデルでは、訓練に使ったデータではなく未知のデータに対して性能(汎化性能)を発揮する必要があります。しかし、訓練データの数に対してモデルの表現力が高すぎると、過剰に訓練データに適合(オーバーフィッティング、または過学習)してしまい、この汎化性能がうまく得られません。
ここでは、この過剰な適合を回避するために正則化(Regularization)という方法を用います。
正則化を導入したロジスティック回帰の目的関数は次のようになります。

$\small{\displaystyle{J(\theta) = \frac{1}{m} \sum_{i=1}^m [ -y^{(i)} \log(h_{\theta} (x^{(i)})) - (1-y^{(i)}) \log(1-h_{\theta}(x^{(i)}))] + \frac{\lambda}{2m} \sum_{j=1}^n \theta_j^2}}$

( $\theta_0$ 以外の)パラメータ $\theta$ の値が極端な重みを取ることに対してペナルティを課すことで過学習を抑制する、という方法です。ここで、λ(ラムダ)は正則化の強度を決定するパラメータです。
λが小さ過ぎると過学習はあまり抑制されず良い汎化性能が得られません。逆にλが大きければ強く抑制されますが、大き過ぎると今度は訓練データにも適合できなくなり、結局汎化性能は悪くなります。要するに、モデルの複雑さと単純さはトレードオフの関係にあり、ちょうど良いほどほどの地点を探るということになります。
続いて、偏微分です。

$\small{\displaystyle{ \frac{\partial}{\partial \theta_0} J(\theta) = \frac{1}{m} \sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)}) x_j^{(i)}}} ~~~~~ (j = 0)$

$\small{\displaystyle{ \frac{\partial }{\partial \theta_j}J(\theta) = \biggl( \frac{1}{m} \sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)}) x_j^{(i)} \biggr) + \frac{\lambda}{m} \theta_j}} ~~~~~ (j \ge 1)$

下のコードは、2つのテスト結果を持つ製造工場のマイクロチップが品質テストに合格するかどうかを予測するというものです。やりたいことはだいたい先のプログラムと同じ(合格 or 不合格を分類)ですが、今度のモデルは2次元の特徴量を6次の項(28次元)まで拡張することにより複雑な決定境界を持ち、非線形分離することが可能となります。

import numpy as np
import matplotlib.pyplot as plt
 
# 正則化を導入したロジスティック回帰
 
def map_feature(x1, x2, degree=6):
    m = len(x1)
    out = np.ones((m, 1))
    for i in range(1, degree + 1):
        for j in range(0, i + 1):
            temp = (x1 ** (i - j) * x2 ** j).reshape((m, 1))
            out = np.hstack((out, temp))
    return out
 
def sigmoid(x):
    return 1 / (1 + np.exp(-x))
 
# 正則化を導入
def compute_cost(theta, x, y, lam):
    m = len(x)
    y_ = y.reshape(m, 1)
    h = sigmoid(np.dot(x, theta))
    j = 1 / m * np.sum(-y_ * np.log(h) - (1 - y_) * np.log(1 - h)) \
        + lam / (2 * m) * np.sum(theta[1:] ** 2)
    reg = lam / m * theta[1:]
    grad = 1 / m * np.dot(x.T, h - y_)
    grad[1:] += reg
    return j, grad
 
def gradient_descent(x, y, theta, alpha, num_iters, lam):
    j_history = np.zeros(num_iters)
    for i in range(0, num_iters):
        j_history[i], grad = compute_cost(theta, x, y, lam)
        theta = theta - alpha * grad
    return j_history, theta
 
 
if __name__ == '__main__':
    data = np.loadtxt('./data/ex2data2.txt', delimiter=',', dtype=float)
    X = data[:, 0:2]
    y = data[:, 2]
 
    # 訓練データのプロット
    pos = X[data[:, 2] == 1]  # positive
    neg = X[data[:, 2] == 0]  # negative
    plt.scatter(pos[:, 0], pos[:, 1], color='red', marker='x', label='y=1')
    plt.scatter(neg[:, 0], neg[:, 1], color='blue', marker='o', label='y=0')
    plt.xlabel('Microchip Test 1')
    plt.ylabel('Microchip Test 2')
    plt.title('Plot of training data')
    plt.legend()
    # plt.show()
 
    # 各データポイントからより多くの特徴量を作成する.
    # 2次元→28次元のベクトルに変換する. これによって分類器はより複雑な決定境界を持つことができるが,
    X = map_feature(X[:, 0], X[:, 1])
 
    m, n = X.shape      # m:トレーニングデータ数, n:特徴量の数
    initial_theta = np.zeros((n, 1))
    alpha = 0.1    # 学習率
    num_iters = 3000000
    lam = 1    # ラムダ 0, 1, 100で試す
 
    cost, theta = gradient_descent(X, y, initial_theta, alpha, num_iters, lam)
    print("theta: ", theta)
    print("cost: ", cost[-1])
 
    # 決定境界の描画
    test1 = np.linspace(-1, 1.5, 100)
    test2 = np.linspace(-1, 1.5, 100)
    t1, t2 = np.meshgrid(test1, test2)
    z = np.zeros((len(t1), len(t2)))
 
    for i in range(len(test1)):
        for j in range(len(test2)):
            x1 = np.array([t1[i, j]])
            x2 = np.array([t2[i, j]])
            z[i, j] = np.dot(map_feature(x1, x2), theta)
 
    plt.contour(test1, test2, z, levels=[0])
    # グラフのタイトル. λに応じて入れ替えた.
    plt.title('Training data with decision boundary ($λ=1$)')
    # plt.title('No regularization (Overfitting) ($λ=0$)')
    # plt.title('Too much regularization (Underfitting) ($λ=100$)')
    plt.ylim([-0.8, 1.2])
    plt.legend()
 
    # イテレーション毎のコストをプロット
    plt.figure(2)
    plt.plot(cost)
    plt.xlabel("iteration")
    plt.ylabel("J(theta)")
    plt.show()
 
    # トレーニングデータにおける正解率
    predict = np.where(sigmoid(np.dot(X, theta)) >= 0.5, 1, 0)
    print("Train Accuracy: ", np.mean(y.reshape(m, 1) == predict))

データをmap_feature関数に通すことで、入力された2つの特徴量を $x_1, x_2, x_1^2, x_2^2, x_1x_2, ...$ などからなる、より多くの特徴量を持つ配列に変換しています。
多くの特徴量を持つことでより表現力のあるモデルになりますが、これによってオーバーフィッティングもしやすくなります。
決定境界の変化を確認するために正則化パラメータλを1、0(正則化なし)、100の3種類で試しています。
f:id:enokisaute:20200614193249p:plain
λ=1です。訓練データの正解率は83%。決定境界の複雑さが適度に抑制されています。