혼합모델 군집화(mixture model clustering)의 혼합 모델에는 어떠한 분포도 가능하지만 일반적으로 정규분포(Gaussian normal distribution)을 가장 많이 사용하며, 다음번 포스팅에서는 가우시안 혼합 모델(Gaussian misture model, 줄여서 GMM) 군집화를 다루겠습니다. .

 

혼합모델 군집화(mixture model clustering)에 들어가기에 앞서, 이번 포스팅에서는 정규분포(Gaussian normal distribution)에 대해서 먼저 몸풀기 차원에서 짚어보고 넘아가겠습니다.  준비운동 안하고 혼합모델 군집화에 바로 투입됐다가는 자칫 머리에 쥐가 나거나 OTL 하게 되는 불상사가 예상되기 때문입니다. ^^; 

 

이번 포스팅은 쉬어 가는 코너, 한번에 다 몰라도 좌절하지 않기 코너, 생각날 때 찾아보고 복기하는 코너로 여기면 좋겠습니다.

 

기계학습, 다변량 통계분석 책을 보다 보면 수식, 표기법 때문에 많이 어려워했던 기억이 있을 겁니다.  가장 기본이 되고 많이 사용하는 개념, 수식, 표기법을 중심으로 소개해보겠습니다.

 

 

1) 모집단(Population) vs. 표본집단(Sample)

 

알고자 하는 '전체 객체, 성분, 관측치를 포함하는 집단'을 모집단이라고 하며, 보통은 시간, 비용 상의 이슈로 인해서 모수(parameter)가 알려지지 않은(unknown) 경우가 많습니다.

 

모집단의 특성을 알기 위해 시간과 비용을 감당할 수 있는 범위 내에서 모집단으로 부터 모집단을 잘 대표할 수 있도록 표본(sample)을 추출해서 분석을 하게 됩니다.

 

 

모집단과 표본집단을 알았으니, 각 집단의 기대값과 공분산, 상관행렬, 상관계수는 어떻게 표기하는지 비교해서 살펴보겠습니다.

 

 

2) 모집단(Population)의 기대값과 공분산, 상관행렬, 상관계수

 

성분들이 확률변수로 이루어진 p x 1 확률벡터 X에 대한 평균 벡터(mean vector), 공분산 행렬(covariance matrix), 상관행렬(correlation matrix), 상관계수(correlation coefficient) 는 아래와 같습니다. 

 

