2. 최대 가능도 추정

이번 장에서는 실제 데이터에서 얻어낸 분포를 정규분포에 적합<fit> 시키는 방법인 최대 가능도 추정<Maximum Likelihood Estimation, MLE>에 대해 다룬다. 

2.1 생성모델 개요

이제야 생성모델에 대해 간단하게나마 이야기를 꺼낸다.

2.1.1 생성모델 목표

생성모델이란 특정 데이터x의 확률 분포 p(x)를 모델링(Metric으로 표현)한 다음 그 집단에서 선택된 것 같은 유사 데이터를 새롭게 생성하는 것이다. 좀 쉽게 표현하면 실제 데이터와 매우 유사한 새로운 데이터를 생성해내는 것이다. 여기서 중요한 것은 실제 데이터와 매우 유사하려면 실제 데이터 집단의 특징을 잘 모델링 해야 한다는 것이다. 이것을 위해 이전 장에서 확률분포의 매개변수를 공부하고 계산했던 것이다. 집단의 특징을 모델링하기 위해서는 매개변수를 적절히 설정해야 하고 이를 추론<inference>라고 한다. 

2.1.2 모집단과 샘플

모집단은 대상의 전체 집합을 의미하는데 통계에서 문제는 모집단을 대상으로 측정이 사실상 불가능하기 때문에 우리는 샘플을 사용하여 측정을 한 뒤 이것을 통해 모집단을 추정한다. 여기서 생성모델의 맥락에서 보면 모집단은 샘플을 뒷받침하는 확률분포가 된다. 따라서 모집단에서 샘플을 추출하여 특징을 모델링하고 이것이 모집단을 잘 대변할 수 있도록 우리는 매개변수를 잘 추론해야 한다. 간단히 말하면 다음과 같다.

 

1. 모델링: 모집단 분포를 매개변수로 조정 가능한 확률분포로 비슷하게 표현할 수 있다고 가정한다. 다시말하면 적절한 확률 분포를 정의하는 과정이다. 

2. 매개변수 추론:: 모델링 결과로 만들어진 확률분포가 샘플 데이터에 부합하도록 매개변수를 추론한다. 

 

1번에서는 매개변수를 통해 제어할 수 있는 확률분포를 설정하고, 2단계에서 샘플과 가장 잘 부합하도록 매개변수를 조정한다. 이 2단계에서 우리는 최대 가능도 추정 기법을 사용한다.

2.2 실제 데이터로 생성 모델 구현

간단한 생성 모델을 구현해보자. 데이터셋은 SOCR 데이터로 1993년도 당시 홍콩인들 18세의 가상의 키 데이터 25000개이다. 

import numpy as np
import matplotlib.pyplot as plt
import os

path =  os.path.join(os.getcwd(), "height.txt")
xs = np.loadtxt(path)
print(xs.shape)

plt.hist(xs, bins='auto', density=True)
plt.xlabel('Height')
plt.ylabel('Probability Density')
plt.show()

다음과 같이 정규분포의 모양을 잘 따르는 분포도인 것을 알 수 있다. 이제 키의 확률분포를 정규분포로 모델링해보자.

2.2.2 정규분포를 따르는 생성모델

정규분포를 따르는 생성 모델을 구현하는데 여기서 해야 할 일은 위에서와 마찬가지로 2가지다.

1. 모델링: 키 데이터가 정규분포라고 가정((적절한 확률분포인데 여기서는 정규분포가 가장 적절하다.))

2. 매개변수 추론: 샘플을 기반으로 정규분포의 매개변수를 추론

해당 작업으로 생성모델을 만든다. 모델링은 이미 위 코드를 통해 정규분포를 따른가는 것을 직접 확인하였다. 따라서 2번인 매개변수를 추론하는 것만 하면 된다. 

 

여기서 사용되는 기법이 최대 가능도 추정이다. 

mu = np.mean(xs)
sigma = np.std(xs)
print(mu, sigma)

위 코드만 추가하면 위 데이터의 매개변수인 기댓값과 표준편차 값을 구할 수 있다. 간단히 확인해보면 다음과 같다.

import numpy as np
import matplotlib.pyplot as plt
import os

path =  os.path.join(os.getcwd(), "height.txt")
xs = np.loadtxt(path)
print(xs.shape)

def normal(x, mu=0, sigma=1):
    y = 1 / (np.sqrt(2*np.pi)*sigma) * np.exp(-(x - mu)**2 / (2*sigma**2))
    return y

x = np.linspace(150, 190, 1000)
mu = np.mean(xs)
sigma = np.std(xs)
y = normal(x, mu, sigma)

