최대가능도 추정법

#1. 최대가능도 추정법의 원리

치명타 확률이 $p$인 무기로 일정시간 동안 고정된 타깃을 공격한다고 하자. 타깃을 1번 공격하고 치명타가 발생하기를 기대하는 것은 1번의 베르누이 시행과 같다. 타깃을 매번 공격할 때마다, 치명타의 발생 여부는 0 또는 1의 값을 가지는 확률변수이며, 모수가 $p$인 베르누이 분포를 따른다.

공격 종료후 총 피해량을 측정했더니 꽤 높게 나왔다고 하자. 이때, $p$의 값을 높다고 추정하는 것이 합리적이다. $p$의 값이 낮다면, 이만한 피해를 줄 수 없기 때문이다. 반면, 예상보다 낮은 피해량이 측정되었다면 $p$의 값이 낮다고 추정하는 것이 합리적이다.

모수를 추정하는 방법 중 최대가능도 추정법은 바로 이런 원리에 기반한 방법이다. 모집단 $X$로부터 추출한 표본의 관측값이 $x_1,\cdots,x_n$이라고 하자. 경험에 의거해 이런 관측값을 가지는 것은 모집단의 모수가 $\theta$이기 때문일 것이라고 생각한다. $\theta$가 아니라면, 이런 관측값을 가질리가 없기 때문이다.

#2. 최대가능도 추정량

이제 이 모수 $\theta$를 표본의 관측값 $x_1,\cdots,x_n$으로 추정하려고 한다. 표본의 관측값은 이미 일어난, 연구자가 이미 관측한 사실이다. 모수가 $\theta$라는 조건 하에 이런 표본을 얻을 확률을 가능한 크게 만드는 추정량이 합리적인 추정량이다.

$$ \hat{\theta}=\arg\max P(X_1=x_1,\cdots,X_n=x_n\mid \theta) $$

식에서 $P(X_1=x_1,\cdots,X_n=x_n\mid \theta)$는 베이즈 추정에서 가능도(Likelihood; 似然)라고 부른다. 따라서 이 식을 만족하는 추정량을 최대가능도 추정량(Maximum Likelihood Estimator, MLE; 最大似然估计量)이라고 한다. 가능도를 최대로 만드는 추정량이라는 뜻이다. 최대가능도 추정량을 찾는 방법이 바로 최대가능도 추정법이다.

note: 베이즈 추정은 가능도가 아닌 사후확률을 최대화하는 모수추정이라고 볼 수 있다.

#3. 가능도 함수

가능도의 식이 매우 길어서 간단히 할 필요가 있다. 그리고 가능도의 식에서 $X_i=x_i$는 이미 일어난 사실이고, 미지수는 $\theta$이다. 따라서 $\theta$에 관한 함수임을 강조하기 위해, 다음과 같이 가능도 함수(Likelihood Function; 似然函数)를 정의한다.

$$ L(\theta):=P(X_1=x_1,\cdots,X_n=x_n\mid \theta) $$

수리통계학의 단순무작위 추출 가정 덕분에 표본의 개체 $X_i$는 서로 독립이며, 결합확률은 개별확률의 곱으로 나타낼 수 있다. 따라서 가능도함수를 다음과 같이 변형할 수 있다.

$$ L(\theta)=\prod_{i=1}^{n}P(X_i=x_i\mid \theta) $$

한편 $P(X_i=x_i\mid \theta)$는 $\theta$라는 조건 하에 확률변수 $X_i$가 $x_i$라는 값을 가질 확률질량이다. 먼저 모집단 $X$가 이산형일 때, 가능도 함수는 확률질량함수 $f(x_i\mid\theta)$의 곱이다.

$$ L(\theta)=\prod_{i=1}^{n}f(x_i\mid\theta) $$

모집단 $X$가 연속형이라면, 확률질량을 표현하기 위해 먼저 $X_i$가 $x_i$ 근처의 값을 가질 확률밀도 $f(x_i\mid\theta)$를 생각한 후, 거기에 미소구간의 길이 $dx_i$를 곱해야 한다.

