[Boostcamp 선택 과제 - 3] Maximum Likelihood Estimation (MLE)

우선 최대 가능도 추정에 대한 계산과 내용은 이전에 작성한 통계론을 확인하면된다.

최대가능도추정의 실습

$x_0 = 1$라는 샘플을 하나 가지고 있고 $x_0 = 1$일 때, 어떤 평균(모수)에 대해서 가장 큰 가능도를 갖는지 실습하는 문제이다. 그리고 $\mu = -1, 0, 1$ 이렇게 세 가지의 경우 정규분포에 대해 그래프를 그려보는 것이다.

Import

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.stats import norm #연속확률밀도함수를 계산하는 함수이다.

plot

plt.figure(figsize=(10,6))

x = np.linspace(-5, 5, 100) # x축 정의

p1 = sp.stats.norm(loc=-1).pdf(1) # 모수가 -1일 때 1에서의 가능도
p2 = sp.stats.norm(loc=0).pdf(1) # 모수가 0일 때 1에서의 가능도
p3 = sp.stats.norm(loc=1).pdf(1) # 모수가 1일 때 1에서의 가능도

# 산점도
plt.scatter(1, p1, s=100, c='blue', marker='v', 
         label=r"$N(x_1;\mu=-1)$={:.2f}".format(np.round(p1, 2)))
plt.scatter(1, p2, s=100, c='orange', marker='^', 
         label=r"$N(x_1;\mu=0)$={:.2f}".format(np.round(p2, 2)))
plt.scatter(1, p3, s=100, c='green', marker='s', 
         label=r"$N(x_1;\mu=1)$={:.2f}".format(np.round(p3, 2)))

plt.plot(x, norm(loc=-1).pdf(x), ls="-.") # 모수가 -1일 때 x값에 따른 정규분포 그래프
plt.plot(x, norm(loc=0).pdf(x), ls="--") # 모수가 0일 때 x값에 따른 정규분포 그래프
plt.plot(x, norm(loc=1).pdf(x), ls="-") # 모수가 1일 때 x값에 따른 정규분포 그래프

plt.scatter(1, 0, s=100, c='k')
plt.vlines(1, -0.09, 0.45, linestyle=":")
plt.hlines(0, -5, 5, colors='gray', linestyle="-")
plt.text(1-0.3, -0.15, "$x_0=1$")

plt.xlabel("x")
plt.ylabel("likelihood")
plt.legend()
plt.title("MLE for population mean")
plt.show()

MLE 계산

x0=1

# 모수에 따라 x가 1일 때 가능도를 계산한다.
print('mu=-1: likelihood at x_0=1 is {:.4f}'.format(norm.pdf(x0, -1)))
print('mu=0: likelihood at x_0=1 is {:.4f}'.format(norm.pdf(x0, 0)))
print('mu=1: likelihood at x_0=1 is {:.4f}'.format(norm.pdf(x0, 1)))

[출력]
mu=-1: likelihood at x_0=1 is 0.0540
mu=0: likelihood at x_0=1 is 0.2420
mu=1: likelihood at x_0=1 is 0.3989

따라서 모수가 1일 때 최대 가능도 약 0.4가 나오는 것을 확인 할 수 있다.

정규분포의 가능도 함수

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

plt.figure(figsize=(12,10))

# 기댓값 모수 mu일 때 0에서의 가능도를 계산
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(311)

# 기댓값 모수 = 0, 분산 모수 = 1
mu = 0
variance = 1

sigma = math.sqrt(variance) #표준편차

plt.plot(mus, stats.norm.pdf(mus, mu, sigma))
plt.title("Probability Density Function of Normal $P(x; \mu=0, \sigma^2=1)$")
plt.xlabel("$x$")
plt.show()

plt.figure(figsize=(12,10))
plt.subplot(312)
plt.plot(mus, likelihood_mu)
plt.title("Likelihood $L(\mu, \sigma^2=1; x=0)$")
plt.xlabel("$\mu$")
plt.show()

# 분산 모수 값에 따른 0에서의 가능도 계산
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.figure(figsize=(12,10))

plt.subplot(313)
plt.plot(sigma2s, likelihood_sigma2)
plt.title("Likelihood $L(\mu=0, \sigma; x=0)$")
plt.xlabel("$\sigma^2$")
plt.show()

 

 

3차원 공간에서의 시각화

MU, SIGMA2 = np.meshgrid(mus, sigma2s)
L = np.exp(-MU ** 2 / (2 * SIGMA2)) / np.sqrt(2 * np.pi * SIGMA2)

fig = plt.figure(figsize=(10,10))
ax = fig.gca(projection='3d')
ax.plot_surface(MU, SIGMA2, L, linewidth=0.1)
plt.xlabel('$\mu$')
plt.ylabel('$\sigma^2$')
plt.title('likelihood $L(\mu, \sigma^2)$ in 3D')
plt.show()

베르누이분포의 최대가능도

모수가 $\mu$인 베르누이 분포의 확률질량함수는 다음과 같다.
$$p(x; \mu) = Bern(x; \mu) = \mu^x(1 - \mu)^{1-x}$$
표본 데이터가 모두 독립일 때 전체 확률질량함수는 각각의 확률질량함수의 곱과 같으므로 다음과 같은 식이 만들어진다.

$$\prod_{i=1}^n \mu^x_i(1 - \mu)^{1-x_i}$$
이 때의 모수 $\mu$를 추정해보는 문제이다.

 

일단 log를 씌워서 정리하면 다음과 같은 식이 될 것이다.
$$\sum_{i=1}^n x_i \log \mu + (1 - x_i) \log{(1 - \mu)}$$

그리고 여기서 성공 횟수와 실패 획수를 각각 $N_0$, $N_1$이라고 표기하고 대입하면 $\sum_{i=1}^n N_1 \log \mu + N_0 \log{(1 - \mu)}$가 될 것이다.

 

이제 정리된 식을 미분하면 $\frac{N_1}{\mu} - \frac{N_0}{1-\mu} = 0$이 될 것이고 최종적으로 정리하면 $\mu = \frac{N_1}{N}$이 된다.($N_1 + N_0 = N$이기 때문)