혼합모델 군집화(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 으로 그려보겠습니다.


먼저 "mvtnorm" 패키지의 dmvnorm() 함수로 다변량 정규분포 함수를 정의해서 3차원 그래프를 그리는 방법은 아래와 같습니다. 



# Multi-variable Normal Distribution

install.packages("mvtnorm")

library(mvtnorm)


x <- seq(from = -5, to = 5, by = 0.2)

y <- x


gaussian_func <- function(x, y) {

  dmvnorm(cbind(x, y))

}


persp(x, y, outer(x, y, gaussian_func), theta = 30, phi = 10)


 




이번에는 위 포스팅에서 소개한 다변량 정규분포 공식을 이용해서 직접 gaussian_func() 를 정의해보고, 이 함수를 사용해서 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)'에 대해서 알아보겠습니다.

 

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

 

728x90
반응형
Posted by Rfriend
,

다수의 변수간 상관관계를 파악하려고 할 때, 회귀분석에서 종속변수와 독립변수간 선형관계를 파악하거나 독립변수간 다중공선성을 파악하려고 할 때 사용하는 분석 기법이 상관계수 행렬이며, 시각화 방법이 산점도 행렬과 상관계수 행렬 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

 

-----------

 

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

 

728x90
반응형
Posted by Rfriend
,