ロジスティック回帰
ロジスティック回帰とは
今、観測されたデータ x が2つのクラス{C1,C2}のいずれかに属するとします。 データ x がクラスC1に属する確率(事後確率)は、ベイズの定理を用いて、次のように書くことができます。
(1)
データ x がクラスC2に属する確率は、p(C2|x)=1-p(C1|x) となります。2つのクラス{C1,C2}に属するデータが等しい共分散行列を持つガウス分布に従うと仮定します。すなわち、
(2)
(2)式の条件付き確率は、一般に尤度と呼ばれます。これを用いて、(1)式の事後確率は、次式のように書けます[1]。
(3)
ここで、次式はシグモイド関数と呼ばれます。
(4)
(3)式の回帰式の係数は次式で与えられます。
(5) (6)
図1は、尤度p(x|C1) (青)、p(x|C2) (緑)に相当するヒスとグラムと、これに対応する事後確率p(C1|x) (青)、p(C2|x) (緑)の例です。実際の応用では、学習サンプルから尤度のパラメタ(μ1, μ2, Σ)を学習し、回帰係数(5)、(6)を求めておきます。クラスが未知のデータ x が入ってきたら、事後確率p(C1|x) (青)、p(C2|x) を計算して、x がどちらのクラスに属するかの確率値を求めます。
図1:尤度(ヒストグラム)と事後確率(曲線)の例
参考文献
[1] C. M. Bishop, "Pattern recognition and machine learning," Springer, 2006
サンプルプログラム
logisticRegression2D.py
このサンプルでは、sklearnライブラリのlinear_model.LogisticRegression()を用いて、2次元のデータについて、上述のロジスティック回帰の学習と予測を行う例題を示します。
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 8 18:43:56 2016
@author: asano
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from sklearn import linear_model
from mpl_toolkits.mplot3d import Axes3D
FigSize = (8,6)
(gxmin,gxmax) = (-4,10)
(gymin,gymax) = (-4,10)
# データの生成
mu1 = np.array([1.0, 1.0])
cov1 = np.array([[1.0,0.0],[0.0, 1.0]])
n1 = 1000
x1 = np.random.multivariate_normal(mu1, cov1, n1)
mu2 = np.array([4.0, 3.0])
cov2 = np.array([[1.5,0.0],[0.0, 2.0]])
n2 = 1000
x2= np.random.multivariate_normal(mu2, cov2, n2)
# 散布図
plt.figure(figsize=FigSize)
plt.scatter(x1[:,0], x1[:,1],c='r')
plt.scatter(x2[:,0], x2[:,1],c='b')
plt.xlabel('x1')
plt.ylabel('x2')
plt.xlim(gxmin,gxmax)
plt.ylim(gymin,gymax)
# 学習
X = np.concatenate((x1,x2),axis=0)
y = np.ones(n1+n2)
y[n1:n1+n2] = -1
clf = linear_model.LogisticRegression(C=1e5)
clf.fit(X, y)
score = clf.score(X, y)
yHat = clf.predict(X)
# 識別結果の表示
plt.figure(figsize=FigSize)
plt.plot(range(n1+n2),y,'b')
plt.plot(range(n1+n2),yHat,'ro')
plt.ylim(-1.2,1.2)
plt.title('Accuracy: ' + str(score) )
# 確率密度関数
def sigmoid(x):
return 1 / (1 + np.exp(-x))
(ngx,ngy)=(100,100)
gx = np.linspace(gxmin,gxmax,ngx)
gy = np.linspace(gymin,gymax,ngy)
p1 = np.zeros((ngx,ngy))
for m in range(ngx):
for n in range(ngy):
x = np.array([gx[m], gy[n]])
p1[m,n] = sigmoid( np.dot(x, clf.coef_[0]) + clf.intercept_)
plt.figure(figsize=FigSize)
(gxm,gym) = np.meshgrid(gx,gy)
plt.pcolor(gxm, gym, p1)
plt.xlabel('x1')
plt.ylabel('x2')
fig = plt.figure(figsize=FigSize)
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(gxm, gym, p1, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.view_init(30, 20)
plt.xlabel('x1')
plt.ylabel('x2')
fig.colorbar(surf, shrink=0.5, aspect=5)
図2:データの散布図
図3:識別結果と識別率
図4:確率密度関数 p(C1|x)
図5:確率密度関数 p(C1|x) (3次元表示)