다운로드
작성자: admin 작성일시: 2016-05-30 00:55:09 조회수: 3172 다운로드: 274
카테고리: 기초 수학 태그목록:

MLE 모수 추정의 예

베르누이 분포의 모수 추정

모수가 $\theta$인 베르누이 분포의 확률 질량 함수는 다음과 같다.

$$ P(x ; \theta ) = \text{Bern}(x ; \theta ) = \theta^x (1 - \theta)^{1-x}$$

그런데 $N$ 번의 반복 시행으로 표본 데이터가 $x_1, \cdots, x_N$ 가 있는 경우에는 모두 독립이므로 전체 확률 밀도 함수는 각각의 확률 질량 함수의 곱과 같다.

$$ L(\theta ; x_1, \cdots, x_N) = P(x_1, \cdots, x_N;\theta) = \prod_{i=1}^N \theta^{x_i} (1 - \theta)^{1-x_i} $$

미분을 쉽게 하기 위해 로그 변환을 한 Log-Likelihood 를 구하면 다음과 같다.

$$ \begin{eqnarray*} \log L &=& \log P(x_1, \cdots, x_N;\theta) \\ &=& \sum_{i=1}^N \big\{ {x_i} \log\theta + (1-x_i)\log(1 - \theta) \big\} \\ &=& \sum_{i=1}^N {x_i} \log\theta + \left( N-\sum_{i=1}^N x_i \right) \log( 1 - \theta ) \\ \end{eqnarray*} $$

$x = 1$(성공) 또는 $x= 0$ (실패) 이므로 성공 횟수는 다음과 같이 $N_1$이라고 표기할 수 있다.

$$ N_1 = \sum_{i=1}^N {x_i} $$

따라서 Log-Likelihood는 다음과 같아진다.

$$ \begin{eqnarray*} \log L &=& N_1 \log\theta + (N-N_1) \log(1 - \theta) \\ \end{eqnarray*} $$

이 목적 함수를 모수로 미분한 값이 0이 되게 하는 모수 값을 구하면 다음과 같다.

$$ \begin{eqnarray*} \dfrac{\partial \log L}{\partial \theta} &=& \dfrac{\partial}{\partial \theta} \big\{ N_1 \log\theta + (N-N_1) \log(1 - \theta) \big\} = 0\\ &=& \dfrac{N_1}{\theta} - \dfrac{N-N_1}{1-\theta} = 0 \\ \end{eqnarray*} $$$$ \dfrac{N_1}{\theta} = \dfrac{N-N_1}{1-\theta} $$$$ \dfrac{1-\theta}{\theta} = \dfrac{N-N_1}{N_1} $$$$ \dfrac{1}{\theta} - 1 = \dfrac{N}{N_1} - 1 $$$$ \theta= \dfrac{N_1}{N} $$
In:
np.random.seed(0)
theta0 = 0.6
x = sp.stats.bernoulli(theta0).rvs(1000)
N0, N1 = np.bincount(x, minlength=2)
N = N0 + N1
theta = N1/N
theta
Out:
0.60999999999999999

카테고리 분포의 모수 추정

모수가 $\theta = (\theta_1, \cdots, \theta_K)$인 카테고리 분포의 확률 질량 함수는 다음과 같다.

$$ P(x ; \theta_1, \cdots, \theta_K ) = \text{Cat}(x ; \theta_1, \cdots, \theta_K) = \prod_{k=1}^K \theta_k^{x_k} $$$$ \sum_{k=1}^K \theta_k = 1 $$

그런데 $N$ 번의 반복 시행으로 표본 데이터가 $x_1, \cdots, x_i, \cdots, x_N$ 가 있는 경우에는 모두 독립이므로 전체 확률 밀도 함수는 각각의 확률 질량 함수의 곱과 같다.

