機械学習アルゴリズム実装シリーズ vol 2[ロジスティック回帰編]

こんにちは。ブレインズテクノロジー樋口です。
機械学習アルゴリズムを一から実装するシリーズ2回目。今回のお題はロジスティック回帰(Logistic Regression)の実装。
前回と同様あまり難しい説明はなしに実装していくつもりですが、もっと詳しいこと教えんかい( ﾟдﾟ )ｸﾜｯ!!となったら以下のような専門書をどうぞ。
machinelearningmastery.com

ロジスティック回帰とは

ロジスティック回帰(Logistic Regression)は分類問題を解決するための手法で教師あり学習に分類されます。
元々は統計の分野で発展してきたもので、応用先としては
・e-mailのスパム判定
・病理標本の悪性判定
・トランザクションの正常性判定

などに使われているようです。
# 名前には回帰(Regression)とついてるんですけど分類なんですね。。。

ロジスティック回帰[理論編]

前回、線形回帰ではモデルを以下のような数式で表しました。
$h_θ(x)=θ_0+θ_1x$
ロジスティック回帰の場合は上記の式を更にシグモイド関数という関数に入れます。式で表わすと以下のようになります。
$h_θ(x)=g(θ^TX)$ -①
$g_θ(z)=\frac{1}{1+e^{-z}}$
※ちみなにシグモイド関数は以下のような波形を描きます。
f:id:brains-tech:20160307110209p:plain
シグモイド関数 - Wikipedia
jp.mathworks.com

①により、データを入力すると0~1の間の値を返してきて、それがあるクラスに当てはまる確率となります。
例えば、ある特徴量 $x_0, x_1$ を入力した時、出力された値が $h_θ(X)=0.7$ だった場合、それが当該クラスである確率は70%であるといった具合です。

では、①の式が適切な値を出力してくれるようにθの値を学習していくためにはどうすればよいでしょうか。

前回と同じく、目的関数を定めた後その目的関数を最小とするような $θ$ の値を求めていきます。

