[SciPy] 2. Density Estimation (Application to Scatter Plot)

두번째로 작성하는 Scipy 주제입니다.
사실 제가 저번 시간 선형 회귀로 Scipy 주제를 시작한 계기가
바로 오늘의 Density Estimation 입니다.
_{따라서 3편은 (아직) 예정이 없다는 얘기 ㅋㅋ}

이번에도 먼저 아래 Scatter Plot 산점도로부터 시작해 보겠습니다.

이 그림에 대한 자세한 설명은 역시 생략합니다!
다만 저번과 마찬가지로 X축, Y축, 그리고 Color 이렇게 세가지 변수로
각각의 데이타, 즉 한 점이 정의되었음을 알려드립니다.
그리고 검은 Dash Line은 제가 비교를 위해 임의로 그려넣은 것이라 신경 안쓰셔도 됩니다.

Scatter Plot 산점도가 직관적인 모양새를 파악하는 데 뛰어난 장점이 있습니다.
그런데 문제는 2차원 평면 상의 어느 범위에 너무 많은 데이타가 몰릴 경우
점들이 겹쳐지면서 착시 효과가 나타날 수 있습니다.
그래서 자료의 Density 밀도를 주의해서 볼 필요가 있습니다 .

사실 자료의 밀도가 가장 중요한 분석 요소라면,
처음부터 Scatter Plot 산점도를 그릴 게 아니라 2-D Histogram을 그리면 되죠.
그런데 Scatter Plot 산점도와 히스토그램은 목적이 상이하고,
더 다양한 정보를 하나의 평면에 표시하기 위해
굳이 밀도를 계산하여 Contour를 위에 겹치도록 하겠습니다.

자료의 밀도 계산은
scipy --> stats --> kde에 있는 gaussian_kde를 이용할 예정입니다.
kde는 kernel-density estimation의 약자이며 wiki page에 기술된 바로는

In statistics, kernel density estimation (KDE) is a non-parametric way to estimate the probability density function (PDF) of a random variable.

즉, 무작위 변수의 확률밀도함수(PDF)를 근사하는 방법론 입니다. ~~한글이 더 어려운 느낌~~ 그리고 "gaussizn_kde"는 짐작하시듯 가우시안 정규분포를 가정하고 밀도를 근사하는 방식입니다.

간략화된 코드는 다음과 같습니다.

from scipy.stats import kde

### 인풋 데이타를 가지고 gaussian_kde 클래스를 정의합니다. (밀도 함수가 계산됩니다.)
##* xdata와 ydata 모두 여기서는 1-D 입니다. 
k = kde.gaussian_kde([xdata,ydata])

### 계산된 밀도 함수를 표현할  범위의 격자를 만듭니다. 
##* mgrid에 대해서는 부연설명1 참조
xi, yi = np.mgrid[xdata.min():xdata.max():300j,ydata.min():ydata.max():300j]

### 근사된 밀도함수 값을 위에서 정의된 격자 값에 따라 계산합니다.
##* 입력 자료의 형태에 대해서는 부연설명2 참조 
zi=k(np.vstack([xi.ravel(),yi.ravel()]))  

### Contour plot을 그립니다. 
### 옵션으로 5개의 Contour만 그리도록 지정되었습니다. 
ax.contour(xi,yi,zi.reshape(xi.shape),5,colors='k',linewidths=1.)

그래서 위의 코드를 추가한 결과 위의 Scatter Plot은 다음과 같이 바뀌었습니다.

비슷한 모양 (분포)인 줄 알았는데, 밀도는 확연히 다르군요!

부연설명 0
전과 마찬가지로 혹시 자료에 Missing Value 혹은 Undef 값이 있을 경우,
확실히 제거한 후에 밀도를 계산하는 게 좋습니다.
Numpy에서는 다음과 같이 "조건부 인덱싱"을 이용하면 쉽게 제거할 수 있습니다.

idx=np.logical_and(xdata>=0., ydata>=0)
xdatatmp=xdata[idx]
ydatatmp=ydata[idx]

부연설명 1
Numpy의 mgrid는 함수가 아닙니다! 그래서 끝에 "괄호()"가 아니라 Array의 일부를 뜻하는 "[ ]"가 쓰였습니다.
숫자 끝에 "j"가 붙으면 원래는 Complex number 복소수를 뜻하는데, mgrid의 경우에는 특별히 격자 개수로 해석됩니다. 위 예제에서 끝에 300j가 나오므로, 주어진 범위를 300등분 한 격자가 만들어집니다. 그래서 위의 xi와 yi는 각각 [300,300] 크기의 Array입니다.

부연설명 2
가우시안 근사 함수의 값을 주어진 격자 값으로 계산할 때 필요한 형태는 [# of dims, # of data] 입니다. 그래서 위에서는 Numpy의 vstack을 이용하여 x값과 y값을 붙여 하나의 Array로 만들었습니다. 예를 들어 길이 10인 2개의 1-D Array x를 vstack을 이용하여 붙이면 (e.g., np.vstack((x,x)) 입력은 튜플!) [2,10]인 Array가 만들어집니다.

부연설명 3
Numpy에서, ravel() flatten() reshape(-1) 어느 것을 써도 위의 밀도 계산 및 그리기 기능에는 전혀 차이가 없습니다. 다만, flatten()과 그 외 다른 2개의 차이점은, flatten()는 copy를 만들고, 다른 2개는 형태가 그렇게 보이도록 바꿔줄 뿐입니다.

[SciPy] 2. Density Estimation (Application to Scatter Plot)

Matplotlib List

F2PY List

Scipy List