$$ L(\theta)=\prod_{i=1}^{n}f(x_i\mid\theta)\space dx_i=\prod_{i=1}^{n}f(x_i\mid\theta)\prod_{i=1}^{n}dx_i $$

그런데 어차피 $\prod_{i=1}^{n}dx_i$는 $\theta$와는 관련 없는 상수이기 때문에, 가능도 함수를 최대화 하는데는 영향이 없다. 따라서 따로 떼어버리고 가능도 함수를 이산형과 마찬가지의 형태로 생각하기로 한다.

#4. 로그 가능도

이제 가능도 함수 $L(\theta)=\prod_{i=1}^{n}f(x_i\mid\theta)$가 정해졌으므로, 남은 것은 이 함수를 최대화하는 $\hat\theta=\arg\max L(\theta)$를 찾는 것이다.

함수를 최대화하는 값을 찾는 법은 함수의 미분이 0이라는 방정식을 세운 뒤, 푸는 것이다. 따라서 $L(\theta)=\prod_{i=1}^{n}f(x_i\mid\theta)$를 미분하는 것을 생각하게 된다.

그런데 누율생성함수의 탄생계기에서도 알 수 있듯이, 곱셈은 미분친화적이지 않다. 가능도 함수도 적률생성함수에 로그를 씌워 누율생성함수를 만든것처럼 로그를 씌워서 해결하면 된다. 이를 로그 가능도(Log Likelihood; 对数似然)라고 한다.

로그 함수는 단조증가 함수이므로 가능도를 최대화하는 문제는 곧 로그 가능도를 최대화하는 문제로 바꿀 수 있다. 따라서 MLE는 최종적으로 로그 가능도를 최대화하는 추정량으로 정의할 수 있고

$$ \hat{\theta}=\arg\max\ln L(\theta) $$

이 추정량을 구하고 싶으면 다음과 같은 방정식을 풀면 된다.

$$ \frac{d}{d\theta}\ln L(\theta)=0 $$

모집단의 모수가 하나라면 모수는 $\theta$와 같이 나타낸다. 그런데 모수가 만약 여러개($k$개)라면 모수와 MLE를 다음과 같이 벡터의 형태로 나타내야 하고,

$$ \boldsymbol{\theta}=(\theta_1,\cdots,\theta_k),\space\boldsymbol{\hat{\theta}}=\arg\max\ln L(\boldsymbol{\theta}) $$

MLE는 다음과 같은 연립방정식을 풀어서 구할 수 있다.

