기본 콘텐츠로 건너뛰기

[matplotlib] 등고선(Contour)

공분산의미(ing)

공분산이란

다음은 일정한 평균과 분산을 가진 정규분포에서 미지의 값을 신뢰할 수 있는 구간과 평균, 표준편차의 관계를 나타내면 다음과 같습니다.

그림 1. 표준 정규분포 상에서 평균과 표준편차의 관계

위 그림에서 빨간 점선, 파란 점선, 녹색 점선은 분포의 전확률 중에서 누적확률 68,3%, 95.4%, 99.7%를 차지하는 경계를 나타내며 각 구간의 value(값)은 다음과 같습니다.
red 점선 : CDF=68.3%, μ±σ
blue 점선 : CDF=95.4%, μ±2σ
green 점선 : CDF=99.7%, μ±3σ

위 그림의 python 코드는 다음과 같습니다.

plt.plot(x, rv.pdf(x), lw=2, label='pdf')
plt.xlabel("value")
plt.ylabel("Prob")
lineref=[mu, rv.ppf(0.683), rv.ppf(1-0.683), rv.ppf(0.954), rv.ppf(1-0.954), 
         rv.ppf(0.997), rv.ppf(1-0.997)]
cl=['black', 'red', 'red', 'blue', 'blue', 'green', 'green']

#plt.axvline(lineref[1], linestyle='--', color=cl[1])
for i in range(len(lineref)):
    plt.axvline(lineref[i], linestyle='--', color=cl[i])


plt.text(-1, 0.3, "mu+-sigma", color="red")
plt.text(-1, 0.2, "mu+-2*sigma",color="blue")
plt.text(-1, 0.05,"mu+-3*sigma", color="green")
plt.show()

위 데이터의 분산은 다음과 같이 계산됩니다.
σ2=i=1n(xiμ)2n1=E[(XE(X))(XE(X))]=σ(x,x)

그러나 분산은 단지 데이터의 퍼짐정도를 나타냅니다.
다음 그래프는 두 데이터 그룹의 경향을 보입니다.
그림 2. 데이터 셋 x, y의 산점도 
x=np.linspace(-3, 3, 100)
y=x+np.random.randn(100)
plt.scatter(x, y)
plt.xlim(-5, 5)
plt.ylim(np.min(y)-2,np.max(y)+2)
plt.show()

위 데이터들에서 x축에 따른 분산은 σ(x,x) 이며 y 축에 따른 분산은 σ(y,y) 로 나타낼 수 있습니다. 그러나 이 분산 정도는 각각 하나의 축에 대한 데이터의 퍼짐정도만을 나타내는 것으로 위 x와 y를 대응시킨 상태에서의 퍼짐 정도를 나타내지는 못합니다. 그림 2의 경우는 x의 증가에 따른 y의 증가가 명확해 보입니다. 즉, 그 둘의 관계는 양의 상관성을 보입니다. 
다음은 x, y의 상관계수 행렬을 나타낸 것입니다. 

np.corrcoef(x, y)
Out[68]: 
array([[ 1.        ,  0.87420057],
       [ 0.87420057,  1.        ]])

위의 상관계수는 다음과 같이 계산됩니다. 
r=(xx¯)(yy¯)(xx¯)2(yy¯)2=σ(x,y)σ(x,x)σ(y,y)
위 식에서 분모의 경우는 x, y의 각 분산이며 분자는 공분산이 됩니다.  이 관계에서 즉, 변수 두개에서 계산할 수 있는 분산의 경우의 수는 4가지가 됩니다. 
σ(x,x),σ(y,y),σ(x,y),σ(y,x)
x가 증가하면 y가 증가하는 양의 상관성을 가진다면 그 역의 관계 역시 성립되므로 다음의 관계가 성립되지요.
σ(x,y)=σ(y,x)
다음과 같이 공분산 행렬을 생성할 수 있으며 항상 대각원소들 을 중심으로 대칭이되는 대칭행렬이 됩니다. 
Σ=[σ(x,x)σ(x,y)σ(y,x)σ(y,y)]
이와 같이 그림 2와 같이 2차원의 공분산 행렬은  2×2 공분산 행렬이 생성됩니다. 유사하게 3차원의 입체 그림으로 나타나는 경우는 3×3 공분산 행렬이 생성되며 일반화하면 N 차원인 경우는 N×N의 공분산 행렬이 생성됩니다. 

그림 2의 데이터에서 각 데이터의 분산과 공분산을 계산해 보면 다음과 같습니다. 

np.var(x, ddof=1)
Out[80]: 3.0915212733394553

np.var(y, ddof=1)
Out[81]: 4.4333637534151134

np.sum((x-np.mean(x))*(y-np.mean(y)))/(len(x)-1)
Out[77]: 3.2364126231653638

위의 값은 np.cov()를 사용하여 공분산 행렬을 생성할 수있습니다. 
Out[78]: 
array([[ 3.09152127,  3.23641262],
       [ 3.23641262,  4.43336375]])

다음 그림3은 각 두 데이터의 상관성에 따른 공분산의 변화를 보여줍니다. 
그림 3. 대응되는 두 데이터 간의 상관성에 따른 공분산의 변화 (그림의 legend는 공분한행렬)

y1=np.flip(y, 0)
y3=np.random.randn(100)/100
x1=np.random.randn(100)/100