plt.hist(xs, bins='auto', density=True)
plt.plot(x, y)
plt.xlabel('Height')
plt.ylabel('Probability Density')
plt.show()

 

다음과 같이 코드를 실행해보면 키 데이터와 정규분포가 잘 맞아 떨어지는 것을 알 수 있다. 

2.3 최대 가능도 추정 이론

앞에서는 데이터를 통해 코드로 계산을 하여 정규분포의 매개변수들을 추론하였다. 여기서 사용된 식은 다음과 같다.

 

$\hat{\mu} = \frac{1}{N}\sum_{n=1}^{N}x^{(n)}$

 

$\hat{\sigma} = \sqrt{\frac{1}{N}\sum_{n=1}^{N}(x^{(n)} - \hat{\mu})^{2}}$

 

위 식의 성립함은 최대 가능도 추정 이론을 통해 확인이 가능하다. 

2.3.1 가능도 최대화

매개변수 $\theta$에 의해 모양이 결정되는 확률분포가 있다고 가정했을 때, 데이터 x를 얻을 수 있는 확률 밀도는 다음과 같다.

 

$p(x; \theta)$

 

다음 샘플(x1, x2, x3,  ....)을 얻은 경우를 가정하자. 각 데이터는 확률 분포에 따라 독립적으로 생성된다고 가정했을 때 이 때 샘플 D를 얻을 수 있는 확률 분포는 각각의 확률밀도를 전부 곱한 값이다.

 

$p(D; \theta) = \prod_{n=1}^{N}p(x^{(n)}; \theta)$

 

각각이 독립시행이기 때문에 단순 곱으로 계산이 가능하다. 이러한 위 식을 쎄타를 인수로 받는 함수 그 자체로 간주하여 다음과 같은 식도 정의할 수 있다.

 

$L(\theta) = p(D; \theta)$

 

여기서 L(()) 함수를 가능도 또는 가능도 함수 라고 한다. 가능도란 특정 매개변수(위 식의 경우 쎄타)에 대해 샘플 D가 일어날 가능성을 확률밀도로 표현해준다. 보통은 그냥 L()을 따로 정의하지 않고 p()로만 표현하기도 한다. 

 

최대 가능도 추정은 위에서 배운 가능도를 최대화 하는 매개변수 쎄타 값을 찾는 기법이다. 이것을 사용하는 이유는 매개변수가 이 특정 쎄타값일 때 샘플이 관찰될 확률이 가장 높기 때문이다. 다시 말하면 해당 매개변수를 가지는 모델이 샘플에 가장 잘 맞다는 의미이다. 실제로는 일반적인 가능도 p()를 계산하지 않고 로그 가능도인 log p(())를 사용하는 경우가 많다. 계산이 용이하기 때문이다. 그리고 log함수는 단조 증가 함수이기 때문에 최대 가능도 쎄타 값이 그냥 p(())에서의 최대 가능도 값과 동일하다. <여기서 log 함수의 밑은 자연상수 e를 사용한다.>

2.3.2 미분을 이용한 최댓값 찾기

최대 가능도를 찾아야 하기 때문에 함수의 최댓값을 찾아야 하는데 고등학교를 나왔으면 다 아는 미분을 사용하여 최댓값을 찾을 것이다. 대충 다 아는 내용이니 걍 넘어가겠다. 아무튼 수식을 풀어 해를 구하는 것을 해석적으로 풀 수 있다 또는 닫힌 해를 구할 수 있다 라고 표현한다. 

2.3.3 정규분포의 최대 가능도 추정

정규분포의 확률 밀도는 이전 1장에서 배웠다시피 다음 식을 따른다.

 

$N(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x - \mu)^{2}}{2\sigma^{2}})$

 

여기서 N개의 관측 데이터로 구성된 샘플 D = {x1, x2, x3 ...}가 있다고 할 경우 확률분포는 다음과 같이 표현할 수 있다.

 

$p(D; \mu, \sigma) = \prod_{n=1}^{N}\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x^{(n)} - \mu)^{2}}{2\sigma^{2}})$

 

여기서 이대로 식을 쓰지 않고 위에서 말했다시피 로그 가능도를 계산한다고 하면 여기에 log를 취한 다음 계산을 한다. 

 

$log p(D; \mu, \sigma) = log \prod_{n=1}^{N}\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x^{(n)} - \mu)^{2}}{2\sigma^{2}}) $