$$ L(\theta_1, \cdots, \theta_K ; x_1, \cdots, x_i, \cdots, x_N) = P(x_1, \cdots, x_i, \cdots, x_N;\theta_1, \cdots, \theta_K) = \prod_{i=1}^N \prod_{k=1}^K \theta_k^{x_{i,k}} $$

위 식에서 $x_{i,k}$는 $i$번째 시행 결과인 $x_i$의 $k$번째 원소를 뜻한다.

미분을 쉽게 하기 위해 로그 변환을 한 Log-Likelihood 를 구하면 다음과 같다.

$$ \begin{eqnarray*} \log L &=& \log P(x_1, \cdots, x_N;\theta_1, \cdots, \theta_K) \\ &=& \sum_{i=1}^N \sum_{k=1}^K \left( {x_{i,k}} \log\theta_k \right) \\ &=& \sum_{k=1}^K \sum_{i=1}^N \left( \log\theta_k \cdot {x_{i,k}}\right) \\ &=& \sum_{k=1}^K \left( \log\theta_k \left( \sum_{i=1}^N {x_{i,k}} \right) \right) \end{eqnarray*} $$

$x_k$가 나온 횟수를 $N_k$이라고 표기하자.

$$N_k = \sum_{i=1}^N {x_{i,k}}$$

그러면 Log-Likelihood 가 다음과 같아지며 이 함수를 최대화하는 모수의 값을 찾아야 한다.

$$ \begin{eqnarray*} \log L &=& \sum_{k=1}^K \left( \log\theta_k \cdot N_k \right) \end{eqnarray*} $$

그런데 모수는 다음과 같은 제한 조건을 만족해야만 한다.

$$ \sum_{k=1}^K \theta_k = 1 $$

따라서 라그랑주 승수법을 사용하여 Log-Likelihood 를 다음과 같이 변형한 후 모수로 미분한 값이 0이 되는 값을 구한다.

$$ \begin{eqnarray*} \dfrac{\partial \log L}{\partial \theta_k} &=& \dfrac{\partial}{\partial \theta_k} \left\{ \sum_{k=1}^K \log\theta_k N_k + \lambda \left(1- \sum_{k=1}^K \theta_k\right) \right\} = 0 & \;\;\; (k=1, \cdots, K) \\ \dfrac{\partial \log L}{\partial \lambda} &=& \dfrac{\partial}{\partial \lambda} \left\{ \sum_{k=1}^K \log\theta_k N_k + \lambda \left(1- \sum_{k=1}^K \theta_k \right) \right\} = 0 & \\ \end{eqnarray*} $$

이를 풀면 다음과 같이 모수를 추정할 수 있다.

$$ \dfrac{N_1}{\theta_1} = \dfrac{N_2}{\theta_2} = \cdots = \dfrac{N_K}{\theta_K} = \lambda $$$$ N_k = \lambda \theta_k $$$$ \sum_{k=1}^K N_k = \lambda \sum_{k=1}^K \theta_k = \lambda = N $$$$ \theta_k = \dfrac{N_k}{N} $$
In:
np.random.seed(0)
theta0 = np.array([0.1, 0.3, 0.6])
x = np.random.choice(np.arange(3), 1000, p=theta0)
N0, N1, N2 = np.bincount(x, minlength=3)
N = N0 + N1 + N2
theta = np.array([N0, N1, N2]) / N
theta
Out:
array([ 0.098,  0.317,  0.585])

정규 분포의 모수 추정

가우시안 정규 분포의 확률 밀도 함수는 다음과 같다.

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

그런데 $N$ 번의 반복 시행으로 표본 데이터가 $x_1, \cdots, x_N$ 가 있는 경우에는 모두 독립이므로 전체 확률 밀도 함수는 각각의 확류 밀도 함수의 곱과 같다.

$$ L(\theta;x_1, \cdots, x_N) = p(x_1, \cdots, x_N;\theta) = \prod_{i=1}^N \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp \left(-\dfrac{(x_i-\mu)^2}{2\sigma^2}\right)$$

