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