다운로드
작성자: admin 작성일시: 2016-04-29 19:23:57 조회수: 5390 다운로드: 366
카테고리: 기초 수학 태그목록:

최대 가능도 모수 추정의 원리

가능도 함수

확률변수 $X$에 대한 확률모형은 확률밀도함수 $p(x;\theta)$(또는 확률질량함수 $P(x;\theta)$)에 의해 정의된다. 여기에서 $x$는 확률변수가 가질 수 있는 실수값이고 $\theta$는 확률밀도함수 즉, 확률모형의 모수(parameter) 집합을 대표하는 기호이다.

예를 들어 가우시안 정규 분포 확률변수 $X$의 확률밀도함수는 다음과 같다. 이 식에서 $\theta = (\mu, \sigma^2)$ 이다.

$$ p(x; \theta) = p(x; \mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left({-\frac{(x-\mu)^2}{2\sigma^2}}\right) $$

함수의 관점에서 일반적으로 $\theta$는 고정된 값, 즉 상수 계수이고 $x$를 변수(variable)로 가정한다. 즉, 이미 확률변수 모형은 고정되어 있고 주어진 실수 입력값에 대해 그 실수값이 나올 상대적 가능성을 출력하는 것이 $x$를 변수로 가지는 확률밀도함수이다.

그러나 반대로 추정 문제에서는 $x$ 즉, 이미 실현된 표본값은 알고 있지만 모수 $\theta$를 모르고 있다. 이 경우에는 확률밀도함수라는 수식에서 $x$를 이미 결정되어 버린 상수 계수로 놓고 $\theta$를 찾아내야 할 변수로 생각할 수 있다. 물론 $f(x, \theta)$의 의미 자체는 변함없이 주어진 $x$가 나올 수 있는 확률밀도 값이다. 이를 가능도 함수라고 한다.

가능도 함수를 수식으로 나타내면 수식 자체는 확률밀도함수의 수식과 다르지 않다. 모수를 상수가 아닌 확률변수로 보는 경우에는 조건부 확률밀도함수로 쓸 수도 있다.

$$ \mathcal{L}(\theta;x) = p(x ; \theta) = p(x \vert \theta) $$

확률밀도함수는 모든 표본값에 대해 적분하면 전체 면적이 1이 되지만,

$$ \int_{-\infty}^{\infty} p(x; \theta) dx = 1 $$

가능도 함수는 확률밀도함수처럼 적분하였을 때 1이 된다는 보장이 없다.

$$ \int_{-\infty}^{\infty} p(x; \theta) d\theta \neq 1 $$
가능도
  • 가능도 함수
    • 확률밀도함수를 랜덤변수의 값 $x$의 함수가 아닌 파라미터 $\theta$의 함수로 보는 것
    • 확률 분포로부터 특정한 표본 값 $x$가 발생하였을 때, 이 표본 값 $x$가 나오게 하는 파라미터 $\theta$의 가능성
    • 확률 분포로부터 특정한 표본 값 $x$가 발생하였을 때, 표본 값 $x$와 변수 $\theta$에서의는 확률(밀도함수)
가능도 함수와 확률밀도함수
  • 확률밀도함수 $f(x; \theta) $
    • $\theta$ 값을 이미 알고 있음
    • $\theta$는 상수, $x$는 변수
    • $\theta$가 이미 정해져 있는 상황에서의 $x$ 값의 상대적 확률
    • 적분하면 전체 면적은 항상 1
  • 가능도 함수 $L(\theta) = p(x|\theta)$
    • $x$가 이미 발생. 값을 이미 알고 있음
    • $x$는 상수, $\theta$는 변수
    • $x$가 이미 정해져 있는 상황에서의 $\theta$ 값의 상대적 확률
    • 적분하면 전체 면적이 1이 아닐 수 있다.

예를 들어 정규 분포 함수에서 기댓값 모수와 분산 모수를 입력 변수로 가지는 가능도 함수를 그리면 각각 다음과 같다. 기댓값 모수를 입력 변수로 가지는 가능도 함수의 모양이 확률밀도함수와 같은 모양인 것은 ($x$와 $\mu$를 바꾸어도 식이 같아지는) 정규 분포의 확률밀도함수가 가지는 특별한 성질 때문이며 일반적이지 않다.

In [1]:
def likelihood_mu(mu):
    return sp.stats.norm(loc=mu).pdf(0)

mus = np.linspace(-5, 5, 1000)
likelihood_mu = [likelihood_mu(m) for m in mus]


plt.subplot(211)
plt.plot(mus, likelihood_mu)
plt.title("가능도 함수 $L(\mu)=N(x;\mu,\sigma^2)$")
plt.xlabel("$\mu$")
plt.show()

def likelihood_sigma2(sigma2):
    return sp.stats.norm(scale=np.sqrt(sigma2)).pdf(0)

sigma2s = np.linspace(0.1, 10, 1000)
likelihood_sigma2 = [likelihood_sigma2(s) for s in sigma2s]

plt.subplot(212)
plt.plot(sigma2s, likelihood_sigma2)
plt.title("가능도 함수 $L(\sigma^2)=N(x;\mu,\sigma^2)$")
plt.xlabel("$\sigma^2$")
plt.show()

최대 가능도 추정

최대 가능도 추정(MLE: Maximum Likelihood Estimation) 방법은 주어진 표본 $x$에 대해 가능도를 가장 크게 해 주는 모수 $\theta$를 찾는 방법이다.

예를 들어 정규 분포를 가지는 확률변수의 분산 $\sigma^2=1$은 알고 있으나 평균 $\mu$를 모르고 있어 이를 추정해야 하는 문제를 생각해 보자. 확률변수의 표본은 하나 $x_1=1$를 가지고 있다고 하자. 이 경우 어떤 $\mu$ 값이 가장 가능성(가능도)이 커 보이는가? 다음 그림에는 $\mu=-1$, $\mu=0$, $\mu=1$,세가지 후보를 제시한다. 이 세가지 $\mu$ 값에 대해 $1$이 나올 확률밀도의 값이 바로 가능도이다.

In [2]:
x = np.linspace(-5, 5, 100)

p1 = sp.stats.norm(loc=-1).pdf(1)
p2 = sp.stats.norm(loc=0).pdf(1)
p3 = sp.stats.norm(loc=1).pdf(1)

plt.scatter(1, p1, s=100, c='r')
plt.scatter(1, p2, s=100, c='b')
plt.scatter(1, p3, s=100, c='g')

plt.plot(x, sp.stats.norm(loc=0).pdf(x), ls="-.", 
         label=r"$N(x;\mu=-1)$={}".format(np.round(p1, 2)))
plt.plot(x, sp.stats.norm(loc=-1).pdf(x), ls="--", 
         label=r"$N(x;\mu=0)$={}".format(np.round(p2, 2)))
plt.plot(x, sp.stats.norm(loc=1).pdf(x), ls="-", 
         label=r"$N(x;\mu=1)$={}".format(np.round(p3, 2)))
plt.scatter(1, 0, s=100, c='k')
plt.vlines(1, -0.09, 0.45, linestyle=":")
plt.text(1-0.3, -0.15, "$x_1=1$")

plt.legend()
plt.title("최대 가능도 추정의 원리")
plt.show()
  • $N(x; \mu=-1)$이라는 확률변수에서 $x=1$이 나올 가능도(확률밀도)는 0.05이다.
  • $N(x; \mu=0)$이라는 확률변수에서 $x=1$이 나올 가능도(확률밀도)는 0.24이다.
  • $N(x; \mu=1)$이라는 확률변수에서 $x=1$이 나올 가능도(확률밀도)는 0.4이다.

어떤 확률변수를 고르는 것이 합리적인가? 당연히 가장 큰 가능도를 가진 확률변수를 선택해야 한다. 그림에서 볼 수 있듯이 $\mu=1$일 경우의 가능도가 가장 크다. 따라서 최대 가능도 추정에 의한 추정값은 $\hat\mu_{\text{ML}}=1$이다.

일반적으로는 추정을 위해 확보하고 있는 확률변수 표본의 수가 하나가 아니라 복수개 $\{ x_1, x_2, \cdots x_N \}$ 이므로 가능도도 복수 표본값에 대한 결합 확률밀도 $p_{X_1, X_2, \cdots, X_N}(x_1, x_2, \cdots, x_N ; \theta)$ 에서 구해야 한다. 표본 데이터 $x_1, x_2, \cdots x_N$는 같은 확률 분포에서 나온 독립적인 값들이므로 결합 확률밀도함수는 다음처럼 곱으로 표현된다.

$$ L(\theta; x_1, \ldots, x_N) = L(\theta; \{ x_i \}_{i=1}^N) = \prod_{i=1}^N p(x_i; \theta) $$

예를 들어 가우시안 정규 분포로 부터 얻은 표본값이 $x_1 = 1, x_2 = 0, x_3 = -3$이었다면 이 경우의 가능도 함수는 다음과 같다.

$$ \begin{eqnarray} L(\theta; x_1, x_2, x_3) &=& \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left({-\frac{(1-\mu)^2}{2\sigma^2}}\right) \cdot \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left({-\frac{(0-\mu)^2}{2\sigma^2}}\right) \cdot \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left({-\frac{(-3-\mu)^2}{2\sigma^2}}\right) \\ &=& \dfrac{1}{(2\pi\sigma^2)^{\frac{3}{2}}} \exp\left({-\frac{\mu^2 + (1-\mu)^2 + (-3-\mu)^2}{2\sigma^2}}\right) \\ &=& \dfrac{1}{(2\pi\sigma^2)^{\frac{3}{2}}} \exp\left({-\frac{3\mu^2+4\mu+10}{2\sigma^2}}\right) \cdot \end{eqnarray} $$
In [3]:
mu = np.linspace(-3, 3, 1000)
sigma2 = 1

def likelihood(mu):
    return (2 * np.pi * sigma2) ** (3 / 2) * np.exp(-(3 * mu ** 2 + 4 * mu + 10) / (2 * sigma2))

li = likelihood(mu)

plt.plot(mu, li)
plt.xlabel(r"$\mu$")
plt.title("데이터가 1, 0, -3이 나온 경우의 정규 분포의 가능도 함수")
plt.show()

최대 가능도 추정의 구현

실제로 최대 가능도 추정 방법을 사용하려면 가능도가 최대가 되는 $\theta$를 수치적으로 계산해야 한다. 즉 수치적 최적화(numerical optimization) 문제가 된다.

$$ \hat\theta_{\text{ML}} = \arg \max_{\theta} L(\theta; \{x_i\}) $$

앞에서 예를 든 것처럼 정규분포에서 데이터가 $\{1, 0, -3\}$이 나왔을 경우에는 $\mu=-2/3$일 때 최대값을 가진다.

일반적으로 가능도를 직접 사용하는 것이 아니라 다음과 같은 이유로 로그 변환한 로그 가능도 함수 ${LL} = \log{{L}}$를 사용하는 경우가 많다.

$$ \hat\theta_{\text{ML}} = \arg \max_{\theta} \log{L}(\theta; \{x_i\}) $$

이유는 다음과 같다.

  1. 로그 변환에 의해서는 최대값의 위치가 변치 않는다
  2. 반복시행으로 인한 복수 표본 데이터인 경우 결합 확률밀도함수가 동일한 함수의 곱으로 나타나는 경우가 많은데 이 때 로그 변환에 의해 곱셈이 덧셈이 되어 계산이 단순해진다.

앞에서 든 예에서도 로그 변환을 하면 최대값의 위치가 $-2/3$라는 것을 쉽게 구할 수 있다.

$$ \begin{eqnarray} \log L(\mu; x_1, x_2, x_3) &=& \log \left( \dfrac{1}{(2\pi\sigma^2)^{\frac{3}{2}}} \exp\left({-\frac{3\mu^2+4\mu+10}{2\sigma^2}}\right) \right) \\ &=& \log \left( \dfrac{1}{(2\pi\sigma^2)^{\frac{3}{2}}} \right) -\frac{3\mu^2+4\mu+10}{2\sigma^2} \\ &=& \log \left( \dfrac{1}{(2\pi\sigma^2)^{\frac{3}{2}}} \right) -\frac{3\left(\mu+\frac{2}{3}\right)^2+\frac{26}{4}}{2\sigma^2} \\ \end{eqnarray} $$

연습 문제 10.2.1

베르누이 분포 확률변수로부터 다음과 같은 표본을 얻었다. 이 확률변수의 모수 $\mu$를 최대 가능도 추정법을 사용하여 구하라.

$$ \{ 1, 0, 1, 1 \} $$

연습 문제 10.2.2

$K=4$인 카테고리 분포 확률변수로부터 다음과 같은 표본을 얻었다. 이 확률변수의 모수 $\mu$를 최대 가능도 추정법을 사용하여 구하라.

$$ \{ 1, 4, 1, 2, 4, 2, 3, 4 \} $$

연습 문제 10.2.3

2차원 다변수 가우시안 정규 분포 확률변수로부터 다음과 같은 표본을 얻었다. 이 확률변수의 모수 $\mu$, $\Sigma$를 최대 가능도 추정법을 사용하여 구하라.

$$ x_1 = (-1, 2), \; x_2 = (0, 2) , \; x_3 = (-2, -1) $$

최대 가능도 추정의 신뢰 구간

가능도 함수가 최대가 되는 값 즉 모수 추정치를 둘러싼 구간을 신뢰 구간(confidence interval)이라고 한다. 가우시안 정규 분포의 기댓값 모수에 대한 가능도 함수는 확률밀도함수처럼 적분하면 1이 되므로 이러한 경우에는 확률값을 가진 신뢰 구간을 계산할 수 있다.

질문/덧글

likelihood_sigma2 함수 chan*** 2018년 10월 30일 11:36 오전

함수 likelihood_sigma2(theta)에서 theta는 분산이니까
sp.stats.norm(**scale=theta**).pdf(0)에서 sp.stats.norm(**scale=np.sqrt(theta)**).pdf(0) 바뀌어야 하지 않을까요?

답변: likelihood_sigma2 함수 관리자 2018년 11월 1일 11:20 오전

지적 감사합니다. 수정하였습니다.