$= log\prod_{n=1}^{N} \frac{1}{\sqrt{2\pi}\sigma} + log \prod_{n=1}^{N} exp(-\frac{(x^{(n)} - \mu)^{2}}{2\sigma^{2}}) $

$= log(\frac{1}{\sqrt{2\pi}\sigma})^{N} + \sum_{n=1}^{N}log exp(-\frac{(x^{(n)}-\mu)^{2}}{2\sigma^{2}})$

$= log(\frac{1}{\sqrt{2\pi}\sigma})^{N} + \sum_{n=1}^{N} -\frac{(x^{(n)}-\mu)^{2}}{2\sigma^{2}} $ <여기서 log는 밑이 e이고 exp는 e의 지수를 뜻하기에 log를 없앨 수 있다.>

$= -\frac{N}{2}log2\pi\sigma^{2} - \frac{1}{2\sigma^{2}}\sum_{n=1}^{N}(x^{(n)}-\mu) ^{2}$

 

이렇게 된다. 이제 식은 깔끔하게 정리 되었으니 미분을 하면 된다. 여기서 매개변수가 뮤와 시그마가 있는데 각각 구해야 하므로 편미분을 이용한다. 편미분은 간단하게 설명하면 두 변수 중 하나에 대해서만 미분하고 나머지 변수는 상수 취급을 하여 미분하는 방식이다.

 

$L(\mu, \sigma) = log p(D; \mu, \sigma)$

일 때 미분 식은 다음과 같다 <학교 혹은 학원에서 배운 겉미분 속미분 그거...>

 

$\frac{\partial L}{\partial \mu} = \frac{1}{\sigma^{2}} \sum_{n=1}^{N}(x^{(n)}-\mu)$

 

위 식에서 미분 값 즉 $\frac{\partial L}{\partial \mu}$ 가 0이 되는 지점이 기댓값의 최댓값이 된다. 식을 전개하면 다음과 같다. 

$\frac{1}{\sigma^{2}}\sum_{n=1}^{N}(x^{(n)}-\mu) = 0$  <앞에 식은 시그마에 대한 식으로 사실상 상수이므로 미분하면 없어지니 걍 날림>

$\Leftrightarrow \sum_{n=1}^{N}(x^{(n)} - \mu) = 0$

$\Leftrightarrow \sum_{n=1}^{N}x^{(n)} = \sum_{n=1}^{N}\mu$

$\Leftrightarrow \sum_{n=1}^{N}x^{(n)} = N\mu$

$\Leftrightarrow frac{1}{N}\sum_{n=1}^{N} = \mu$

 

이렇게 해서 뮤 값이 저 위의 값일 때 최댓값을 가지는 것을 알 수 있었다. 이 최댓값을 가지는 뮤 값을 다음과 같이 표현할 수 있다.

 

$\hat{\mu} = \frac{1}{N}\sum_{n=1}^{N}x^{(n)}$

 

다음은 표준편차인 시그마에 대한 편미분 차례이다. 여기서 구할 것은 뮤 값이 최댓값을 가지는 뮤 일 때 즉 뮤 햇일 때의 시그마 최댓값이므로 뮤에 뮤햇을 대입한 다음 미분을 진행해주면 간단하게 풀 수 있다.

 

$log p(D; \mu = \hat{\mu}, \sigma) = -\frac{1}{2N}log2\pi\sigma^{2} - \frac{1}{2\sigma^{2}}\sum_{n-1}^{N}(x^{(n)}-\hat{\mu})^{2}$

 

일단 앞에 식부터 미분해보자

 

$\frac{\partial}{\partial \sigma}log(2\pi\sigma^{2}) = \frac{\partial}{partial\sigma}log(2\pi) + log(\sigma^{2})$

$= \frac{\partial}{\partial\sigma}log(2\pi) + \frac{\partial}{\partial\sigma}log(\sigma^{2})$

$= \frac{\partial}{\partial\sigma}log(\sigma^{2})$

$= \frac{\partial}{\partial\sigma}2log\sigma$

$= \frac{2}{\sigma}$

 

그럼 뒤 식과 이어서 전체 미분을 해보면

 

$- \frac{N}{2}\frac{2}{\sigma} - \frac{\partial}{\partial\sigma}(\frac{1}{2\sigma^{2}}\sum_{n=1}^{N}(x^{(n)} - \hat{\mu})^{2})$

$= -\frac{N}{\sigma}+\frac{1}{\sigma^{3}}\sum_{n=1}^{N}(x^{(n)} - \hat{\mu})^{2}$

 

위 값이 0인 시그마일 때 최댓값을 가지기에 다음과 같은 식을 짤 수 있다.

 

