[data analysis] 주성분분석(principal component analysis)

주성분분석

관련된 내용

공분산행렬의 고유벡터와 그에 대응하는 고유값을 고려하여 변수행렬의 차원을 감소시키면서 지정된 주축(고유벡터)을 기반으로 존재하는 축소된 차원의 변환 데이터를 생성할 수 있습니다. 이러한 주축 즉, 고유벡터들을 주성분(principal component)라고 하며 내림차순으로 정렬된 각각의 고유값에 대응하는 고유벡터를 제1 주성분 성분, 제 2 주성분 등과 같이 구별합니다.

주성분 분석의 목적은 식 1과 같이 공분산 행렬의 고유벡터(p)와의 내적으로 변환된 축소된 차원의 벡터 Xp = Y를 결정하는 것입니다.

\begin{align}A&=y^Ty\\\tag{식 1}&=(Xp)^T(Xp)\\&=p^TX^TXp\\&=p^T\Sigma P \end{align}

import numpy as np
import pandas as pd
from scipy import stats
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
import yfinance as yf
from sklearn import feature_selection
import matplotlib.pyplot as plt
plt.rcParams['font.family'] ='NanumGothic'
plt.rcParams['axes.unicode_minus'] =False

예 1)

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

	0	1	2
0	2382.78	2611.84	2136.20
1	2611.84	3160.47	2553.90
2	2136.20	2553.90	2650.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]]

d, p=la.eig(S)
print(d.round(2))

[7635.72  125.24  433.  ]

d_idx=np.argsort(d)[::-1]
print(d_idx)

[0 2 1]

위 내림차순으로 정렬된 순서의 고유값에 대응하는 고유벡터가 각각 제1, 2, 3의 주성분이 됩니다.

pcom1=p[:,0]
print(pcom1.round(2))

[-0.54 -0.63 -0.56]

pcom2=p[:,2]
print(pcom2.round(2))

[-0.46 -0.33  0.82]

pcom3=p[:,1]
print(pcom3.round(2))

[-0.7   0.7  -0.11])

예 2)

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

	0	1	2	3	4
Weight	120	125	125	135	145
Height	61	60	64	68	72

이 자료의 주성분을 결정합니다.

X=np.array([[120, 125, 125, 135, 145], [61, 60, 64, 68, 72]])
print(X)

[[120 125 125 135 145]
 [ 61  60  64  68  72]]

#공분산 행렬
covM=np.cov(X)
print(covM)

[[100.   47.5]
 [ 47.5  25. ]]

제 1 주성분에 의한 변환된 변수벡터 Y=Xp₁. p₁는 가장 큰 고유값에 대응하는 고유벡터입니다.

d, P=la.eig(covM)
d_SortIdx=np.argsort(d)[::-1]
print(d_SortIdx)

[0 1]

print(d.round(2))

[123.02   1.98]

print(P.round(2))

[[ 0.9  -0.44]
 [ 0.44  0.9 ]]

X는 2×5 형태의 행렬이며 제1 주성분인 P[0,:]은 2차원 벡터입니다. 이들의 행렬곱을 위해서는 다음 코드와 같이 전치된 주성분(1×2)과 설명변수행렬인 X의 순서로 연산해야 합니다.

y1=P[:,0].reshape(1,-1)@X
print(y1.round(2))

[[134.59 138.65 140.4  151.14 161.88]]

위 결과와 같이 주성분에 의해 1행만 존재하는 새로운 변수를 생성할 수 있습니다. 즉, weight, height 두 변수로 구성된 2차원인 행렬 형태인 설명변수는 벡터(1차원)인 새로운 변수 y로 대체 될 수 있습니다. 결과적으로 자료는 2차원에서 1차원으로 차원이 축소됩니다. 위의 분석은 python의 모듈 sklearn의 하부모듈 decomposition의 PCA() 클래스로 실행될 수 있습니다. 통계등의 모든 데이터 분석에서 변수의 값들은 열(column) 기준으로 정리됩니다. 그러므로 이 예의 자료 X는 행기준으로 작성된 것으로 이 함수에 적용하기 위해서는 X를 전치하여야 합니다.

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

from sklearn.decomposition import PCA

X1=X.T
pca=PCA().fit(X1)
print(pca.components_.round(2)) #고유벡터(=d)

[[ 0.9   0.44]
 [ 0.44 -0.9 ]]

print(pca.explained_variance_.round(2)) #고유값

[123.02   1.98]

전체 분산 중 주성분에 의한 분산의 비는 속성 .explained_variance_ratio_로 확인할 수 있으며 이 자료에서 제1 주성분에 의한 분산은 전체의 약 98%가 됩니다.

print(pca.explained_variance_ratio_.round(2))

[0.98 0.02]

#이 모형에 사용된 평균
print(pca.mean_.round(2))

[130.  65.]

예 3)

