주어진 데이터에 대해 가장 적합한 확률 분포 찾기 - 확률분포의 추정

이번에는 확률분포를 추정하는 판↘타↗스틱한 일을 해 보려고 해요. 계속 확률분포를 따졌으니, 당연히 이건 호기심이 생길 수밖에 없는 주제라고 생각합니다. 

확률분포를 추정한다는 것은 데이터만 가지고 짠!하고 확률분포를 찾아낸다는 의미는 아니고, 사실 특정한 확률분포를 가정하고, 관측값들을 이용해서 parameter(모수)를 찾아낸 후에 그 모수로 모델링한 분포와 실제 데이터가 얼마나 비슷한지 찾아내는 것이라고 생각하면 쉽습니다. 보통 parameter는 θ라고 표시하고요, μ, σ를 의미합니다. 왜 θ라고 하는지 잘은 모르겠습니다만, 참으로 묘한 통계입니다. 어쨌든 θ를 찾아내는 일인데, θ를 찾아낸 후 관측한 데이터와 가정한 확률분포가 얼마나 잘 정합되는가를 찾는 과정입니다. 당연히 특정한 확률분포를 어떻게 가정하는가? 하는 질문이 나오는 것이 당연합니다. 

그래서, 앞서 확률에서 여러 가지 경우의 확률 분포를 보긴 했습니다만, 그것만 보고는 과연 알아낼 수 있을까? 하는 조금 불안한 감이 없지 않아 있죠. 그런 접근법이 3가지 정도 있습니다. 

① 어떤 종류의 현상을 관찰한 것인지에 따라 우리가 이미 알고 있는 확률 분포를 가정할 수 있습니다. 이것은 이미 각각의 확률 분포가 어떤 현상을 설명하는지 알아봤으니까 조금 감이 있겠죠.

② 또는 데이터의 히스토그램을 그려보고 알고 있는 분포와 가장 비슷한 확률 분포를 이거 아닐까? 하고 가정할 수 있습니다. 확률 분포를 많이 알아둬야겠군요. 앞서 얘기한 확률분포의 족보가 조금은 도움이 될 거라 생각합니다. 

③ 마지막으로 이래도 저래도 어럽다 하면, 데이터의 타입을 보고 확률 분포를 가정하기도 합니다. 

예를 들어,

이런 식으로 가정하기 시작하면 꽤나 괜찮은 접근을 할 수 있습니다. 

그래도 좀 불안한 감이 없지 않아 있는데, 이 불안감을 불식시키기 위해 쓰는 실무적인 방법이 있습니다. 이제까지 말한 것과 크게 다르지는 않은데요,  

여러 가지 가능성 있는 분포를 가정하고, 그중 가장 적정한 분포를 Kolmogorov-Smirnov(콜모고로프-스미르노프) test를 이용하여 찾아내는 방법이 있습니다. 물론 SPSS나 MiniTab에서도 이런 기능이 있을 텐데, Python을 이용하면 조금은 있어 보이니까, 이 방법을 사용해 보도록 하겠습니다. 

감마분포 데이터를 따로 구하기 어려우니까, 감마분포의 데이터를 한번 만들어 보겠습니다. 

shape, scale = 2., 2.  # mean=4, std=2*sqrt(2)
data = np.random.gamma(shape, scale, 100) 

※ 생성되는 데이터가 그때 그때 달라서 아래의 fitting 결과와 다를 수 있습니다.

이렇게 하면 100개의 감마분포를 따르는 데이터를 만들어 낼 수 있습니다. 손쉬운 세상이군요. 

그런 다음, 분포를 한번 그려보겠습니다. 

음.. 감마랑 비슷해 보이나요. 뭐 그럭저럭 비슷한 것 같습니다. 이거 분포를 보니까, Gaussian, Gamma, Beta, Exponential, t, χ²가 비슷할 것 같은데요. 이 중에서 한번 찾아보겠습니다. 

아이디어는 다음과 같이 간단합니다. 
ⓐ 미리 준비된 후보군의 본포의 특징을 확인하고,
ⓑ Observation의 특징을 맞춰봅니다.
ⓒ 각 후보군과의 비교 (콜모고로프-스미르노프 검증) 후 p_value를 계산합니다.
ⓓ 그중 p_value가 가장 큰 값을 선택하는 방법입니다. 

Kolmogorov-Smirnov 검정의 원리는 별거 없습니다. 우리가 가정한 확률분포의 누적확률분포(cdf)와 실제 데이터의 최대 차이를 통계량으로 사용하여 그것을 보고 확률적으로 얼마나 그런 일이 벌어질 것인지를 p_value 산출하는 검정 방법입니다. 

- 갑자기 p_value를 얘기해서 미안합니다만, p_value에 대해서는 더 자세하게 얘기하겠지만, 여기에서는 p_value가 크다 함은 더 비슷하다는 의미 - 즉, 최대거리가 작다는 - 로 일단 받아들여주면 좋겠습니다. 여기에서 p_value를 설명하려면 원래 설명하려던 것을 설명하지 못할 지경에 이르니까, 조금만 참아주세요. -

그런데, Kolmogorov-Smirnov(콜모고로프-스미르노프) 검정을 구현하려면 너무 힘드니까, scipy에 준비되어 있는 패키지를 이용해서 가장 적합한 분포를 찾아내 보겠습니다. 