plt.subplot(2,2,1)
plt.scatter(x, y)
plt.xlim(-5, 5)
plt.ylim(np.min(y)-2,np.max(y)+2)
plt.text(0.5, -5, np.round(np.cov(x, y), 2))
plt.subplot(2,2,2)
plt.scatter(x, y1)
plt.xlim(-5, 5)
plt.ylim(np.min(y)-2,np.max(y)+2)
plt.text(0.5, 2, np.round(np.cov(x, y1), 2))
plt.subplot(2,2,3)
plt.scatter(x, y3)
plt.xlim(-5, 5)
plt.ylim(np.min(y)-2,np.max(y)+2)
plt.text(0.5, -5, np.round(np.cov(x, y3), 2))
plt.subplot(2,2,4)
plt.scatter(x1, y)
plt.xlim(-5, 5)
plt.ylim(np.min(y)-2,np.max(y)+2)
plt.text(0.5, -5, np.round(np.cov(x1, y), 2))

결과적으로 분산은 각 데이터 셋의 퍼짐(spread)을 나타내고 공분산은 두 데이터 그룹간의 퍼짐뿐만 아니라 방향성(orientation)을 함께 나타낸다고 할 수 있습니다. 

공분산의 고유분해(eigendecomposition)

공분산 행렬이 위의 데이터들을 변형 시킬 선형적 연산에 의해 어떻게 해석될 수 있는가를 봅시다. (고유값, 고유벡터, 고유분해 참조)
공분산은 데이터들 사이의 퍼짐뿐만 아니라 방향을 나타낸다는 것은 공분산 행렬로 부터 가장 큰 퍼짐의 방향을 지시하는 벡터를 발견할 수 있으며 그 방향에서의 퍼짐의 강도는 동일한 방향의 분산과 같을 것입니다.
이 벡터를 v로 정의한다면 이 벡터 위에 데이터들(D)을 투영시키면 vTD가 되며 그 투영된 데이터의 분산은 vTΣv이 됩니다.
이 계산에서 가장 큰 분산의 방향을 지시하는 v를 찾는 것이므로 투영된 데이터의 공분산 행렬 vTΣv 중에서 가능한 가장 큰 값을 갖는 원소를 선택해야 합니다. 표준화된 단위벡터인 v에 관해 vTΣv를 최대화 하는 것은 Rayleigh Quotient라 명명되는 것으로서 공식화할 수 있습니다. 이 최대화는 Σ의 가장 큰 고유벡터와 같은 v로 설정하는 것에 의해 얻어질 수 있습니다. 이 언급의 결과는 다음과 같습니다.
공분산 행렬의 가장 큰 고유벡터는 그 데이터의 가장 큰 방향을 지시하며 그 벡터의 크기는 대응되는 고유값과 같습니다. 두번째로 큰 고유벡터는 가장 큰 고유벡터에 항상 직교하고 그 데이터의 두번째로 큰 퍼짐의 방향을 지시합니다.

vTΣv으로 부터 v을 계산하기 위해 다음식을 유도할 수 있습니다.
Σv=λx
vΣ의 고유벡터이고 λ는 대응하는 고유값입니다.
그림 3의 각 데이터에 적용해 보면 각각의 고유값과 고유벡터는 다음과 같습니다.
np.linalg.eig(np.cov(x, y))
Out[246]:
(array([ 0.45716888,  6.63341568]), array([[-0.75727825, -0.65309238],
        [ 0.65309238, -0.75727825]]))
np.linalg.eig(np.cov(x, y1))
Out[247]:
(array([ 0.41454293,  7.00670062]), array([[-0.77065861,  0.63724823],
        [-0.63724823, -0.77065861]]))
np.linalg.eig(np.cov(x, y3))
Out[248]:
(array([  3.09152129e+00,   9.33124710e-05]),
 array([[  9.99999998e-01,   6.66444911e-05],
        [ -6.66444911e-05,   9.99999998e-01]]))
np.linalg.eig(np.cov(x1, y))
Out[249]:
(array([  7.46810051e-05,   4.32972227e+00]),
 array([[ -1.00000000e+00,  -2.46096962e-06],
        [  2.46096962e-06,  -1.00000000e+00]]))

그림 3의 왼쪽 위 그래프에서 큰 고유값에 해당하는 고유벡터의 방향을



그림 4. 공분산행렬의 고유분해 결과 

그림 4에서 붉은 선은 위 데이터에서 넓게 퍼진 방향을 나타냅니다. 이 방향은 큰 고유값에 대응하는 벡터의 방향과 같습니다. 또한 위 데이터가 빨간 선으로 투영된다면 그 투영된 데이터의 분산은 다음식으로 계산됩니다. 
vTΣv
np.dot(np.dot(ve1[:,1].T, np.cov(x, y)), ve1[:,1])
Out[267]: 6.6334156754030547
이 값은 빨간 선을 나타내는 고유벡터에 대응되는 고유값과 같습니다.  그러므로 각 벡터에 대응되는 고유값은 그 벡터의 방향으로 투영된다면 계산되는 분산과 같습니다. 
만약에 공분산 행렬에서 공분산이 0이라면 어떤 축과 동일하게 위치하므로 투영된 것이나 원 데이터들이나 동일하기 때문에 생성되는 고유값이 그대로 원 데이터의 분산이 되겠지요. 

선형변형으로서 공분산(Covariance matrix as a linear transformation)







댓글

이 블로그의 인기 게시물

[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 실수...