要因をさぐる〜主成分分析
<課題1>[br]都市の中に区があります。[br]区のうちのいくつかの調査結果があるとします。[br]現在は、子育て、介護、防犯、健康の4つの軸の評価に[br]なっています。[br]次のデータから区の新しい評価軸をさがしてみよう。[br][table][tr][td][/td][td]子育て[/td][td]介護[/td][td]防犯[/td][td]健康[/td][/tr][tr][td]区1[/td][td]6.0[/td][td]2.7[/td][td]2.1[/td][td]6.5[/td][/tr][tr][td]区2[/td][td]6.7[/td][td]2.2[/td][td]2.5[/td][td]6.4[/td][/tr][tr][td]区3[/td][td]5.9[/td][td]3.0[/td][td]1.9[/td][td]6.5[/td][/tr][tr][td]区4[/td][td]6.6[/td][td]5.0[/td][td]2.4[/td][td]6.8[/td][/tr][tr][td]区5[/td][td]7.0[/td][td]2.6[/td][td]2.5[/td][td]6.0[/td][/tr][tr][td]区6[/td][td]8.1[/td][td]4.2[/td][td]4.5[/td][td]6.3[/td][/tr][tr][td]区7[/td][td]7.2[/td][td]2.8[/td][td]2.5[/td][td]6.2[/td][/tr][tr][td]区8[/td][td]7.0[/td][td]3.1[/td][td]2.5[/td][td]5.9[/td][/tr][tr][td]区9[/td][td]7.2[/td][td]2.5[/td][td]2.1[/td][td]5.4[/td][/tr][tr][td]区10[/td][td]5.6[/td][td]2.2[/td][td]1.8[/td][td]4.5[/td][/tr][/table]
[b][size=150][color=#0000ff]<分散を大きくする重み付けは固有ベクトルになる>[/color][/size][/b][br]評価軸4つに対する重みづけw1,w2,w3,w4を評価の値x1,x2,x3,x4にかけると、[br]新評価値yi=Σwjxij(i=1..6,j=1..4)になる。[br]だから、平均y^=1/nΣyi=Σwj(1/nxij)=Σwj(xj^),[br][b] 分散[/b]ρ=1/nΣ(yi-m)[sup]2[/sup]=1/nΣ(Σwjxij-Σwj(xj)^)[sup]2[/sup]=1/nΣ(Σwj(xij-xj^))[sup]2[/sup]=1/nΣ(ΣΣwjwk(xij-xj^)(xik-xk^))[br] =ΣΣwjwk(1/n(xij-xj^)(xik-xk^))=ΣΣwjwksjk=Σwj(Σwksjk)[br]={w1...w6}{{w1s11+w2s12+...+w6s16},{w1s21+w2s22+....w6s26}........{w1s61+w2s62+....w6s66}}[br]={w1...w6}{ sij行列}{w1....w6}[sup]t[/sup]= =[b]w[sup]t[/sup] S w[br][/b]Sが標準分散共分散行列と言われ、sij=sjiだから対称行列になる。データから求められる定数行列になる。[br]分散[b]w[sup]t[/sup] S wが最大になるように重み変数ベクトルを動かそう。wが単位ベクトルなら1-wt w=0[br][/b]ラグランジュ関数L(w,μ)=wt S w + μ(1- wt w)をwとμで偏微分すると、=0[br]2Sw-2μw=0, 1-wt w=0から、wt S w=wt(μw)=μ(wt w)=μとなる。[br]S[b]w[/b]=μ[b]wはμが行列Sの固有値となるという意味だから、行列Sの固有値の最大値をμとし[br][/b]固有ベクトルをwにすればよい。このwで重み付けすると、分散が最大になる。[br][br][b][size=150][color=#0000ff]<固有値の大きい順の固有ベクトルを重みにしたのが主成分分析>[br][/color][/size][/b]上記が主成分を1にしたときのwの求め方になる。[br]主成分を2つ以上にするには重みベクトルw1,w2、…が直交するとして、同じように計算すればよい。Sの固有値に対する固有ベクトルは直交するので、固有値の大きい方から第1主成分、第2主成分[br]、、、とすればよいね。[br]第1主成分の重みベクトルw1での評価値をx,[br]第2主成分の重みベクトルw2での評価値をyとすると、2次元にデータを分布することができる。[br]そこから、逆に主成分を意味づけしたりすることもできる。[br]固有値を比例配分することで、主成分ごとに新評価の寄与率も求められるね。[br]
[b][size=150]<julia>[/size][/b][br][b][color=#38761d]using [/color]LinearAlgebra[/b][br]A=[6.0 2.7 2.1 6.5;[br]6.7 2.2 2.5 6.4;[br]5.9 3.0 1.9 6.5;[br]6.6 5.0 2.4 6.8;[br]7.0 2.6 2.5 6.0;[br]8.1 4.2 4.5 6.3;[br]7.2 2.8 2.5 6.2;[br]7.0 3.1 2.5 5.9;[br]7.2 2.5 2.1 5.4;[br]5.6 2.2 1.8 4.5][br]At=transpose(A)[br]x1=A[:,1][br]x2=A[:,2][br]x3=A[:,3][br]x4=A[:,4][br]function ave(x)[br] return sum(x)/length(x)[br]end[br]aves=[ave(x1),ave(x2),ave(x3),ave(x4)][br]s11= sum([(x-aves[1])^2 for x in x1])/length(x1)[br]s12= sum([(x-aves[1])*(y-aves[2]) for (x,y) in zip(x1,x2)])/length(x1)[br]s13= sum([(x-aves[1])*(y-aves[3]) for (x,y) in zip(x1,x3)])/length(x1)[br]s14= sum([(x-aves[1])*(y-aves[4]) for (x,y) in zip(x1,x4)])/length(x1)[br]s1=[s11,s12,s13,s14][br]#このように定義どおりに計算を続けていけば共分散行列は求められる。[br]#地道もいいけれど。。。[br]s1[br][color=#38761d]#======================================[br]4-element Vector{Float64}:[br] 0.4981[br] 0.2120999999999999[br] 0.4145999999999999[br] 0.09950000000000005[/color]
[b][size=150]<julia>2[/size][/b][br][color=#38761d]#統計パッケージを使っていきなり共分散行列を求めよう。[br][/color][b][color=#38761d]using [/color]LinearAlgebra[br][color=#38761d]using[/color] Statistics[br][/b]A=[6.0 2.7 2.1 6.5;[br]6.7 2.2 2.5 6.4;[br]5.9 3.0 1.9 6.5;[br]6.6 5.0 2.4 6.8;[br]7.0 2.6 2.5 6.0;[br]8.1 4.2 4.5 6.3;[br]7.2 2.8 2.5 6.2;[br]7.0 3.1 2.5 5.9;[br]7.2 2.5 2.1 5.4;[br]5.6 2.2 1.8 4.5][br][color=#0000ff]#次の1行で共分散行列を求めることができる。すばらしい。[br][/color][b]S=cov(A,corrected=false) [color=#38761d]# 4列を軸にする、行は軸にしないからFalse。[/color][br][/b][color=#38761d]#=================================[br]4×4 Matrix{Float64}:[br] 0.4981 0.2121 0.4146 0.0995[br] 0.2121 0.7261 0.3086 0.2925[br] 0.4146 0.3086 0.5176 0.132[br] 0.0995 0.2925 0.132 0.4025[br][br][/color][b]eigvals(S)[br][/b][color=#38761d]#=================================[br]4-element Vector{Float64}:[br] 0.08424486156172269[br] 0.22766837981082833[br] 0.5177726308968318[br] 1.3146141277306178[br][/color][br][b]ws=-eigvecs(S)[/b][color=#38761d][br]# 第一主成分の重み(4列目)がすべて負になってしまったから正にした。[br]#=================================[br]4×4 Matrix{Float64}:[br] -0.664542 -0.151997 -0.557937 0.473273[br] -0.151736 0.536886 0.556615 0.615556[br] 0.73034 0.0235598 -0.419361 0.538688[br] 0.0443031 -0.829515 0.450581 0.326985[br][br][/color][b]x1=A*ws[:,4] [color=#134f5c]#固有値1.314..の固有ベクトルを重みにして新評価値にする。[/color][br]x2=A*ws[:,3][b] [color=#38761d]#固有値0.517..の固有ベクトルを重みにして新評価値にする。[br][/color][/b][/b][color=#38761d]# 新評価値x1,x2で2次元表示[br][b]using[/b][/color][b] PyPlot[/b][br][b]plt.scatter(x1,x2)[br][/b]n = [1,2,3,4,5,6,7,8,9,10][br][b][color=#38761d]for[/color] (i, txt) [color=#38761d]in[/color] enumerate(n)[br] plt.annotate(txt, (x1[i], x2[i]))[br][color=#38761d]end[/color][/b]
[b]<Python>[br][/b][color=#38761d]#PythonのNumpyには統計関数が入っている!すばらしい。[br][/color][b][color=#38761d]import [/color]numpy as np[br][color=#38761d]from[/color] numpy.linalg import eig [/b][br]A=np.array([br][[6.0,2.7,2.1,6.5],[br][6.7,2.2,2.5,6.4],[br][5.9,3.0,1.9,6.5],[br][6.6,5.0,2.4,6.8],[br][7.0,2.6,2.5,6.0],[br][8.1,4.2,4.5,6.3],[br][7.2,2.8,2.5,6.2],[br][7.0,3.1,2.5,5.9],[br][7.2,2.5,2.1,5.4],[br][5.6,2.2,1.8,4.5]])[br][br][b]S=np.cov(A, rowvar=False) [color=#38761d]# 4列を軸にする、行は軸にしないからFalse。[/color][br]eigvals, eigvects= eig(S)[br]eigvals[br][/b][color=#38761d]#===================================[br][/color]array([1.46068236, 0.57530292, 0.0936054 , 0.25296487])[br][br][color=#38761d]#固有値が大きい順にならんでくれている。[br][/color][color=#38761d]# 新評価値x1,x2で2次元表示[br][b]import [/b][/color][b]matplotlib.pyplot as plt[br][br]x1=A@eigvecs[:,0][br]x2=A@eigvecs[:,1][br]plt.scatter(x1,x2)[br][/b]n = [1,2,3,4,5,6,7,8,9,10][color=#38761d][br][b]for[/b][/color][b] (i, txt) [color=#38761d]in[/color] enumerate(n):[br] plt.annotate(txt, (x1[i], x2[i]))[br][/b][br]# 結果のグラフはpythonもjuliaと同じ[br]
[b]<geogebra>[br][/b]A={{6.0,2.7,2.1,6.5},[br]{6.7,2.2,2.5,6.4},[br]{5.9,3.0,1.9,6.5},[br]{6.6,5.0,2.4,6.8},[br]{7.0,2.6,2.5,6.0},[br]{8.1,4.2,4.5,6.3},[br]{7.2,2.8,2.5,6.2},[br]{7.0,3.1,2.5,5.9},[br]{7.2,2.5,2.1,5.4},[br]{5.6,2.2,1.8,4.5}}[br][color=#38761d]# 4行のリストとして扱うために転置する。[br][/color]B=Transpose(A)[br][color=#38761d]# geogebraには共分散行列を求める関数が見当たらないのd,[br]# 4つのリストを2重のfor文のようにネストさせて共分散関数を連続して使う。[br][/color]S=Sequence(Sequence(Covariance(B(i),B(j)),j,1,4),i,1,4)[br]# ここで、「表示メニュー」から「CAS」を選び入力する。[br]# 変数:=式 の形の代入になるので注意!!![br][b]Eigenvalues(S)[br]vecs:=Transpose(Eigenvectors(S)) [color=#38761d]#要素が行指定になるため転置する。[br][/color]w1:=Element(vecs,4)[br]w2:=Element(vecs,3)[br]x1:= A w1[br]x2:= A w2 [br][color=#38761d]# ここまで来たら、「CAS」は非表示にできる。[br][/color][/b][color=#38761d]# 新評価値x1,x2で2次元表示[br][/color][b]zip((x,y),x,x1,y,y1)[/b]