미분을 쉽게 하기 위해 로그 변환을 한 Log-Likelihood 를 구하면 다음과 같다. 여기에서 상수 부분은 모아서 $C$ 로 표기하였다.

$$ \begin{eqnarray*} \log L &=& \log p(x_1, \cdots, x_N;\theta) \\ &=& \sum_{i=1}^N \left\{ -\dfrac{1}{2}\log(2\pi\sigma^2) - \dfrac{(x_i-\mu)^2}{2\sigma^2} \right\} \\ &=& -\dfrac{N}{2} \log(2\pi\sigma^2) - \dfrac{1}{2\sigma^2}\sum_{i=1}^N (x_i-\mu)^2 \end{eqnarray*} $$

이 확률 밀도 함수가 최대가 되는 모수 값을 찾기 위해서는 각각의 모수로 미분한 값이 0이 되어야 한다.

$$ \begin{eqnarray*} \dfrac{\partial \log L}{\partial \mu} &=& \dfrac{\partial}{\partial \mu} \left\{ \dfrac{N}{2} \log(2\pi\sigma^2) + \dfrac{1}{2\sigma^2}\sum_{i=1}^N (x_i-\mu)^2 \right\} = 0 \\ \dfrac{\partial \log L}{\partial \sigma^2} &=& \dfrac{\partial}{\partial \sigma^2} \left\{ \dfrac{N}{2} \log(2\pi\sigma^2) + \dfrac{1}{2\sigma^2}\sum_{i=1}^N (x_i-\mu)^2 \right\} = 0\\ \end{eqnarray*} $$

이 두 식을 풀면 주어진 데이터 표본에 대해 모수의 likelihood를 가장 크게 하는 모수의 값을 구할 수 있다. 먼저 $\mu$에 대한 미분을 정리하면 다음과 같다.

$$ \dfrac{2}{2\sigma^2}\sum_{i=1}^N (x_i-\mu) = 0 $$$$ N \mu = \sum_{i=1}^N x_i $$$$ \mu = \dfrac{1}{N}\sum_{i=1}^N x_i = \bar{x} $$

다음으로 $\sigma^2$에 대한 미분을 정리하면 다음과 같다.

$$ \dfrac{N}{2\sigma^2} - \dfrac{1}{2(\sigma^2)^2}\sum_{i=1}^N (x_i-\mu)^2 = 0 $$$$ \sigma^2 = \dfrac{1}{N}\sum_{i=1}^N (x_i-\mu)^2 = \dfrac{1}{N}\sum_{i=1}^N (x_i-\bar{x})^2 = s^2 $$
In:
np.random.seed(0)
mu0 = 1
sigma0 = 2
x = sp.stats.norm(mu0, sigma0).rvs(1000)
xbar = x.mean()
s2 = x.std(ddof=1)
xbar, s2
Out:
(0.90948658501960922, 1.9750540913890255)

다변수 정규 분포의 모수 추정

다변수 가우시안 정규 분포의 확률 밀도 함수는 다음과 같다.

$$ p(x ; \theta ) = (x ; \mu, \Sigma) = \dfrac{1}{(2\pi)^{D/2} |\Sigma|^{1/2}} \exp \left( -\dfrac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right) $$

그런데 $N$번의 반복 시행으로 표본 데이터가 $x_1, \cdots, x_N$ 가 있는 경우에는 모두 독립이므로 전체 확률 밀도 함수는 각각의 확류 밀도 함수의 곱과 같다.

$$ L(\theta;x_1, \cdots, x_N) = p(x_1, \cdots, x_N;\theta) = \prod_{i=1}^N \dfrac{1}{(2\pi)^{D/2} |\Sigma|^{1/2}} \exp \left( -\dfrac{1}{2} (x_i-\mu)^T \Sigma^{-1} (x_i-\mu) \right)$$

미분을 쉽게 하기 위해 로그 변환을 한 Log-Likelihood 를 구하면 다음과 같다. 여기에서 상수 부분은 모아서 $C$로 표기하였다.

