フィッシャーの判別分析について-フィッシャーのあやめを例に-
来週勉強会でやる内容の一部を書きます。
実演は当日やります.....
そういえば初めてコメントをもらいました。
嬉しかったです。
フィッシャーの判別分析の一通りの流れ
1. 群内分散共分散行列W、群間分散・共分散行列Bを得る
2. を求める。
3. Xの固有値λ(1), λ(2), ..... , λ(n)とそれぞれの固有値ベクトルl^(1),l^(2), ..... , l^(n)を求める
4. 判別関数
が求まる。
フィッシャーの線形判別しました。
詳しくはPRMLのp185~190までを読んでください。
データは「フィッシャーのあやめ」という有名なデータ。
現代の推測統計学の基礎を築いた統計学者のうちの一人のフィッシャーさんが
使った”あやめ”（実はアンダーソンさんが計測したものらしい）のデータを使います。
ちなみにデータ（後で整形したものを使います。）は
1. ガクの長さ
2. ガクの幅
3. 花弁の長さ
4. 花弁の幅
5. 花の種類（これについては、群と見なします。）
と5つの特徴量を持っています。
このデータを使って、花の種類を数理的なモデルを使って分類したいと思います。
分類をする際に人間にとって分かりやすいのは作図です。
ここでも作図をしたいのですが、
この四次元のデータを使って作図することはできないので、次元を減らさなければなりません。
そこで、特徴量を(1)ようにのベクトルを四つ置きます。
ベクトルを用いて以下のような線形関数
を求めることで、
この線形関数が特徴量から得た値からサンプルがどの群に属するのか予測（分類）するために[tex; L]で”分離”できるように”を選ぶのです。
この線型関数は群間分散・共分散行列、全体分散・共分散行列、群間分散・共分散を用いて作ることができます。
ここではは群間分散・共分散行列で、
は全体分散・共分散行列となる。さらに、
は群間分散・共分散で定義されます。ちなみに、
は
と
で表すことができるので今回はこの関係を使うことにします。
そして、今回はなんやかんや省略して、は以下の行列
の固有値ベクトルの成分を用いて表現することができるのです。
すなわち、最大固有値()と、その次に大きい値の固有値(
)の固有ベクトル成分は先に述べた線型関数
の成分として
と
と扱うことができます。今回は二つですが足りない場合は、その次の大きい値の固有値の固有ベクトルを用いて判別関数を作ればいいわけです。
さて、以下はコードです。
#iris.txtは一番下にあります。
from numpy import * #データの読み込み(150,5) data=loadtxt("iris.txt") #data=data.T #Virginicaのデータ Vi=array([data[i] for i in range(0,50)]) Vi=Vi.T Vi=array([Vi[i] for i in range(1,5)]) #Versicolorのデータ Ve=array([data[i] for i in range(50,100)]) Ve=Ve.T Ve=array([Ve[i] for i in range(1,5)]) #Setonaのデータ Se=array([data[i] for i in range(100,150)]) Se=Se.T Se=array([Se[i] for i in range(1,5)]) #全体のデータ data=data.T #(4,150) data=array([data[i] for i in range(1,5)]) W=zeros((4,4)) #Seの群内分散 m1=array([[mean(i) for i in Se]]) m1=m1.T W+=dot(Se-m1,transpose(Se-m1)) #Veの群内分散 m2=array([[mean(i) for i in Ve]]) m2=m2.T W+=dot(Ve-m2,transpose(Ve-m2)) #Viの群内分散 m3=array([[mean(i) for i in Vi]]) m3=m3.T W+=dot(Vi-m3,transpose(Vi-m3)) #全体分散共分散行列 M=array([[mean(i) for i in data]]).T T=dot(data-M,transpose(data-M)) #クラス間分散共分散 B=T-W #Sw^-1*St X=dot(linalg.inv(W),B) #固有値と固有値ベクトル la,v = linalg.eig(X) #グラフ作成 from pylab import * a=array([ sum(i) for i in v.T[0]*Vi.T]) b=array([ sum(i) for i in v.T[1]*Vi.T]) plot(a,b,'y*',label='Virginica') a=array([ sum(i) for i in v.T[0]*Ve.T]) b=array([ sum(i) for i in v.T[1]*Ve.T]) plot(a,b,'b+',label='Versicolor') a=array([ sum(i) for i in v.T[0]*Se.T]) b=array([ sum(i) for i in v.T[1]*Se.T]) plot(a,b,'ro',label='Setosa') #ラベルの追加 title('Discriminant Analysis') #タイトル xlabel('L(1)') # X 軸 ylabel('L(2)') # Y 軸 # 凡例 legend() # 描画 show()
最終的に以下の図のようになります。
というわけで超雑ですが、今週の8/28日に発表後にはちゃんと書き直します。
#iris.txt
#一番左から番号、ガクの長さ、ガクの幅、花弁の長さ、花弁の幅
#さららに1~50までがセトーナ(setona)、51~100までがヴェルシカラー(versicolor)
#101~150までがヴァージニカ(virginica)です。
1 5.1 3.5 1.4 0.2 2 4.9 3.0 1.4 0.2 3 4.7 3.2 1.3 0.2 4 4.6 3.1 1.5 0.2 5 5.0 3.6 1.4 0.2 6 5.4 3.9 1.7 0.4 7 4.6 3.4 1.4 0.3 8 5.0 3.4 1.5 0.2 9 4.4 2.9 1.4 0.2 10 4.9 3.1 1.5 0.1 11 5.4 3.7 1.5 0.2 12 4.8 3.4 1.6 0.2 13 4.8 3.0 1.4 0.1 14 4.3 3.0 1.1 0.1 15 5.8 4.0 1.2 0.2 16 5.7 4.4 1.5 0.4 17 5.4 3.9 1.3 0.4 18 5.1 3.5 1.4 0.3 19 5.7 3.8 1.7 0.3 20 5.1 3.8 1.5 0.3 21 5.4 3.4 1.7 0.2 22 5.1 3.7 1.5 0.4 23 4.6 3.6 1.0 0.2 24 5.1 3.3 1.7 0.5 25 4.8 3.4 1.9 0.2 26 5.0 3.0 1.6 0.2 27 5.0 3.4 1.6 0.4 28 5.2 3.5 1.5 0.2 29 5.2 3.4 1.4 0.2 30 4.7 3.2 1.6 0.2 31 4.8 3.1 1.6 0.2 32 5.4 3.4 1.5 0.4 33 5.2 4.1 1.5 0.1 34 5.5 4.2 1.4 0.2 35 4.9 3.1 1.5 0.2 36 5.0 3.2 1.2 0.2 37 5.5 3.5 1.3 0.2 38 4.9 3.6 1.4 0.1 39 4.4 3.0 1.3 0.2 40 5.1 3.4 1.5 0.2 41 5.0 3.5 1.3 0.3 42 4.5 2.3 1.3 0.3 43 4.4 3.2 1.3 0.2 44 5.0 3.5 1.6 0.6 45 5.1 3.8 1.9 0.4 46 4.8 3.0 1.4 0.3 47 5.1 3.8 1.6 0.2 48 4.6 3.2 1.4 0.2 49 5.3 3.7 1.5 0.2 50 5.0 3.3 1.4 0.2 51 7.0 3.2 4.7 1.4 52 6.4 3.2 4.5 1.5 53 6.9 3.1 4.9 1.5 54 5.5 2.3 4.0 1.3 55 6.5 2.8 4.6 1.5 56 5.7 2.8 4.5 1.3 57 6.3 3.3 4.7 1.6 58 4.9 2.4 3.3 1.0 59 6.6 2.9 4.6 1.3 60 5.2 2.7 3.9 1.4 61 5.0 2.0 3.5 1.0 62 5.9 3.0 4.2 1.5 63 6.0 2.2 4.0 1.0 64 6.1 2.9 4.7 1.4 65 5.6 2.9 3.6 1.3 66 6.7 3.1 4.4 1.4 67 5.6 3.0 4.5 1.5 68 5.8 2.7 4.1 1.0 69 6.2 2.2 4.5 1.5 70 5.6 2.5 3.9 1.1 71 5.9 3.2 4.8 1.8 72 6.1 2.8 4.0 1.3 73 6.3 2.5 4.9 1.5 74 6.1 2.8 4.7 1.2 75 6.4 2.9 4.3 1.3 76 6.6 3.0 4.4 1.4 77 6.8 2.8 4.8 1.4 78 6.7 3.0 5.0 1.7 79 6.0 2.9 4.5 1.5 80 5.7 2.6 3.5 1.0 81 5.5 2.4 3.8 1.1 82 5.5 2.4 3.7 1.0 83 5.8 2.7 3.9 1.2 84 6.0 2.7 5.1 1.6 85 5.4 3.0 4.5 1.5 86 6.0 3.4 4.5 1.6 87 6.7 3.1 4.7 1.5 88 6.3 2.3 4.4 1.3 89 5.6 3.0 4.1 1.3 90 5.5 2.5 4.0 1.3 91 5.5 2.6 4.4 1.2 92 6.1 3.0 4.6 1.4 93 5.8 2.6 4.0 1.2 94 5.0 2.3 3.3 1.0 95 5.6 2.7 4.2 1.3 96 5.7 3.0 4.2 1.2 97 5.7 2.9 4.2 1.3 98 6.2 2.9 4.3 1.3 99 5.1 2.5 3.0 1.1 100 5.7 2.8 4.1 1.3 101 6.3 3.3 6.0 2.5 102 5.8 2.7 5.1 1.9 103 7.1 3.0 5.9 2.1 104 6.3 2.9 5.6 1.8 105 6.5 3.0 5.8 2.2 106 7.6 3.0 6.6 2.1 107 4.9 2.5 4.5 1.7 108 7.3 2.9 6.3 1.8 109 6.7 2.5 5.8 1.8 110 7.2 3.6 6.1 2.5 111 6.5 3.2 5.1 2.0 112 6.4 2.7 5.3 1.9 113 6.8 3.0 5.5 2.1 114 5.7 2.5 5.0 2.0 115 5.8 2.8 5.1 2.4 116 6.4 3.2 5.3 2.3 117 6.5 3.0 5.5 1.8 118 7.7 3.8 6.7 2.2 119 7.7 2.6 6.9 2.3 120 6.0 2.2 5.0 1.5 121 6.9 3.2 5.7 2.3 122 5.6 2.8 4.9 2.0 123 7.7 2.8 6.7 2.0 124 6.3 2.7 4.9 1.8 125 6.7 3.3 5.7 2.1 126 7.2 3.2 6.0 1.8 127 6.2 2.8 4.8 1.8 128 6.1 3.0 4.9 1.8 129 6.4 2.8 5.6 2.1 130 7.2 3.0 5.8 1.6 131 7.4 2.8 6.1 1.9 132 7.9 3.8 6.4 2.0 133 6.4 2.8 5.6 2.2 134 6.3 2.8 5.1 1.5 135 6.1 2.6 5.6 1.4 136 7.7 3.0 6.1 2.3 137 6.3 3.4 5.6 2.4 138 6.4 3.1 5.5 1.8 139 6.0 3.0 4.8 1.8 140 6.9 3.1 5.4 2.1 141 6.7 3.1 5.6 2.4 142 6.9 3.1 5.1 2.3 143 5.8 2.7 5.1 1.9 144 6.8 3.2 5.9 2.3 145 6.7 3.3 5.7 2.5 146 6.7 3.0 5.2 2.3 147 6.3 2.5 5.0 1.9 148 6.5 3.0 5.2 2.0 149 6.2 3.4 5.4 2.3 150 5.9 3.0 5.1 1.8
