기본 콘텐츠로 건너뛰기

[matplotlib] 등고선(Contour)

[data analysis] 직교행렬과 주성분

직교행렬과 주성분(Principal Component)

식 1의 행렬 X를 평균-편차 형태로 가정합니다.

(식 1)X=[x1x2xn]

주성분분석의 목적은 식 2와 같이 X=PY 형태로 변환가능한 p×p차원의 직교 행렬 P를 발견하는 것입니다.

(식 2)[x1x2xp]=[u1u2up][y1y2yp]P=[u1u2up]

식 2에서 새로운 변수 y1,y2,,yp는 상관성이 없고 내림차순으로 정렬된 형태입니다. Y는 P에 대해 X의 좌표벡터가 됩니다. P가 직교행렬이므로 식 3이 성립합니다(정규직교의 특징 참조).

(식 3)P1PY=PTPY=PTXY=PTX

식 3으로 부터 새로운 변수인 Y 역시 평균-편차 형태이므로 Y의 공분산 행렬을 유도 할 수 있습니다(식 4).

(식 4)YYT=(PTX)(PTX)T=PTXXTP=PTSP

S는 변수 X에 대한 공분산 행렬로서 대칭행렬입니다. 식 4에서 YYT=A라 하고 S에 대해 정리하면 식 5와 같습니다.

(식 5)S=PAPT

S가 대칭행렬이고 P가 직교행렬이므로 위 식은 스펙트럴 분해와 동일한 형태가 됩니다. 그러므로 행렬 A=YYT는 행렬 S의 고유값들을 내림차순으로 정렬한 값들을 대각요소로 하는 대각행렬이 됩니다. 결과적으로 새로운 변수의 공분산 행렬은 원래의 변수의 공분산 행렬의 고유값에 의한 대각행렬이 됩니다.

공분산 행렬 S의 고유값에 대응하는 단위 고유벡터들 u1,u2,,up은 자료의 주성분(principal component)라고 합니다. 고유값들을 내림차순으로 정렬하기 때문에 제1 주성분은 S의 가장 큰 고유값에 대응하는 고유벡터, 두번째 주성분은 두번째로 큰 고유값에 대응하는 고유벡터가 됩니다. 즉, P는 S의 단위고유벡터들을 성분으로 하는 고유행렬이 되므로 제1 주성분에 의한 새로운 변수는 식 6과 같이 나타낼 수 있습니다.

(식 6)y1=u1TX

import numpy as np
import numpy.linalg as la
import pandas as pd
from sympy import *
np.set_printoptions(precision=4, suppress=True)

예 1)

어떤 데이터의 공분산 행렬 S가 다음과 같다면 그 자료의 주성분?

S=[2382.782611.842136.22611.843160.472553.92136.22553.92650.71]

S=np.array([[2382.78, 2611.84, 2136.20],[2611.84, 3160.47, 2553.90],[2136.20, 2553.90, 2650.71]])
print(S)
[[2382.78 2611.84 2136.2 ]
 [2611.84 3160.47 2553.9 ]
 [2136.2  2553.9  2650.71]]

S를 식 5에 대입하기 위해 고유행렬과 고유벡터를 계산하기 위해 np.eig() 함수를 적용합니다.

d, p = la.eig(S)
print(d.round(2))
[7635.72  125.24  433.  ]
print(p.round(2))
[[-0.54 -0.7  -0.46]
 [-0.63  0.7  -0.33]
 [-0.56 -0.11  0.82]]

식 5의 A 즉, 새로운 변수들의 공분산 행렬(YYT)은 다음의 결과와 같이 S의 고유값을 대각요소로 하는 대각행렬입니다.

D=np.diag(d)  #새로운 변수의 공분산 행렬
print(D.round(2))
[[7635.72    0.      0.  ]
 [   0.    125.24    0.  ]
 [   0.      0.    433.  ]]

S의 고유값 중 가장 큰 7635.72에 대응하는 고유벡터 p[:,0]이 제1 주성분이 됩니다. 즉, 행렬 p의 각 열이 주성분이 됩니다.

