ロジスティック回帰を実装してみる

ロジスティック回帰

2クラス分類のロジスティック回帰をPythonで実装して見たいと思います。

C1, C2の2クラス分類についてクラスC1の事後確率は以下のように表せられます。

この時P(C1|x)は何かしらの事象xが発生した後にそれがC1である確率になります。次にp(x|C1)はクラスC1に分類されるという条件付きででxの事象に含まれる確率です。それからP(C1)はC1に分類される確率です。p(x|C2), P(C2)についても同様です。

この時以下のように置いたら

事後確率は次のように表すことができます

このσ(a)のことをシグモイド関数と呼びます。aが-∞ ~ ∞の範囲の値をとる時シグモイド関数に渡した結果は0 ~ 1の値を取るので確率を表すモデルとして使用できます。

ロジスティック回帰では一般化線形モデルの結果をシグモイド関数に渡して使用します。

例えばアヤメの品種がsetosaかvirginicaかをSepalLengthCmとSepalWidthCmでプロットすると以下のようになるので

シグモイド関数に渡す線形モデルはb + w1x1 + w2x2のように表すことができます。 シグモイド関数の結果が0.8であればC1に分類される確率は0.8でC2に分類される確率は0.2と考えることができます。

線形モデルであればあとは誤差関数から勾配を求めてパラメータ更新すれば良いのですが、ロジスティック回帰でも同様に勾配を求めてパラメータ更新を繰り返して正しいモデルを手に入れれるようにします。

ロジスティック回帰の誤差関数には尤度関数を使います。N回の試行に基づく尤度関数は以下のように表せます。

このπiはi番目の事象がクラスC1に分類される確率で、tiはi番目の事象がC1なら1でC2なら0とすることで尤度関数は分類が外れる確率の結果を掛け合わせたものと考えることができ、尤度関数の結果が小さい程よいモデルだと考えることができます。

ただプログラムの計算で尤度関数をそのまま扱うと、0より小さい数をなんども掛け合わせることになるので正しい結果を求められなくなります。なのでプログラム中でも利用できるように負の対数で表します。

この時piは以下のように表すことができるので

これを負の対数尤度関数に代入する

ここまで求めたらあとは勾配を求めると以下のようになります。

勾配の計算方法がわかったのであとはpythonで実装してみます。
アヤメのcsvデータを読み込んで負の対数尤度が小さくなるかを確認します。

%matplotlib inline
import math
import pandas as pd
import warnings
warnings.filterwarnings("ignore")
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
sns.set(style="white", color_codes=True)

iris = pd.read_csv("./Iris.csv") # the iris dataset is now a Pandas DataFrame


test_data1 = pd.DataFrame(index=iris_testdata.index, columns=[])
test_data1['bias'] = 1
test_data1['SepalLengthCm'] = iris_testdata['SepalLengthCm']
test_data1['SepalWidthCm'] = iris_testdata['SepalWidthCm']
test_data1.head(3)

ans_data1 = pd.DataFrame(index=iris_testdata.index, columns=[])
ans_data1['ans'] = iris_testdata[['ans']]
ans_data1.head(3)

learn_rate = 0.3
current_model = np.array([0.0, 0.0, 0.0])
calc = pd.DataFrame(index=iris_testdata.index, columns=[])
calc['y'] = np.dot(test_data1, current_model)
calc['pi'] = math.e ** calc['y'] / ( 1 + math.e ** calc['y'])
calc.tail(3)

for i in range(300):
    calc = pd.DataFrame(index=iris_testdata.index, columns=[])
    calc['y'] = np.dot(test_data1, current_model)
    calc['pi'] = math.e ** calc['y'] / ( 1 + math.e ** calc['y'])
    calc['sig'] = 1 / ( 1 + math.e ** calc['y'])
    calc['udo'] = - ans_data1['ans']  * list(map(lambda pi:math.log(pi), calc['pi']) ) - (1 - ans_data1['ans'])  * list(map(lambda pi:math.log(1 - pi), calc['pi']) )
    grad_base = (calc['pi'] - ans_data1['ans'])
    sepLenGrad = (test_data1['SepalLengthCm'] * grad_base/len(test_data1)).sum()
    sepWidGrad = (test_data1['SepalWidthCm'] * grad_base/len(test_data1)).sum()
    current_model[1] -= sepLenGrad * learn_rate
    current_model[2] -= sepWidGrad * learn_rate
    calc['y'] = np.dot(test_data1, current_model)
    calc['pi'] = math.e ** calc['y'] / ( 1 + math.e ** calc['y'])
    current_model[0] -= (calc['pi'] - ans_data1['ans']).sum() /len(test_data1) * learn_rate
print('対数尤度:' + str(calc['udo'].sum()))
current_model

対数尤度:7.942095114838584
array([-1.12827656,  2.91918166, -4.85369212])

最急降下法でパラメータ更新をするとこのような実装になるかと思います。