공분산 행렬(covariance matrix)을 'Σ' 로 표기한다는 점 유념하시기 바랍니다.  다변량 정규분포 확률분포를 표기할 때 'Σ'가 나옵니다.  (* 주의 : 모두 summation 하라는 뜻이 아님 -.-").

 

 

 

위의 모집단(population)에 대한 표기법을 아래의 표본집단(sample)에 대한 표기법과 비교를 해가면서 유심히 살펴보시고 기억해놓으면 좋겠습니다. 

(강의 듣다보면 혼동해서 잘못쓰는 강사나 교재도 가끔 있음 -_-;) 

 

 

 

(3) 표본집단의 기대값과 공분산, 상관행렬, 상관계수

 

모집단으로부터 n개의 확률표본 X1, X2, ..., Xn에 대해 각 n x p 확률벡터 Xi = (Xi1, Xi2, ..., Xip)' , i = 1, 2, ..., n 의 표본 평균 벡터(sample mean vector), 표본 공분산 행렬(sample covariance matrix), 표본 상관행렬(sample correlation matrix), 표본 상관계수(sample correlation coefficient)는 아래와 같습니다.

 

표본 상관행렬의 우상(rigth upper)과 좌하(left down)의 값은 상호 대칭으로 같습니다. (가령, r12r21은 같은 값을 가짐)

 

 

 


 

 

정규분포는 변수의 개수(1개, 2개, p개)에 따라서 '일변량 정규분포 (univariate normal distribution)', '이변량 정규분포 (bivariate normal distribution)', 다변량 정규분포 (multivariate normal distribution)'로 구분할 수 있습니다.

 

일변량/이변량/다변량 정규분포 확률밀도함수(probability density function)의 표기법과 수식, 그리고 그래프를 살펴보도록 하겠습니다.

 

 

(4) 일변량 정규분포 (univariate normal distribution) vs. 다변량 정규분포 (multivariate normal distribution)

 

일변량 정규분포의 확률변수 X에 대한 확률밀도함수(probability density function) 식은 기억하기가 좀 복잡하고 어려운데요, 자꾸 보다보면 어느샌가 익숙해지고, 더 자주 보다보면 자기도 모르게 기억될 날이 올겁니다. (기억을 했다가도... 인간은 위대한 망각의 동물인지라 두어달만 안보면 다시... ㅋㅋ)

 

 

위의 표기에서 일변량 정규분포 대비 다변량 정규분포의 확률밀도함수의 어디가 서로 다른지 유심히 살펴보시기 바랍니다.  특히 다변량 정규분포의 확률밀도함수 표기법이 통계를 전공하지 않았거나, 선형대수 기본을 모르면 매우 생소하고 어렵게만 느껴질 것 같습니다. (왜 저런거냐고 묻는 댓글 질문은 정중히 사양해염... -_-*)

 

 

 

(5) 이변량 정규분포 (bivariate normal distribution)

 

3차원 이상 넘어가면 사람 머리로 인식하기가 곤란하니, 시각화로 살펴보기에 그나마 가능한 2개의 변수를 가지는 2차원의 '이변량 정규분포(bivariate normal distribution)'의 확률밀도함수(probability density function)에 대해서 알아보겠습니다.  위의 다변량 정규분포의 확률밀도함수 식에서 p = 2 를 대입하면 됩니다.

 

[ 이변량 정규분포 확률밀도함수 - 식 1 ]

 

 

 

위의 이변량 정규분포의 확률밀도함수를 다시 정리해보면 아래와 같습니다.

 

[ 이변량 정규분포 확률밀도함수 - 식 2 ]

 

 

 

'이변량 정규분포 확률밀도함수   식 2'를 보면 이건 뭐... 외계어가 따로 없습니다.  도대체 이렇게 복잡하면서도 쓰임새가 어마무시 많은 요물단지 정규분포를 정리한 천재들이 누구인지 궁금하시지요?

 

위키피디아를 찾아보니, 위대한 수학자 가우스(carl Friedrich Gauss)가 정규분포를 발견했고(그래서 영어로 Gaussian Normal Distribution 이라고 함), 라플라스(Marquis de Laplace)가 그 유명하고 중요하며 시험문제로 반드시 나오는 '중심극한의 정리(Central Limit Theorem)'를 증명하였습니다. 피어슨(Karl Pearson)이 20세기 들어서 '정규분포(Normal Distribution)'이라고 명명하는 것을 확산시키는데 일조했다고 나오네요.  가우스, 라플라스, 피어슨... 역시 천재 거물들 맞군요.

(Normal distribution 이외의 분포는 그럼 비정상(Abnormal) 이라는 소리? -_-?)

 

 

 

 

(* source : https://en.wikipedia.org/wiki/Normal_distribution)

 

 

 

수식으로만 봐서는 저게 뭔가 싶으실텐데요, R의 base 패키지에 내장되어 있는 persp() 함수를 가지고 위의 '이변량 정규분포 - 식 2'를 사용하여 3D plot 으로 그려보겠습니다.

 

persp() 의 3D plot에 manipulate 패키지를 접목해서 동적(interactive plotting)으로 상하(phi), 좌우(theta) 로 돌려가면서 살펴보겠습니다. 

 

아래 첫번째 그래프의 노란색으로 표시해 놓은 단추를 누르면 interactive plot 의 slider 가 나타납니다.  theta(좌~우), phi(상~하)를 조절해가면서 3-D 그래프를 돌려서 보시면 '아, 이변량 정규분포 확률밀도함수가 이런거구나. 느낌 오네~' 하실겁니다.

 

 

##---------------------------------------
## (Gaussian) Mixture Model Clustering
##---------------------------------------

 

# bivariate nomral probability density function
mu1 <- 0   # setting the expected value of x1
mu2 <- 0   # setting the expected value of x2

 

s11 <- 3   # setting the variance of x1
s12 <- 4  # setting the covariance between x1 and x2
s22 <- 3   # setting the variance of x2

 

rho12 <- s12/(sqrt(s11)*sqrt(s22)) # setting the correlation coefficient between x1 and x2
rho12

 

x1 <- seq(-5, 5, length = 50) # generating the vector series x1
x2 <- seq(-5, 5, length = 50) # generating the vector series x2


# multivariate normal density function
gaussian_func <- function(x1, x2){
  term1 <- 1/(2*pi*sqrt(s11*s22*(1-rho12^2)))
  term2 <- -1/(2*(1-rho12^2))
  term3 <- (x1 - mu1)^2/s11
  term4 <- (x2 - mu2)^2/s22
  term5 <- 2*rho12*((x1 - mu1)*(x2 - mu2))/(sqrt(s11)*sqrt(s22))
  term1*exp(term2*(term3 + term4 - term5))
}


# calculating the density values
z_score <- outer(x1, x2, gaussian_func)
head(z_score)

 

 

##------------
# 3-D normal distribution density plot : persp() of {base} package
# with interactive plotting with Manipulate package in RStudio

 

install.packages("manipulate")
library(manipulate)

 

manipulate(persp(x1, x2, z_score,
                 theta = theta_x, # theta gives the azimuthal viewing direction
                 phi = phi_x, # phi gives the colatitude viewing direction
                 main = "Two dimensional Normal Distribution"), 
           theta_x = slider(10, 90, initial = 35), 
           phi_x = slider(10, 90, initial = 10))

 


 

# click the button at yellow circle position

 

[ theta : 35,   phi : 10 ]

 


 

[ theta : 60,   phi : 10 ]

 


[ theta : 80,   phi : 10 ]


[ theta : 90,   phi : 10 ]


 

[ theta : 35,   phi : 30 ]


[ theta : 35,   phi: 50 ]


 

[ theta : 35,   phi: 70 ]


 

[ theta : 35,   phi: 90 ]

 

 

 

 

마지막 그림이 phi = 90 으로서, 꼭대기 상단에서 비행기 타고 내려다본 그림인데요, 이걸 2-D 의 등고선 그림(contour plot)으로 그려보면 아래와 같습니다.

 

 

# 2-D contour plot
contour(x1, x2, z_score, xlab = "x1", ylab = "x2",
        main = "Bivariate(2 dimensional) Normal Distribution : contour plot")

 

 

 

 

 

참고로, 일변량 정규분포 그래프는 http://rfriend.tistory.com/102 의 이전 포스팅을 참고하세요.

 

[Reference]

1. 김재희, 'R 다변량 통계분석', 교우사, 2011

2. Wikipedia :  https://en.wikipedia.org/wiki/Normal_distribution

 

 

이번 포스팅에서 살펴본 '다변량 정규분포'에 대한 기본기를 바탕으로 해서, 다음번 포스팅에서는 '가우시언 혼합 분포 군집화 (Gaussian Mixture Model Clustering)'에 대해서 알아보겠습니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요.

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 심교훈 2019.04.01 11:10 신고  댓글주소  수정/삭제  댓글쓰기

    정말 유용한 글 감사합니다!

다수의 변수간 상관관계를 파악하려고 할 때, 회귀분석에서 종속변수와 독립변수간 선형관계를 파악하거나 독립변수간 다중공선성을 파악하려고 할 때 사용하는 분석 기법이 상관계수 행렬이며, 시각화 방법이 산점도 행렬과 상관계수 행렬 Plot (correlation matrix plot) 입니다.

 

이전 포스팅에서 ggplot2의 geom_point() 산점도를 다루었으며,

 

다음 포스팅에서는 Base Graphics 패키지의 pairs() 함수를 사용한 산점도 행렬을 소개하였고,

 

이번 포스팅에서는 상관계수 행렬 Plot을 중심으로 해서 corrplot 패키지 사용법을 알아보겠습니다. 

 

 

예제로 사용한 데이터는 뉴욕의 1973년도 공기의 질을 측정한 airquality 데이터셋의 Ozone, Solar.R, Wind, Temp 4개의 변수가 되겠습니다.

 

 

> str(airquality)
'data.frame':	153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
> 

> # Month, Day는 빼기

> airquality_1 <- airquality[,c(1:4)]
> 
> str(airquality_1)
'data.frame':	153 obs. of  4 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...

 

 

 

 

상관계수 분석을 할 때 결측값이 있으면 NA 값이 나오게 되므로 사전에 결측값 처리하는 것이 필요합니다.  Ozone과 Solar.R이 결측값이 각각 37개, 7개 있다보니 아래처럼 상관계수가 NA가 나왔습니다.

 

> # 결측값 확인
> sum(is.na(airquality_1$Ozone)) # 37
[1] 37
> sum(is.na(airquality_1$Solar.R)) # 7
[1] 7
> sum(is.na(airquality_1$Wind)) # 0
[1] 0
> sum(is.na(airquality_1$Temp)) # 0
[1] 0
> # 결측값 있는 상태에서 상관계수 계산했을 때
> cor(airquality_1)
        Ozone Solar.R       Wind       Temp
Ozone       1      NA         NA         NA
Solar.R    NA       1         NA         NA
Wind       NA      NA  1.0000000 -0.4579879
Temp       NA      NA -0.4579879  1.0000000

 

 

 

 

na.omit() 함수를 사용하여 결측값이 있는 행 전체를 삭제한 후에 상관계수를 구해보면 아래와 같습니다.  corrplot 패키지의 corrplot() 함수는 상관계수 행렬 데이터셋을 가지고 그래프를 그리므로 아래처럼 결측값을 제거한 후의 데이터셋을 가지고 미리 상관계수 행렬을 계산해두어야 합니다.

 

 

> # 결측값 있는 행 전체 삭제
> airquality_2 <- na.omit(airquality_1)
> str(airquality_2)
'data.frame':	111 obs. of  4 variables:
 $ Ozone  : int  41 36 12 18 23 19 8 16 11 14 ...
 $ Solar.R: int  190 118 149 313 299 99 19 256 290 274 ...
 $ Wind   : num  7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...
 $ Temp   : int  67 72 74 62 65 59 61 69 66 68 ...
 - attr(*, "na.action")=Class 'omit'  Named int [1:42] 5 6 10 11 25 26 27 32 33 34 ...
  .. ..- attr(*, "names")= chr [1:42] "5" "6" "10" "11" ...
> sum(is.na(airquality_2$Ozone)) # 0
[1] 0
> sum(is.na(airquality_2$Solar.R)) # 0
[1] 0
> # 상관계수 계산
> airquality_cor <- cor(airquality_2)
> airquality_cor
             Ozone    Solar.R       Wind       Temp
Ozone    1.0000000  0.3483417 -0.6124966  0.6985414
Solar.R  0.3483417  1.0000000 -0.1271835  0.2940876
Wind    -0.6124966 -0.1271835  1.0000000 -0.4971897
Temp     0.6985414  0.2940876 -0.4971897  1.0000000
 

 

 

 

corrplot 패키지는 별도의 설치 및 호출이 필요한 패키지이므로 아래의 절차를 거칩니다.

 

> install.packages("corrplot")
Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/corrplot_0.73.zip'
Content type 'application/zip' length 2680505 bytes (2.6 MB)
downloaded 2.6 MB

package ‘corrplot’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\Rtmpk1gkRL\downloaded_packages
> library(corrplot) 

 

 

 

산점도 행렬 그림 (scatter matrix plot)을 복습해보자면 아래와 같습니다.

 

> # scatter plot matrix
> plot(airquality_2)

 

 

 

 

 

correlation plot의 method 에는 method = c("circle", "square", "ellipse", "number", "shade", "color", "pie") 등이 있으며, method별로 하나씩 예를 들어보겠습니다.

 

 

> corrplot(airquality_cor, method="circle")
 

 

 

> corrplot(airquality_cor, method="square")

 

 

 

  

> corrplot(airquality_cor, method="ellipse")

 

 

 

 

> corrplot(airquality_cor, method="number")

 

 

 

 

> corrplot(airquality_cor, method="shade")

 

 

 

 

> corrplot(airquality_cor, method="color")

 

 

 

 

> corrplot(airquality_cor, method="pie")

 

 

 

 

마지막으로 mehtod="shade", 상관관계 방향성 제시, 대각선 값 미제시, 상관계수 숫지 검정색으로 해서 추가해서 corrplot을 그려보겠습니다.  order 는 FPC(First Principle Component), hclust(hierarchical clustering), AOE(Angular Order of Engenvectors) 등이 있으며, 정렬 기준을 지정해주면 같은 색깔 끼리 뭉쳐서 보일 수 있도록 정렬을 시켜줘서 보기에, 해석하기에 더 좋게 보여줍니다.

 

> # corrplot
> corrplot(airquality_cor, 
+          method="shade", # 색 입힌 사각형
+          addshade="all", # 상관관계 방향선 제시
+          # shade.col=NA, # 상관관계 방향선 미제시
+          tl.col="red", # 라벨 색 지정
+          tl.srt=30, # 위쪽 라벨 회전 각도
+          diag=FALSE, # 대각선 값 미제시
+          addCoef.col="black", # 상관계수 숫자 색
+          order="FPC" # "FPC": First Principle Component
+                      # "hclust" : hierarchical clustering
+                      # "AOE" : Angular Order of Eigenvectors
+          )

 

 

 

 

많은 도움이 되었기를 바랍니다.

 

-----------

 

(참고) pairs() 함수를 활용한 상관계수 행렬 그리기 ☞ http://rfriend.tistory.com/83

 

-----------

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡' 단추를 꾸욱 눌러주세요.^^

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 스윗쥬스 2016.08.04 17:06 신고  댓글주소  수정/삭제  댓글쓰기

    order Option은 3가지가 있습니다.
    직접 테스트를 했는데도, 어떤 기준으로 정렬하는지 파악이 잘 안되네요.
    각 기준이 어떻게 정렬하는지 알려주실 수 있으실까요..?

    1.FPC(First Principle Component)
    2.hclust(hierarchical clustering)
    3.AOE(Angular Order of Engenvectors)

    • R Friend R_Friend 2016.08.04 23:00 신고  댓글주소  수정/삭제

      스윗쥬스님,

      order option별로 설명을 드리자면,

      1)FPC(First Principle Component)

      주성분분석이 데이터의 분산을 최대로 많이 설명할 수 있는 선형조합을 찾는 것이며, 제1 주성분이 데이터의 정보량을 가장 많이 설명을 합니다. 보통 변수들 간에 상관성이 높을 수록 제1 주성분의 설명력가 높습니다. 이론상으로는 제1 주성분에 대한 변수별 기여도를 가지고 ordering을 하면 서로 상관성이 높은 변수들의 순서대로 정렬이 될겁니다.


      ☞ 주성분분석(PCA) 참고 포스팅 :
      http://rfriend.tistory.com/61


      2)hclust(Hierarchical Clustering) :

      계층적 군집화 방법에는 6가지 세부 알고리즘이 있는데요, 기본 원리는 객체 간 유사성(similarity)을 측정해서 유사한 객체끼리 군집이 하나가 될때까지 반복해서 묶어주는 작업을 합니다. 보통은 Euclidean Distance 를 가지고 유사성(정확히는 비유사성, dissimilarity)를 측정해서 군집화를 해줍니다. 앞서 말씀드린 주성분분석 (혹은 요인분석)은 '변수'에 대해서 '상관성'을 가지고 주성분으로 묶어주는 것이라면, 군집분석은 '객체(관측치)'에 대해서 '유사성'을 가지고 군집으로 묶어주는 것입니다.

      따라서 상관관계분석의 ordering 기준으로 계층적군집분석을 쓴다는게 어떤 매카니즘인지를 저는 정확히 이해를 못하겠습니다. ^^; 죄송합니다.

      보통 군집분석 끝내고 나면 군집별로 특성 파악을 위해 변수들과 cross-tabulation 해가면서 profiling을 하는데요, 각 군집별로 특성을 잘 나타낼 수 있는 변수들이 있다면 그 변수들의 상관성이 높다고 해석할 수도 있을 것 같습니다.

      ☞ 응집형 계층적 군집화 - 단일(최단) 연결법 참고 포스팅 :
      http://rfriend.tistory.com/202


      3) AOE(Angular Order of Eigenvectors) :

      1)번에서 주성분분석을 얘기했었는데요, 주성분 구할 때 성분을 분해하면서 고유값(Eigenvalue)과 고유벡터(Eigenvector)를 사용합니다. 첫번째, 두번째, ... 이 순서대로 고유값이 컸다가 작아지구요, 설명력도 첫번째 고유값, 고유벡터가 제일 크고 점점 작아진다고 생각하시면 됩니다. 1)번의 FPC와 대략 같은 맥락의 order 기준이라고 생각하시면 될거 같습니다.

      ☞ 고유값(eigenvalue), 고유벡터(eigenvector) 참고 포스팅 : http://rfriend.tistory.com/181

    • 스윗쥬스 2016.08.05 16:15 신고  댓글주소  수정/삭제

      답변감사합니다 ~
      덕분에 많이 공부했습니다. ^^

  2. 지나가던학생 2016.11.09 20:26  댓글주소  수정/삭제  댓글쓰기

    저 질문이 있는데... 상관계수 구해서 시각화할때 유의성, 즉 p값은 고려하지 않고, 단순 값만 시각화 된것인가요? 만약 저 매트릭스에서 p값이 유의하지 않게 나온 값은 표시되지 않게 할수 있는 방법은 없을까요

    • R Friend R_Friend 2016.11.10 00:26 신고  댓글주소  수정/삭제

      pairs() 함수와 사용자 정의 함수를 사용하면 (두 변수들 간의) 산점도 행렬 + (한 변수의) 히스토그램 + (두 변수들 간의) 상관계수를 한꺼번에 볼 수 있습니다.

      특히 우측 상단의 상관계수의 절대값 크기에 비례해서 숫자 크기가 달라짐에 따라 상관계수절대값이 큰 값만 눈에 잘 띄며, 별로 관계가 없는 값은 숫자 크기가 매우 작게 나타나는 시각화 효과를 제공합니다.

      아래 링크의 포스팅 참고하세요.

      => http://rfriend.tistory.com/83