예 2)

다음은 5명의 남아의 몸무게와 키에 대한 자료입니다.

0 1 2 3 4
Weight120125125135145
Height61 60646872

자료 내에서 대부분의 변화를 설명하는 새로운 단일한 변수를 발견하기 위한 주성분 분석?

data=np.array([[120, 125, 125, 135, 145], [61, 60, 64, 68, 72]])
print(data)
[[120 125 125 135 145]
 [ 61  60  64  68  72]]
S=np.cov(data, rowvar=True) # 위 data의 공분산 행렬
print(S)
[[100.   47.5]
 [ 47.5  25. ]]
d,P=la.eig(S)
print(d.round(2))
[123.02   1.98]
print(P.round(2)) # data 공분산행렬에 대한 주성분
[[ 0.9  -0.44]
 [ 0.44  0.9 ]]

제1 주성분에 의한 새로운 단일한 변수 y1은 다음과 같이 결정됩니다.

Y=u1TXu1=P[:,0]

weight, height=symbols("w h")
x=Matrix(2,1,[weight, height]); x

[wh]

pc1=Matrix(P[:,0].round(2))
pc1

[0.90.44]

y1=pc1.T*x
y1

[0.44h+0.9w]

주성분에 의해 새로운 변수를 생성할 수 있습니다. 즉, weight, height 두 변수는 하나의 새로운 변수 y로 대체 될 수 있습니다. 결과적으로 자료의 차원이 축소됩니다. 위의 분석은 python의 모듈 sklearn.decomposition.PCA() 클래스로 실행될 수 있습니다.

  • sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False, svd_solver=’auto’, tol=0.0, iterated_power=’auto’, random_state=None)
    • n_components : 산출하고자 하는 주성분의 수로 정수 입니다. 이 값을 정하지 않는다면 모든 주성분은 다음과 같이 지정됩니다. n_components=min(n_sample, n_features), 결국 원 자료의 공분산행렬의 모든 고유값에 대응되는 모든 고유값을 사용합니다.

PCA() 클래스를 예 2에 적용합니다.

from sklearn.decomposition import PCA
data1=data.T
print(data1)
[[120  61]
 [125  60]
 [125  64]
 [135  68]
 [145  72]]
pca=PCA() #PCA 클래스(pca) 생성
pca.fit(data1) #pca 모형생성
pc=pca.components_#고유벡터, 주성분 벡터
print(pc)
[[ 0.89990119  0.43609386]
 [ 0.43609386 -0.89990119]]

n_components를 지정하지 않았으므로 모든 주성분을 고려합니다.

eigvalY=pca.explained_variance_#고유값
print(eigvalY)
[123.01859218   1.98140782]

위 결과에 의한 대각행렬은 weight와 height의 공분산 행렬(XXT)을 기반으로 하는 새로운 변수(Y)의 공분산행렬이 됩니다.

X=PYP1X=PTY=X(식 7)YYT=PTX(PTX)T=PTXXTP=PTSPS=P(YYT)PT

D=np.diag(eigvalY)
print(D)
[[123.01859218   0.        ]
 [  0.           1.98140782]]

위 결과 D는 새로운 변수 Y의 공분산행렬로 Y를 구성하는 여러 벡터들의 분산만 존재하며 그들사이의 공분산은 0입니다. 이 모델의 속성 .explained_variance_ratio_을 사용하여 각각의 주성분이 전체 분산에 영향을 주는 비율을 확인할 수 있습니다.

위 결과는 S=PDPT=P(YYT)PT에서 새로운 변수 Y 의 공분산 행렬로 각 주성분에 대한
분산을 나타냅니다. (D는 대각 행렬외에 모두 0이므로 공분산은 존재하지 않습니다. 즉, 원래의 자료를 직교행렬을 사용하여 변환하는 것으로서 변수들 간의 상관성은 0이 됩니다.)
pca.explained_variance_ratio_
array([0.98414874, 0.01585126])

위의 경우 제 1 주성분에 의한 분산은 전체의 약 98%를 차지합니다.

댓글

이 블로그의 인기 게시물