일정기간의 환율, 다우지수, 나스닥 지수, 필라델피아 반도체 지수, 변동성 지수인 vix, 코스피 지수 자료를 사용하여 KODEX 코스피의 종가를 추정하는회귀모델을 작성합니다. 이 과정에서 설명변수의 차원을 사용된 주성분 수에 의해 축소된 설명변수 기반의 모델의 MSE와 BIC의 변화를 나타내봅니다.

모델 생성을 위한 자료는 다음 코드로 호출합니다.

st=pd.Timestamp(2021,5, 10)
et=pd.Timestamp(2025, 9, 30)
nme={'wd':'KRW=X','dj':"^DJI",'nd':'^IXIC','sox':"^sox",'vix':'^VIX','kos':'^KS11', 'k':'226490.KS'}
da=pd.DataFrame()
nme1=[]
for i, j in nme.items():
    d=yf.download(j,st, et)[["Open","High","Low","Close","Volume"]]
    nme1.append([i+k for k in d.columns])
    da=pd.concat([da, d], axis=1)
da=da.ffill()
da=da.dropna()
da.columns=np.ravel(nme1)
da.shape

(892, 35)

ind=da.drop(columns='kClose').iloc[:-1,:]
de=pd.DataFrame(da['kClose'][1:])
final=da.drop(columns='kClose').iloc[-1:,:]
ind.shape, de.shape, final.shape

((891, 34), (891, 1), (1, 34))

설명변수를 표준화합니다.

scaler=StandardScaler().fit(ind)
X=scaler.transform(ind)
Xfinal=scaler.transform(final)
Xfinal.shape

(1, 34)

자료를 모델 설정에 참여하는 훈련(train)과 검증에 사용할 test 세트로 분리합니다.

Xtr,Xte, ytr, yte=train_test_split(X, de, test_size=0.2, random_state=3)
[i.shape for i in [Xtr, ytr, Xte, yte]]

[(712, 34), (712, 1), (179, 34), (179, 1)]

sklearn.linear_model.LinearRegression() 클래스를 사용한 선형모델에서 훈련과 검정 데이터에서의 R², MSE, BIC 등의 통계량을 계산하고 비교합니다.

m_ols=LinearRegression().fit(Xtr, ytr)
pre_tr=m_ols.predict(Xtr)
mse_tr=mean_squared_error(ytr, pre_tr)
r2_tr=m_ols.score(Xtr, ytr)
print("train set_mse: %.3f, _R2: %.3f" %(mse_tr, r2_tr))

train set_mse: 57133.233, _R2: 0.991

pre_te=m_ols.predict(Xte)
mse_te=mean_squared_error(yte, pre_te)
r2_te=m_ols.score(Xte, yte)
print("teain set_mse: %.3f, _R2: %.3f" %(mse_te, r2_te))

teain set_mse: 66570.592, _R2: 0.991

모델의 통계량 BIC는 식 2에 의해 계산할 수 있습니다(최대우도 추정: AIC와 BIC 참조).

\begin{align}\tag{식 2}\text{BIC}&=n\cdot \ln(\text{MSE})+\ln(n)\cdot k\\\text{MSE}:&\, \text{mean squared error}\\n,\, k:&\, text{sample수, feature 수}\end{align}

shape=Xtr.shape
bic=shape[0]*np.log(mse_tr)+np.log(shape[0])*shape[1]
bic.round(2)

8021.95

다음은 변수 final의 예측값과 실측값입니다.

pred=m_ols.predict(Xfinal) #예측값
pred

array([[26595.49593575]])

da.tail(1)["kClose"]#관찰값

Date
2024-10-08    26530.0
Freq: B, Name: kClose, dtype: float64

위 결과는 훈련과 검정 데이터의 MSE 차이와 모델의 BIC를 나타냅니다. 이 값들을 기준으로 주성분의 수에 따라 변환되는 설명변수에 따른 MSI와 BIC를 비교해 봅니다.

R2={}
mse={}
bic={}
Pre={}
for i in [3, 5, 10, 15, 20, 25, 30, 34]:
    pca=PCA(n_components=i, random_state=3) .fit(Xtr)
    Xtr_pc=pca.transform(Xtr)
    Xte_pc=pca.transform(Xte)
    f=pca.transform(Xfinal)
    m=LinearRegression().fit(Xtr_pc, np.ravel(ytr))
    r2_tr=m.score(Xtr_pc, ytr)
    r2_te=m.score(Xte_pc, yte)
    mse_tr=mean_squared_error(ytr, m.predict(Xtr_pc) )
    mse_te=mean_squared_error(yte, m.predict(Xte_pc) )
    R2[i]=[r2_tr, r2_te]
    mse[i]=[mse_tr,mse_te]
    shape=Xtr_pc.shape
    bic[i]=shape[0]*np.log(mse_tr)+np.log(shape[0])*shape[1]
    Pre[i]=m.predict(pca.transform(Xfinal))

