ピマ・インディアンデータを用いた識別関数

データの準備:散布図を作成しよう!

ピマ・インディアンのデータを利用します。RのMASSパッケージを用いると、200名のピマ・インディアンの女性の糖尿病診断に関する学習データ(Pima.tr)と、332名のテストデータ(Pima.te)が簡単に取得できるため、Rからデータを取得し、そちらをPythonで利用していきます。こちらはRです:

library(MASS)
write.csv(Pima.tr,"pimatr.csv")

次にPythonです。ピマ・インディアンのデータの、gluとbmiに注目していきましょう。糖尿病になった人とそうでない人を色別けし、プロットするにはどうすればよいでしょうか。

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

pimatr = np.array(pd.read_csv('pimatr.csv'))
pimatr = np.array(pimatr)
glu = pimatr[:,2]
bmi = pimatr[:,5]
yesno = (pimatr[:,8]=='Yes')
plt.plot(glu[yesno],bmi[yesno],'o',color="r",label='diabetes')
plt.plot(glu[~yesno],bmi[~yesno],'o',color="b",label='healthy')
plt.xlim([50,210])
plt.ylim([15,50])
plt.xlabel('glu')
plt.ylabel('bmi')
plt.title('Scatter Plot of Data')
plt.legend(loc='lower right',numpoints=1) #numpointsをつけると、マーカーが2つではなく1つだけ表示されます。

f:id:decompose:20160326234719p:plain

統計量を計算しよう!

各種統計量を計算していきましょう。

#class別にデータを準備
class2 = np.array([glu[yesno],bmi[yesno]])
class1 = np.array([glu[~yesno],bmi[~yesno]])
class0 = np.array([glu,bmi])
#各クラスの平均ベクトル
mu1 = np.array([np.mean(class1[0,:]),np.mean(class1[1,:])])
mu2 = np.array([np.mean(class2[0,:]),np.mean(class2[1,:])])
#クラスの事前確率
P1 = len(class1[0,:])*1.0 / ( len(class1[0,:]) + len(class2[0,:]) )
P2 = len(class2[0,:])*1.0 / ( len(class1[0,:]) + len(class2[0,:]) )
#各クラスの分散共分散行列
cov1 = np.cov(class1)
cov2 = np.cov(class2) 

2次識別関数を図示しよう!

ここでは正規分布から導かれる識別関数を考えています。まずは、識別関数に出てくるパーツを予め準備します。

S = np.linalg.inv(cov1) - np.linalg.inv(cov2)
c = mu2.dot( np.linalg.inv(cov2)) - mu1.dot(np.linalg.inv(cov1))
F = mu1.transpose().dot(np.linalg.inv(cov1)).dot(mu1.transpose()) - mu2.dot(np.linalg.inv(cov2)).dot(mu2.transpose()) +np.log(np.linalg.det(cov1)/np.linalg.det(cov2)) -2*np.log(P1/P2)

次に識別関数値を計算します。考え方としては、bluとbmiの値を色々動かしてみた時の関数{f}の値を保存していきます。

x = np.arange(50,210,1)
y = np.arange(15,50,1)
X,Y=np.meshgrid(x,y)
#Xnew Ynewは、ベクトルに直しています。
Xnew = X.reshape(1,35*160)
Ynew = Y.reshape(1,35*160)
#XYの各列に、XY平面上の点が格納されているイメージです。
XY = np.array([Xnew[0,:],Ynew[0,:]])

最後に、識別関数値をそれぞれのXY点に対して計算します。XY点を{x}とすると式は{f(x)=x^T S x+2c^T x +F}です。これがゼロになるところが識別境界です。

f = np.zeros(len(Xnew[0,:])) #fの箱を作っておきます。
for i in range(len(Xnew[0,:])):
    f[i] = XY[:,i].transpose().dot(S).dot(XY[:,i]) + 2 * c.dot(XY[:,i]) + F
Z = np.reshape(f,(35,160))

f:id:decompose:20160329185425p:plain

線形識別関数を図示しよう!

線形識別関数の場合は、2つのクラスの共分散行列は同じと仮定して求めます。そこで、2次識別関数の場合のcov1とcov2の部分を下記のように置き換えます。

covall = P1*np.cov(class1) + P2*np.cov(class2)
cov1 = covall
cov2 = covall

f:id:decompose:20160329185415p:plain

参考

こちらの4章確率モデルと識別関数の中の、4.2確率モデルに基づいています。
www.amazon.co.jp