import scipy
import scipy.stats 

dist_names = ["norm", "gamma", "beta", "expon", "t", "chi2"]
dist_results = []
params = {}
for dist_name in dist_names:
    dist = getattr(scipy.stats, dist_name)  # 분포의 특징을 가져와서
    param = dist.fit(data) # 관측값을 이용하여 가정한 분포의 모수를 추정 
    params[dist_name] = param # 추정한 모수를 Dictionary에 넣고
    Stat, p = stats.kstest(data, dist_name, param) # Kolmogorov-Smirnov 검정을 합니다. 
    dist_results.append((dist_name, p)) # 결과를 list에 넣고요, 

print(dist_results) # 각 분포의 p value를 한 번 보고요. 
best_dist, best_p = (max(dist_results, key=lambda item: item[1])) # 그중에 제일 p값이 큰 것을 찾고요
print("제일 비슷한 분포: %s"%(best_dist))
print("그 분포의 p_value: %f"%(best_p))

결과 : 
[('norm', 0.018913807801915827), 
('gamma', 0.7659252892222164), 
('beta', 0.5599904607881772), 
('expon', 0.013608510683797844), 
('t', 0.19858126886916527), 
('chi2', 0.7658997548318363)]
제일 비슷한 분포: gamma
그 분포의 p_value: 0.765925

결과를 보니, p_value가 가장 큰 Gamma분포가 제일 비슷하군요 - 사실 표본수가 작아서 χ²와도 비슷하다고 생각할 수 있겠네요. 그리고 p_value가 가장 커서 비슷하다고는 했지만, 실제로는 Gamma분포이다라는 가설을 기각할 수 없다는 의미입니다요 - 다행히 Gamma분포라고 Kolmogorov-Smirnov검정이 결론 내주었군요. Gamma 분포가 여러 개의 후보 분포중 가장 데이터를 잘 설명한다고 생각하면 좋겠습니다. - 세상 좋아졌습니다. 사실 표본수가 작아서 모수 예측은 조금 날림입니다만, 표본수가 많아지면 모수도 잘 예측할 것이라 생각합니다. 

모수Parameter θ도 찾아내어 주었는데, 그러면 이 모수는 어떻게 추정했을까요? 이 모수를 추정해야 가장 비슷한 분포를 검정했을 테니까요. 이 θ를 찾아내는 방법이 여러 가지가 있습니다. 이게 유명한 MLE(Maximum Likelyhood Estimation - 최대우도법), Method Of Moment (모멘트 방법), 베이즈 추정법 등이 있습니다.  이 방법은 지금 말하면 정말 머리가 지끈거릴 테니까, 추정을 하면서, 한 번씩 해 보도록 하겠습니다. - 후후 -

마지막으로 모수를 추정한 후에 진짜 그 분포를 잘 따르는지 확인 하는 방법이 Kolmogorov-Smirnov 검정처럼 또 있는데 χ² 검정을 통해서 진짜 그 분포를 잘 따르는지 검정할 수 있습니다. 이것도 검정을 할 때 한번 들여다보면 좋을 거라 생각합니다. 

진지하게 뭔가를 해 본 건 아니지만, 그래도 확률 분포 추정은 이렇게 한다고 생각하고, 뭐 뿌듯하다고 생각하면 조금은 기분이 홀가분할 거라 생각합니다. 

사실 적합한 확률을 고르는 방법이라니, 아무도 발견하지 못한 확률분포가 또 있는 것은 아닐까 하는 의문이 들 때가 많이 있습니다. 당연히 아무도 발견 못한 확률분포가 있을 수 있습니다. 하지만, 정말 정말 이상한 경우가 아니고서는 소개한 확률 분포중에 하나를 따르게 되고, 그렇지 못한 경우에 데이터를 가지고 확률분포를 Estimation 하는 방법은 다음에 얘기할  Non Parametic Estimation을 이용하게 됩니다. 그렇다고 해서 매직키처럼 엄청난 결과가 나오는 것은 아니지만요.

보통 pdf를 표현할 떄, $ f(x;\theta)  $ 이런 식으로 θ가 문득 나올 때가 있는데, θ는 보통 고정된 값을 의미하는데 - 모수는 상수이니까요 -, 다음과 같이 표현이 많이 됩니다. $ f(x;\theta) = f(x; \mu , \sigma ^2) = \dfrac{1}{\sqrt{2\pi \sigma^2}} e^{-\dfrac{(x-\mu^2)}{2\sigma^2}}$ 이런 식입니다. 이럴 때 $\theta = (\mu, \sigma^2)$ 뭐, 결국 모수를 말하는 것입니다. 여기에서 μ가 location 모수, σ가 scale 모수로 볼 수 있는데, μ가 분포의 위치를 σ가 분포의 퍼진 정도를 표현해서 그렇습니다. 

이 방법 이외에도 확률분포 추정 방법은 꽤 많이 있습니다만, 지금 설명한 것보다 훨씬 더 전문적인 방법들이니까, 지금은 그냥 모르는 척 넘어가려고 합니다. 그런 방법들을 더 다루다가, 밤을 새워도 모자랄 테니, 그냥 이 정도에서 음음 그렇군 정도로 넘어가는 것이 정신 건강에 좋을 것 같아서요.