$$ \left\{ \begin{alignat*}{4} \frac{\partial}{\partial\theta_1}\ln L(\boldsymbol{\theta}) &=& 0 \\ \frac{\partial}{\partial\theta_2}\ln L(\boldsymbol{\theta}) &=& 0 \\ \vdots \\ \frac{\partial}{\partial\theta_k}\ln L(\boldsymbol{\theta}) &=& 0 \end{alignat*} \right. $$

#5. 예시

#5.1 베르누이 분포의 최대가능도 추정

위에서 예로든 무기의 치명타 확률 $p$를 추정해보자. 치명타의 발생여부 $X\sim\mathrm{Bern}(p)$이고, 확률질량은

$$ f(x)=p^x(1-p)^{1-x} $$

가능도 함수는

$$ L(p)=\prod_{i=1}^np^{x_i}(1-p)^{1-x_i}=p^{\sum_{i=1}^nx_i}(1-p)^{n-\sum_{i=1}^nx_i} $$

따라서 로그 가능도는

$$ \ln L(p)=(\sum_{i=1}^n x_i)\ln p+(n-\sum_{i=1}^n x_i)\ln(1-p) $$

$\frac{d}{dp}\ln L(p)=0$을 풀면 MLE를 구할 수 있다.

$$ \hat{p}=\frac{1}{n}\sum_{i=1}^{n}x_i=\overline{x} $$

표본평균이다. 지난 글에서 구한 베르누이 분포의 적률 추정량도 표본평균이었다. 따라서 베르누이 분포의 MLE와 MME는 같다.

베르누이 분포에서 표본평균이 의미하는 것은 베르누이 시행의 성공빈도다. 즉, 치명타 발생빈도다. 데이터를 수집한 뒤 치명타 발생 빈도가 높다면, 치명타 확률이 높다고 추정하는 것은 합리적이다.

큰 수의 법칙에 의해, 표본평균 $\overline{X}$는 모평균 $\mu$에 a.s. 수렴한다. 그리고 모평균 $\mu=p$이다. 따라서 MLE는 베르누이 분포의 모수에 a.s. 수렴한다. 다시 말해 따라서 데이터를 충분히 많이 수집한다면, 치명타 확률의 추정량(발생빈도)은 실제 치명타 확률에 가까워진다. 이는 우리의 직관과도 일치한다.

#5.2 정규분포의 최대가능도 추정

모집단 $X\sim\mathrm{N}(\mu,\sigma^2)$의 표본 $X_1,\cdots,X_n$을 추출한 상황이다.

모집단의 확률밀도는

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

로그 가능도는

$$ \ln L(\mu,\sigma^2)=-\frac{n}{2}\ln2\pi-\frac{n}{2}\ln\sigma^2-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2 $$

MLE를 구하기 위한 연립방정식은

$$ \left\{ \begin{alignat*}{4} \frac{\partial}{\partial\mu}\ln L(\mu,\sigma^2) &=& &\frac{1}{\sigma^2}(\sum_{i=1}^nx_i-n\mu) &=& 0 \\ \frac{\partial}{\partial\sigma^2}\ln L(\mu,\sigma^2) &=& &-\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^{n}(x_i-\mu)^2 &=& 0 \end{alignat*} \right. $$

이것을 풀면 다음과 같다.

$$ \begin{split} &\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}x_i=\overline{x} \\ &\hat{\sigma^2}=\frac{1}{n}\sum_{i=1}^{n}(x_i-\overline{x})^2 \end{split} $$

지난글에서 어떤 분포든지 모평균과 모분산의 MME는 표본평균과 표본의 2차 중심적률이라는 것을 확인했다. 지금 정규분포의 MLE도 이와 같다.

#5.3 균등분포의 최대가능도 추정

모집단 $X\sim\mathrm{U}[a,b]$의 표본 $X_1,\cdots,X_n$을 추출한 상황이다.

모집단의 확률밀도함수는 다음과 같다.

$$ f(x)= \begin{cases} \frac{1}{b-a},&a\le x\le b \\ 0,&\mathrm{otherwise} \end{cases} $$

가능도 함수는 다음과 같으며, 불연속점을 포함하게 된다. 따라서 미분법으로 가능도를 최대화하는 지점을 찾는 것은 바람직하지 않으므로 다른 방법을 찾아야 한다.

$$ L(a,b)= \begin{cases} \frac{1}{(b-a)^n},&a\le x_1,\cdots,x_n\le b \\ 0,&\mathrm{otherwise} \end{cases} $$

먼저 부등식 $a\le x_1,\cdots,x_n\le b$는 순서통계량을 이용하여 다음과 같이 간단히 할 수 있다.

$$ a\le x_{(1)},x_{(n)}\le b $$

그리고 이 부등식을 만족하는 $a$와 $b$에 대해, 다음이 성립한다.

$$ L(a,b)=\frac{1}{(b-a)^n}\le\frac{1}{(x_{(n)}-x_{(1)})^n} $$

즉, $a=x_{(1)}$이고 $b=x_{(n)}$일 때, 가능도 $L(a,b)$는 최대가 된다. 따라서 균등분포의 MLE는 다음과 같다.

$$ \begin{split} &\hat{a}=x_{(1)}=\min_{1\le i\le n}x_i \\ &\hat{b}=x_{(n)}=\max_{1\le i\le n}x_i \end{split} $$

note: 이 결과는 직관과 부합한다. 균등분포의 확률밀도가 0이 아닌 값의 구간을 추정할 때, 관찰된 데이터의 최대값과 최소값을 구간의 양 끝으로 잡는 것이다. 그런데 균등분포의 MLE는 지난 글에서 구한 MME와 너무 다르다. 생각해보면 둘 다 나름 이유가 있고 합리적인 추정법인데 다른 결과가 나오는 것을 어떻게 해석해야 할까?