$$ \begin{eqnarray*} \log L &=& \log P(x_1, \cdots, x_N;\theta) \\ &=& \sum_{i=1}^N \left\{ -\log((2\pi)^{D/2} |\Sigma|^{1/2}) - \dfrac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right\} \\ &=& C -\dfrac{N}{2} \log|\Sigma| - \dfrac{1}{2} \sum_i^N (x_i-\mu)^T \Sigma^{-1} (x_i-\mu) \end{eqnarray*} $$

여기에서 기호를 단순하게 하기 위해 precision matrix $\Sigma ^{-1} $를 $ \Lambda$ 로 표시하자.

$$\Lambda = \Sigma^{-1}$$

$$ \begin{eqnarray*} \log L &=& C + \dfrac{N}{2} \log|\Lambda| - \dfrac{1}{2} \sum_i^N(x_i-\mu)^T \Lambda (x_i-\mu) \end{eqnarray*} $$

이 확률 밀도 함수가 최대가 되는 모수 값을 찾기 위해서는 각각의 모수로 미분한 값이 0이 되어야 한다.

$$ \begin{eqnarray} \dfrac{\partial \log L}{\partial \mu} &=& - \dfrac{\partial}{\partial \mu} \sum_{i=1}^N (x_i-\mu)^T \Lambda (x_i-\mu) \\ &=& \sum_{i=1}^N 2\Lambda (x_i - \mu) \\ &=& -2\Lambda \sum_{i=1}^N (x_i - \mu) \\ &=& 0 \end{eqnarray} $$$$ \begin{eqnarray} \dfrac{\partial \log L}{\partial \Lambda} &=& \dfrac{\partial}{\partial \Lambda} \dfrac{N}{2} \log|\Lambda| - \dfrac{\partial}{\partial \Lambda} \dfrac{1}{2} \sum_{i=1}^N (x_i-\mu)^T\Lambda (x_i-\mu)\\ &=& \dfrac{\partial}{\partial \Lambda} \dfrac{N}{2} \log|\Lambda| - \dfrac{\partial}{\partial \Lambda} \dfrac{1}{2} \sum_{i=1}^N \text{tr}( (x_i-\mu)^T\Lambda (x_i-\mu)) \\ &=& \dfrac{\partial}{\partial \Lambda} \dfrac{N}{2} \log|\Lambda| - \dfrac{\partial}{\partial \Lambda} \dfrac{1}{2} \sum_{i=1}^N \text{tr}( (x_i-\mu)(x_i-\mu)^T\Lambda) \\ &=& 0 \end{eqnarray} $$

이 두 식을 풀면 주어진 데이터 표본에 대해 모수의 likelihood를 가장 크게 하는 모수의 값을 구할 수 있다. 첫번째 식을 풀면 다음과 같다.

$$ \sum_{i=1}^N (x_i - \mu) = 0 $$$$ \mu = \dfrac{1}{N}\sum_{i=1}^N x_i = \bar{x} $$

두번째 식을 풀면 다음과 같다.

$$ \dfrac{N}{2} \Lambda^{-T} = \dfrac{1}{2}\sum_{i=1}^N (x_i-\mu)(x_i-\mu)^T $$

$$ \Sigma = \dfrac{1}{N}\sum_{i=1}^N (x_i-\bar{x})(x_i-\bar{x})^T $$

In:
np.random.seed(0)
mu0 = np.array([0, 1])
sigma0 = np.array([[1, 0.2], [0.2, 4]])
x = sp.stats.multivariate_normal(mu0, sigma0).rvs(1000)
xbar = x.mean(axis=0)
S2 = np.cov(x, rowvar=0)
print(xbar)
print(S2)
[-0.0126996   0.95720206]
[[ 0.96100921  0.16283508]
 [ 0.16283508  3.80507694]]

질문/덧글

아직 질문이나 덧글이 없습니다. 첫번째 글을 남겨주세요!