$\frac{N}{\sigma} = \frac{1}{\sigma^{3}}\sum_{n=1}^{N}(x^{(n)}-\hat{\mu})^{2}$

$\Leftrightarrow \sigma^{2} = \frac{1}{N}\sum_{n=1}^{N}(x^{(n)} - \hat{\mu})$

$\therefore \sigma = \sqrt{\frac{1}{N}\sum_{n=1}^{N}(x^{(n)} - \hat{\mu})}$

 

위 식을 계산하여 위에서 했던 매개변수 추론이 가능한 것이다. 

2.4 생성모델의 용도

이제 생성모델로 기존 키 데이터를 통해 생성모델을 구현해보자 먼저 이전에 배웠던 생성모델의 작업은 아래 2개이다.

 

1. 모델링: 주어진 데이터의 분포를 결정 여기서는 정규분포

2. 매개변수 추론: 최대가능도 추정으로 매개변수 추론

 

코드화 하면 다음과 같다.

import numpy as np
import matplotlib.pyplot as plt
import os

path =  os.path.join(os.getcwd(), "height.txt")
xs = np.loadtxt(path)

mu = np.mean(xs)
sigma = np.std(xs)
print(mu, sigma)

 

위 과정은 주어진 데이터로부터 기댓값과 표준편차를 계산한 다음 출력하는 코드이다. 따라서 매개변수 추론을 수행하는 것이고 이것을 다른말로 학습이라고도 한다.

2.4.1 새로운 데이터 생성

이제 학습이 되었으니 새로운 데이터를 생성하는 것도 가능하다. 여기서 사용될 함수는 

np.random.normal(loc=0.0, scale=1.0, size=None)

위 함수이다. 여기서 loc은 평균(기댓값), scale은 표준편차, size는 생성할 값의 갯수이다. 지정하지 않으면 1개가 된다. 이제 데이터를 생성하고 이것또한 그래프로 표현해보는 코드를 작성해보자.

import numpy as np
import matplotlib.pyplot as plt
import os

path =  os.path.join(os.getcwd(), "height.txt")
xs = np.loadtxt(path)

mu = np.mean(xs)
sigma = np.std(xs)
samples = np.random.normal(loc=mu, scale=sigma, size=10000)

plt.hist(xs, bins="auto", density=True, alpha=0.7, label="original")
plt.hist(samples, bins="auto", density=True, alpha=0.7, label="generated")
plt.xlabel("Height")
plt.ylabel("Probability Density")
plt.legend()
plt.show()

위 코드는 hist 함수를 통해 2개의 데이터를 다른 label로 표현을 하고 하나의 그래프에 표현되도록 했다. 

위 결과를 보면 새롭게 생성한 데이터도 실제 데이터의 정규분포를 잘 따르는 것을 알 수 있다. 

위 결과를 통해 알 수 있는 것은 생성 모델이란 것은 데이터의 특징<여기서는 기댓값, 표준편차>을 잘 파악한다면 기존 데이터와 비슷한 형태의 데이터를 생성할 수 있다는 것이다. 

2.4.2 확률 계산

이제 생성 모델로부터 어떤 값이 얼마나 발생하기 쉬운지 확인해볼 차례이다. 확률을 구한다는 의미인데 정규분포는 연속확률분포이므로 단순 넓이 계산을 적분을 통해 해야한다. 정규분포는 식이 정해져있기 때문에 적분을 하면 되고, 해석적으로 적분이 어려운 경우 몬테카를로 방법을 통해 근사적으로 구할 수 있다. 

from scipy.stats import norm

x=1.0
p = norm.cdf(x, loc=0, scale=1)

print(p)

위 코드를 이용하면 기댓값이 0이고 표준편차가 1인 정규분포에서 x가 1이하인 구간에서의 적분 결과를 보여준다. 따라서 위 함수를 이용하여 이전에 사용한 키 데이터에서의 적분을 할 수 있다.  

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
import os

path =  os.path.join(os.getcwd(), "height.txt")
xs = np.loadtxt(path)

mu = np.mean(xs)
sigma = np.std(xs)
samples = np.random.normal(loc=mu, scale=sigma, size=10000)

p = norm.cdf(180, loc=mu, scale=sigma)

print(1-p)

위 코드를 이용한다면 키 데이터에서 180보다 큰 경우에 대한 확률을 구할 수 있다. print할 때 1을 뺀 이유는 p값은 180 이하에 대한 확률이기 때문에 전체 확률인 1에서 p를 빼면 180보다 큰 경우에 대한 확률을 구할 수 있기 때문이다.