[Linear Algebra] 유사변환(Similarity transformation)

유사변환(Similarity transformation) n×n 차원의 정방 행렬 A, B 그리고 가역 행렬 P 사이에 식 1의 관계가 성립하면 행렬 A와 B는 유사행렬(similarity matrix)이 되며 행렬 A를 가역행렬 P와 B로 분해하는 것을 유사 변환(similarity transformation) 이라고 합니다. (1)A=PBP1P1AP=B 식 2는 식 1의 양변에 B의 고유값을 고려한 것입니다. (식 2)BλI=P1APλP1P=P1(APλP)=P1(AλI)P 식 2의 행렬식은 식 3과 같이 정리됩니다. det(BλI)=det(P1(APλP))=det(P1)det((AλI))det(P)=det(P1)det(P)det((AλI))=det(AλI)det(P1)det(P)=det(P1P)=det(I) 유사행렬의 특성 유사행렬인 두 정방행렬 A와 B는 'A ~ B' 와 같...

[sympy] Sympy객체의 표현을 위한 함수들

Sympy객체의 표현을 위한 함수들 General simplify(x): 식 x(sympy 객체)를 간단히 정리 합니다. import numpy as np from sympy import * x=symbols("x") a=sin(x)**2+cos(x)**2 a sin2(x)+cos2(x) simplify(a) 1 simplify(b) x3+x2x1x2+2x+1 simplify(b) x - 1 c=gamma(x)/gamma(x-2) c Γ(x)Γ(x2) simplify(c) (x2)(x1) 위의 예들 중 객체 c의 감마함수(gamma(x))는 확률분포 등 여러 부분에서 사용되는 표현식으로 다음과 같이 정의 됩니다. 감마함수는 음이 아닌 정수를 제외한 모든 수에서 정의됩니다. 식 1과 같이 자연수에서 감마함수는 factorial(!), 부동소수(양의 실수)인 경우 적분을 적용하여 계산합니다. (식 1)Γ(n)={(n1)!n:자연수0xn1exdxn:부동소수 x=symbols('x') gamma(x).subs(x,4) 6 factorial 계산은 math.factorial() 함수를 사용할 수 있습니다. import math math.factorial(3) 6 a=gamma(x).subs(x,4.5) a.evalf(3) 11.6 simpilfy() 함수의 알고리즘은 식에서 공통사항을 찾아 정리하...

sympy.solvers로 방정식해 구하기

sympy.solvers로 방정식해 구하기 대수 방정식을 해를 계산하기 위해 다음 함수를 사용합니다. sympy.solvers.solve(f, *symbols, **flags) f=0, 즉 동차방정식에 대해 지정한 변수의 해를 계산 f : 식 또는 함수 symbols: 식의 해를 계산하기 위한 변수, 변수가 하나인 경우는 생략가능(자동으로 인식) flags: 계산 또는 결과의 방식을 지정하기 위한 인수들 dict=True: {x:3, y:1}같이 사전형식, 기본값 = False set=True :{(x,3),(y,1)}같이 집합형식, 기본값 = False ratioal=True : 실수를 유리수로 반환, 기본값 = False positive=True: 해들 중에 양수만을 반환, 기본값 = False 예 x2=1의 해를 결정합니다. solve() 함수에 적용하기 위해서는 다음과 같이 식의 한쪽이 0이 되는 형태인 동차식으로 구성되어야 합니다. x21=0 import numpy as np from sympy import * x = symbols('x') solve(x**2-1, x) [-1, 1] 위 식은 계산 과정은 다음과 같습니다. x21=0(x+1)(x1)=0x=1or1x4=1의 해를 결정합니다. solve() 함수의 인수 set=True를 지정하였으므로 결과는 집합(set)형으로 반환됩니다. eq=x**4-1 solve(eq, set=True) ([x], {(-1,), (-I,), (1,), (I,)}) 위의 경우 I는 복소수입니다.즉 위 결과의 과정은 다음과 같습니다. x41=(x2+1)(x+1)(x1)=0x=±1,±1=±i,±1 실수...