ロジスティック回帰では目的関数に以下の式を使用します。
$J(θ_0, θ_1)=-\frac{1}{m}\sum_{i=1}^{m} (y^i\log(h_θ(x^i)) + (1-y^i)\log(1-h_θ(x^i))$ -②
※ $xは特徴量, yはラベル$

ここからは前回と同じです。最急降下法を使って、この $J(θ_0, θ_1)$ を最小にする $\theta$ を求めます。
ここで偏微分とかいろいろ出てくるのですが、とりあえず置いておき、つまるところ以下の式を計算していくことになります。
次の $θ_j=$ $θ_j-α\sum_{i=1}^{m} ((h_θ(x^i)-y^i)x^i_j)$ -③

データを見る

以上がかなり簡単な理論のお話です。
移行は実装に入っていきますが、まずは今回のデータを見てみましょう。
以下の様なデータを用意しました。
f:id:brains-tech:20160403195901p:plain

赤、青でプロットされた点があり、それぞれがクラス(赤: class 1, 青: class 0)を表しています。
今回はこの２つのクラスのデータを学習して、新しい点が入力された時そのデータが赤クラスである確率はどれくらいなのか、つまりどちらのクラスなのか分類したいと思います。
※実装がメインですので、シンプルなデータにしてあります。

実装！

ここまで説明してきたことをコードに落とすと以下のようになります。
コメントで説明を入れてありますが、分かりにくいところは少し動かしてみることをオススメします。
[コード]

# -*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt


def gradient_decent(theta, X, y, iter_num, alpha):
    """
    最急降下法メソッド
    与えられたiter_num分だけ最急降下法を適用してJ(θ)を最小とするthetaを求める
    :param theta: モデルのパラメータθ
    :param X: 入力データ(学習に使われるデータ)
    :param y: 入力データに対するラベル
    :param iter_num: 最急降下法の適用回数
    :param alpha: 学習率
    :return: theta, J (学習後のtheta, 学習途中での目的関数の値)
    """

    m = float(len(X))
    J = []
    for i in range(iter_num):
        # 現在のthetaでh_thetaを計算 (
        h = sigmoid(np.inner(theta, X))  # --- ①式
        # 目的関数を可視化のために計算
        J.append(1/m * np.sum((-y*np.log(h) - (1-y)*np.log(1-h))))  # --- ②式
        # 以下で, ③式を計算してthetaを更新していく
        theta1_tmp = 1/m * np.sum((h-y)*X[:, 0])
        theta2_tmp = 1/m * np.sum((h-y)*X[:, 1])
        theta3_tmp = 1/m * np.sum((h-y)*X[:, 2])
        theta = theta - alpha*np.array([theta1_tmp, theta2_tmp, theta3_tmp])
        print "theta: ", theta
    return theta, J


def sigmoid(x):
    """
    与えられたデータをsigmoid関数へ入れて計算結果を返す
    sigmoid関数
    :param x: 入力データ
    :return: sigmoid関数適用後のデータ
    """
    return 1/(1 + np.exp(-x))
    

def main():
    # サンプルデータのロード
    # 前述した散布図のデータ
    sample_data = np.load("./data.npy")

    # ラベルデータのロード
    # 赤: 1, 青: 0 として表現
    y = np.load("./label.npy")

    # 学習データの作成
    x_tmp = np.ones(len(sample_data))
    X = np.c_[x_tmp, sample_data]

    # パラメータの初期化
    # thataは任意の整数
    # alphaは学習率で予めこちらで調整した値が代入されています
    # iter_numの回数により何回最急降下法が回るかが決まります
    theta = np.array([1,1,1])
    alpha = 0.5
    iter_num = 100

    # 最急降下法を使って学習を行う
    theta, J = gradient_decent(theta, X, y, iter_num, alpha)
    plt.plot(J)
    plt.show()


if __name__ == "__main__":
    main()

[実行結果]

$ python lr.py
theta:  [ 0.8319953   0.68474584  0.81423634]
theta:  [ 0.67634834  0.38081935  0.63610057]
theta:  [ 0.53901998  0.09768263  0.47012794]
theta:  [ 0.42658727 -0.1536833   0.32244438]
theta:  [ 0.34236108 -0.36727612  0.19797763]
theta:  [ 0.28391004 -0.54596135  0.09699565]
theta:  [ 0.24549959 -0.69715343  0.01567327]
theta:  [ 0.22146887 -0.82779591 -0.05084846]
theta:  [ 0.20750074 -0.94295183 -0.10667272]
theta:  [ 0.6395262  -2.837131   -1.09546083]
...
...
...
theta:  [ 0.89455876 -3.43082868 -1.50299413]
theta:  [ 0.90015444 -3.44389053 -1.51225235]
theta:  [ 0.90569562 -3.45683949 -1.52143847]
theta:  [ 0.91118314 -3.46967788 -1.53055366]
theta:  [ 0.91661788 -3.48240793 -1.53959907]
theta:  [ 0.92200064 -3.49503179 -1.54857581]
theta:  [ 0.92733226 -3.50755156 -1.55748496]
theta:  [ 0.93261354 -3.51996927 -1.56632758]

[目的関数 $J(θ)$ の推移]
f:id:mhigu:20160314074359p:plain

学習の結果、目的関数の値が小さくなり収束していっていることが分かります。

確認のために学習したモデルに値を入力して結果を観測します。
以下にテスト用に別に作ったスクリプトを記載します。
$x_0=-2, x_1=1$ の場合

# -*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt


def sigmoid(x):
    # 与えられたデータをシグモイド関数へ入れて
    # 計算された結果を返す
    return 1/(1 + np.exp(-x))
    

def test():

    # プロット用にデータをロード
    X = np.load("data.npy")
    label = np.load("label.npy")

    # theta:  [ 0.93261354 -3.51996927 -1.56632758] を利用
    theta = np.array([0.93261354, -3.51996927, -1.56632578])

    # x0=2, x1=1を設定
    x = np.array([1, -2, 1])
    class_1_rate = sigmoid(np.dot(theta, x))
    class_num = -1

    # モデルから返ってきた値が0.5より大きければ赤(class 1)
    # 0.5以下であれば青(class 0) 
    if class_1_rate > 0.5:
        class_num = 1
    elif class_1_rate <= 0.5:
        class_num = 0
    print "[RESULT]"
    print "################"
    print "Model Output: {0}".format(class_1_rate)
    print "Input Data is class {0}".format(class_num)
    print "################"

    # 散布図をプロット
    plt.scatter(X[:,0], X[:,1], c = label, s= 100)
    plt.scatter(x[1], x[2], marker='^', color='g', s=150)
    plt.show()

if __name__ == "__main__":
    test()

[実行結果]

[RESULT]
################
Model Output: 0.998351478053
Input Data is class 1
################

ちゃんと動いているように見えます。。。
散布図を作って見るとしっかり分類できていますね。(緑の三角がinputしたものです)
f:id:brains-tech:20160404000807p:plain
ひとまずは実装完了です。
が、本来であれば、過学習を考慮して罰則項を足してモデルを作ったり、特徴量を取るときに正規化を行ったりと様々工夫が入ります。興味のある方はそのへんまで是非是非調べてみてください。

結び

今回は分類問題を解くための教師あり機械学習アルゴリズム、ロジスティック回帰を実装してきました。
分類問題を解くための機械学習のアルゴリズムには、ロジスティック回帰の他にも、SVM、ニューラルネットワークなどたくさんあります。最近ディープラーニングがすごく話題になってますが、ニューラルネットワークはその走りですよね。次回は、そのニューラルネットワークを今回のロジスティック回帰との違いに着目しながら実装していきたいと思います。

ではでは(((((((((((っ･ω･)っﾌﾞｰﾝ