주성분분석
관련된 내용
공분산행렬의 고유벡터와 그에 대응하는 고유값을 고려하여 변수행렬의 차원을 감소시키면서 지정된 주축(고유벡터)을 기반으로 존재하는 축소된 차원의 변환 데이터를 생성할 수 있습니다. 이러한 주축 즉, 고유벡터들을 주성분(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=Xp1. p1는 가장 큰 고유값에 대응하는 고유벡터입니다.
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() 클래스를 사용한 선형모델에서 훈련과 검정 데이터에서의 R2, 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개로 축소한 모델의 경우 모든 변수들을 사용한 모델과 유사한 결과를 보이며 모델 복잡성에서 개선된다고 할 수 있습니다.
댓글
댓글 쓰기