R2_T=pd.DataFrame.from_dict(R2).T
R2_T.columns=["train_r2", "test_r2"]
mse_T=pd.DataFrame.from_dict(mse).T
mse_T.columns=["train_mse", "test_mse"]
bic_T=pd.DataFrame(bic.values(), index=bic.keys())
bic_T.columns=["bic"]
pd.concat([R2_T, mse_T, bic_T], axis=1).round(3)

	train_r2	test_r2	train_mse	test_mse	bic
3	0.968	0.970	203883.928	213033.191	8724.122
5	0.969	0.971	197169.014	206242.812	8713.414
10	0.986	0.988	87309.999	87207.485	8166.262
15	0.989	0.989	67615.908	76880.491	8017.099
20	0.991	0.990	60564.191	70367.966	7971.520
25	0.991	0.990	57756.695	68137.312	7970.566
30	0.991	0.991	57149.145	66350.432	7995.877
34	0.991	0.991	57133.233	66570.655	8021.951

위 결과에서 주성분의 수를 3개로 한다는 것은 변환한 설명변수가 3임을 의미합니다. 그 주성분 수가 20 이상이 되면서 위의 완전모델(full model)의 통계량(결정계수 등)과 유사한 결과를 보입니다. 이러한 결과는 각 모델의 MSE와 BIC를 나타내는 그림 8.5.2와 함께 고려할 때 더욱 명확해 집니다.

그림 1은 위 결과들을 시각화한 것으로 모델 복잡성을 나타내는 BIC의 경우 주성분 수가 20~ 25 사이에서 선택할 수 있습니다. 이 지점부터 MSE의 변화와 훈련과 검정 데이터의 차이 모두 안정됨을 보여줍니다.

fig, ax=plt.subplots(figsize=(5, 3))
line1, =ax.plot(mse_T.columns, mse_T.iloc[0,:], c="blue",label="mse_tr")
line2, =ax.plot(mse_T.columns, mse_T.iloc[1,:], c="green",label="mse_te")
ax.set_xlabel("n_components")
ax.set_ylabel('mse')
ax2=ax.twinx()
line3, =ax2.plot(bic.keys(), bic.values(), c="red", label="BIC")
ax2.set_ylabel("BIC")
plt.legend(handles=[line1, line2, line3], frameon=False)
plt.show()

다음 결과는 final 데이터에 대한 각 모델에서의 예측치 입니다.

pd.DataFrame(Pre).round(1)

	3	5	10	15	20	25	30	34
0	27039.0	27214.0	26615.0	26409.0	26445.0	26480.0	26551.0	26595.0

위 결과에 의하면 34개의 설명변수를 주성분에 의해 20 ~ 25개로 축소한 모델의 경우 모든 변수들을 사용한 모델과 유사한 결과를 보이며 모델 복잡성에서 개선된다고 할 수 있습니다.

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

유사변환(Similarity transformation) n×n 차원의 정방 행렬 A, B 그리고 가역 행렬 P 사이에 식 1의 관계가 성립하면 행렬 A와 B는 유사행렬(similarity matrix)이 되며 행렬 A를 가역행렬 P와 B로 분해하는 것을 유사 변환(similarity transformation) 이라고 합니다. $$\tag{1} A = PBP^{-1} \Leftrightarrow P^{-1}AP = B $$ 식 2는 식 1의 양변에 B의 고유값을 고려한 것입니다. \begin{align}\tag{식 2} B - \lambda I &= P^{-1}AP – \lambda P^{-1}P\\ &= P^{-1}(AP – \lambda P)\\ &= P^{-1}(A - \lambda I)P \end{align} 식 2의 행렬식은 식 3과 같이 정리됩니다. \begin{align} &\begin{aligned}\textsf{det}(B - \lambda I ) & = \textsf{det}(P^{-1}(AP – \lambda P))\\ &= \textsf{det}(P^{-1}) \textsf{det}((A – \lambda I)) \textsf{det}(P)\\ &= \textsf{det}(P^{-1}) \textsf{det}(P) \textsf{det}((A – \lambda I))\\ &= \textsf{det}(A – \lambda I)\end{aligned}\\ &\begin{aligned}\because \; \textsf{det}(P^{-1}) \textsf{det}(P) &= \textsf{det}(P^{-1}P)\\ &= \textsf{det}(I)\end{aligned}\end{align} 유사행렬의 특성 유사행렬인 두 정방행렬 A와 B는 'A ~ B' 와 같

sons dataStory

이 블로그 검색

[seaborn] seaborn의 함수에 사용하는 속성(인수)들

[data analysis] 주성분분석(principal component analysis)

주성분분석

관련된 내용

태그

댓글

댓글 쓰기

이 블로그의 인기 게시물

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

[matplotlib] 히스토그램(Histogram)

R 미분과 적분