수식이 깨질 경우 새로고침을 눌러주세요.
Estimation (추정) ③¶
베이지안 모수 추정(Bayesian parameter estimation)¶
모수의 값에 해당하는 특정한 하나의 숫자를 계산하는 것이 아니라 모수의 값이 가질 수 있는 모든 가능성, 즉 모수의 분포를 계산하는 작업이다.
이때 계산된 모수의 분포를 표현 방법은 두 가지가 있다.
비모수적(non-parametric) 방법
- 샘플을 제시한 후 히스토그램와 같은 방법으로 임의의 분포를 표현한다.
- MCMC(Markov chain Monte Carlo) 몬테카를로 방법에서 사용한다.
- Likelihood 전체 모양을 구해야 하는 경우 (non-parametric) $$ p(\theta \;|\; x_{1:N}) $$
모수적(parametric) 방법
- 모수의 분포를 잘 알려진 확률 분포 모형을 사용하여 나타낸다.
- 이렇게 하면 모수를 나타내는 확률 분포 수식이 다시 모수(parameter)를 가지게 되는데 이를 hyper-parameter라고도 부른다.
- 모수적 방법은 결국 hypter-parameter의 값을 숫자로 계산하는 작업이 된다.
베이지안 모수 추정의 기본 원리¶
$$ P(\theta) \;\;←^{update}\;\; P(\theta \mid x_{1},\ldots,x_{N}) $$베이지안 모수 추정 방법은 다음 공식을 사용하여 모수의 분포를 갱신(update)하는 작업이다.
$$
p(\theta \mid x_{1},\ldots,x_{N}) = \dfrac{p(x_{1},\ldots,x_{N} \mid \theta) \cdot p(\theta)}{p(x_{1},\ldots,x_{N})} \\
\downarrow
\\
p(\theta \mid x_{1},\ldots,x_{N})\propto p(x_{1},\ldots,x_{N} \mid \theta ) \cdot p(\theta)
$$
사전(Prior) 분포 $$P(\theta)$$
- 사전 분포는 베이지안 추정 작업을 하기 전에 이미 알고 있던 모수의 분포를 뜻한다.
- 아무런 지식이 없는 경우에는 보통 uniform 분포(균일분포)나 표준정규 분포를 사용한다. $$\text{Beta}(1,1) \;\;\; Or \;\;\;\mathcal{N}(0, 1) $$
$$ P(x_{1},\ldots,x_{N} \mid \theta) $$$$ \theta를\; 알고\; 있는\; 상태에서의\; 데이터(x_{1},\ldots,x_{N})가\; 나올\; 조건부\; 확률\; 분포 $$ $$ p(x \;|\; \theta)는\; 문제에서\; 주어지거나\; 특정\; 모형으로\; 가정한다. $$Likelihood 분포
사후(Posterior) 분포
- 구하고자 하는 답
- 샘플에 의해 변화된 모수(theta)의 분포 $$P(\theta \mid x_{1},\ldots,x_{N})$$
$$
데이터\; x_{1},\ldots,x_{N}을\; 알고있는\; 상태에서의\; \theta에\; 대한\; 조건부\; 확률\; 분포이다.\\
\\ 우리가\; 베이지안\; 모수\; 추정을\; 통해 \;구하고자 \;하는 \;것이\; 바로 \;이 \;Posterior다.
$$
MAP (Maximum A Posterior)¶
베이지안 모수 추정의 방법이다.
랜덤 변수 X의 파라미터에 대해 점 추정을 하는 것이 아니고, 분포(distribution) 정보를 제공한다.
모수적 방법(Hyper-parameter)
- theta에 대한 distribution이 몇 개의 parameter로 결정되는 경우 그 parameter를 hyper-parameter라고 한다.
- Likelihood를 최대화하는 하이퍼-파라미터 alpha 값을 구한다. $$ \arg \max_{\alpha} p(\theta \;|\; x_{1:N}, \alpha) $$
이번 포스팅에서는 베이지안 모수적 방법에 대해 알아 보자.¶
Beta 분포를 활용한 베이지안 모수 추정 ①¶
Beta 분포는 값이 바운더리된 것을 추정할 때 사용한다.
- Bernoulli 분포의 parameter는 0 이상 1 이하의 값을 가지는 실수.
- theta라는 1개의 모수 추정 $$ 0\leq \theta \leq 1 \\ \downarrow \\ \theta는\; Beta\; 분포로\; 나타낼\; 수\; 있다. $$
- Beta 분포는 파라미터 a, b를 가진다.
- 따라서 Bernoulli 랜덤 변수의 parameter theta는 hyper-parameter a, b를 가진다.
$$ Posterior = \frac{Likelihood \cdot Prior}{Normalizing} $$$$ 0\leq \theta \leq 1 \rightarrow \;\; ∴\;\;\; \theta \sim Beta(a=1, b=1) $$먼저, 수리통계학으로 풀어보자.
- 가장 단순한 이산 확률 분포인 베르누이 분포의 모수를 베이지안 추정법으로 추정해보자.
$$ \begin{eqnarray*} P(\theta) &=& \frac{1}{B(a,b)}\cdot \theta^{a-1} (1-\theta)^{b-1} \\ &=& \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\cdot \theta^{a-1} (1-\theta)^{b-1} \\ &=& \frac{\Gamma(2)}{\Gamma(1)\Gamma(1)}\cdot \theta^{a-1} (1-\theta)^{b-1} \\ &=& \frac{1!}{0!\cdot0!}\cdot \theta^{a-1} (1-\theta)^{b-1} \\ &=& \theta^{a-1} (1-\theta)^{b-1} \end{eqnarray*} $$(1) Prior
- 우리가 갖고 있는 데이터는 특정 분포를 따를 것으로 가정했다.
여기서는 베르누이 분포
- 이 특정 분포의 파라미터가 존재하는 구간(interval)을 찾자.
베르누이는 0이상 1이하
- 사전 분포는 파라미터가 존재하는 구간을 정의역으로 갖는 함수(분포)이다.
베르누이는 Beta(1,1)
$$
∴\;\;\;prior \;=\; P(\theta) = \bbox[7px, border:2px solid red]{\theta^{a−1}(1−\theta)^{b−1} \;\;\; (a=1, b=1)}
$$
(2) Likelihood
- 데이터는 서로 모두 독립이며, 베르누이 분포를 따르므로 다음과 같이 likelihood를 구하면 된다. $$ P(x_{1},\ldots,x_{N} \mid \theta) = \prod_{i=1}^N \theta^{x_i} (1 - \theta)^{1-x_i} $$
(3) Posterior
- Bayes Rule을 사용하여 사후 분포를 구하면 다음과 같이 갱신된 베타 분포가 된다.
$$
\begin{eqnarray}
P(\theta \mid x_{1},\ldots,x_{N})
&\propto & P(x_{1},\ldots,x_{N} \mid \theta) P(\theta) \\
&=& \prod_{i=1}^N \theta^{x_i} (1 - \theta)^{1-x_i} \cdot \theta^{a−1}(1−\theta)^{b−1} \\
&=& \theta^{\sum_{i=1}^N x_i + a−1} \cdot (1 - \theta)^{\sum_{i=1}^N (1-x_i) + b−1 }, \;\;\;\; (단,\;\sum_{i=1}^N x_i = N_1,\;\; \sum_{i=1}^N (1-x_i) = N_0) \\
&=& \theta^{N_1 + a−1} (1 - \theta)^{N_0 + b−1 } \\
&=& \theta^{a'−1} (1 - \theta)^{b'−1 } \\
\end{eqnarray}
$$
$$ a' = N_1 + a $$$$ b' = N_0 + b $$
- 이렇게 사전 분포와 사후 분포가 같은 확률 분포 모형을 가지게 하는 사전 분포를 conjugate prior 라고 한다.
- 갱신된 하이퍼 모수의 값은 다음과 같다.
아래 코드를 보면서 prior와 posterior를 갱신해보자.¶
In [1]:
# -*- coding: utf-8 -*-
import seaborn as sns
import pandas as pd
import scipy as sp
import matplotlib as mpl
import matplotlib.pylab as plt
import numpy as np
%matplotlib inline
In [2]:
theta0 = 0.6
a0, b0 = 1, 1
print("step 0: mode = unknown")
plt.figure(figsize=(10,6))
xx = np.linspace(0, 1, 1000)
plt.plot(xx, sp.stats.beta(a0, b0).pdf(xx), label="initial");
np.random.seed(0)
x = sp.stats.bernoulli(theta0).rvs(50)
N0, N1 = np.bincount(x, minlength=2)
a1, b1 = a0 + N1, b0 + N0
plt.plot(xx, sp.stats.beta(a1, b1).pdf(xx), label="1st");
print("step 1: mode =", (a1 - 1)/(a1 + b1 - 2))
x = sp.stats.bernoulli(theta0).rvs(50)
N0, N1 = np.bincount(x, minlength=2)
a2, b2 = a1 + N1, b1 + N0
plt.plot(xx, sp.stats.beta(a2, b2).pdf(xx), label="2nd");
print("step 2: mode =", (a2 - 1)/(a2 + b2 - 2))
x = sp.stats.bernoulli(theta0).rvs(50)
N0, N1 = np.bincount(x, minlength=2)
a3, b3 = a2 + N1, b2 + N0
plt.plot(xx, sp.stats.beta(a3, b3).pdf(xx), label="3rd");
print("step 3: mode =", (a3 - 1)/(a3 + b3 - 2))
x = sp.stats.bernoulli(theta0).rvs(50)
N0, N1 = np.bincount(x, minlength=2)
a4, b4 = a3 + N1, b3 + N0
plt.plot(xx, sp.stats.beta(a4, b4).pdf(xx), label="4th");
print("step 4: mode =", (a4 - 1)/(a4 + b4 - 2))
plt.legend()
plt.show()
update가 될수록 베르누이의 파라미터(theta)가 위치하는 분포가 정확해진다는 것을 확인할 수 있다.¶
4th 분포 안에 theta는 존재한다.
conjugate prior
- Prior → Posterior, 사후확률이 업데이트되고, 그 사후확률이 사전확률이 되는 과정이 반복
- Data가 많을 수록 분포를 정확하게 찾아간다.
켤레 사전 분포(Conjugate prior distribution)¶
사전확률과 사후확률이 동일한 분포를 따른다면 계산이 매우 편해지기 때문에 베이즈 통계학에서 많이 사용한다.
| 우도 | 켤레사전분포(Prior) | 사후확률분포(Posterior) | patameter |
|---|---|---|---|
| 이항분포 | 베타분포 Beta(α,β) | 베타분포 Beta(α′,β′) | 성공횟수 : α−1, 실패횟수 : β−1 |
| 다항분포 | 디리클레분포 Dir(α) | 디리클레분포 Dir(α′) | i번째 범주가 나타날 횟수 : $$α_i−1$$ |
Dirichlet 분포를 활용한 베이지안 모수 추정 ②¶
$$ 0\leq \theta_i \leq 1, \;\;\;\; (단,\;\sum_{i=1}^K \theta_i = 1) $$클래스 개수가 K인 카테고리 분포의 모수 벡터를 베이지안 추정법으로 추정해보자.
- 베르누이와 다른 점은 파라미터 벡터를 추정하는 것이다.(모수가 2개 이상)
$$ P(\theta) \propto \prod_{i=1}^K \theta_i^{\alpha_i - 1} \;\;\; (\alpha_i = 1, \; \text{ for all } i) $$(1) Prior
$$ \begin{eqnarray*} P(x_{1},\ldots,x_{N} \mid \theta) &=& P(x_{1} \;|\; \theta) \cdots P(x_{N} \;|\; \theta) \\ &=& \prod_{j=1}^N P(x_{j} \;|\; \theta) \\ &=& \prod_{j=1}^N \prod_{i=1}^K \theta_i^{x_{j,i}} \end{eqnarray*} $$(2) Likelihood 데이터는 모두 독립이며, 카테고리 분포를 따르므로
$$ \begin{eqnarray} P(\theta \mid x_{1},\ldots,x_{N}) &\propto & P(x_{1},\ldots,x_{N} \mid \theta) P(\theta) \\ &=& \prod_{j=1}^N \prod_{i=1}^K \theta_k^{x_{j,i}} \cdot \prod_{i=1}^K \theta_i^{\alpha_i - 1} \\ &=& \prod_{i=1}^K \theta^{\sum_{j=1}^N x_{j,i} + \alpha_i − 1} \;\;\;\; (\sum_{j=1}^N x_{j,i} = N_i,\;\; N_i는\; x_i가\; 나온\; 횟수) \\ &=& \prod_{i=1}^K \theta^{N_i + \alpha_i −1} \\ &=& \prod_{i=1}^K \theta^{\alpha'_i −1} \\ \end{eqnarray} $$(3) Posterior
- Bayes Rule을 사용하여 사후 분포를 구하면 다음과 같다.
- 이 경우에도 conjugate prior 임을 알 수 있다.
- 갱신된 Hyper Parameter는 다음과 같다. $$ \alpha'_k = N_k + \alpha_k $$
다음을 보자.¶
In [3]:
def plot_dirichlet(alpha):
def project(x):
n1 = np.array([1, 0, 0])
n2 = np.array([0, 1, 0])
n3 = np.array([0, 0, 1])
n12 = (n1 + n2)/2
m1 = np.array([1, -1, 0])
m2 = n3 - n12
m1 = m1/np.linalg.norm(m1)
m2 = m2/np.linalg.norm(m2)
return np.dstack([(x-n12).dot(m1), (x-n12).dot(m2)])[0]
def project_reverse(x):
n1 = np.array([1, 0, 0])
n2 = np.array([0, 1, 0])
n3 = np.array([0, 0, 1])
n12 = (n1 + n2)/2
m1 = np.array([1, -1, 0])
m2 = n3 - n12
m1 = m1/np.linalg.norm(m1)
m2 = m2/np.linalg.norm(m2)
return x[:,0][:, np.newaxis] * m1 + x[:,1][:, np.newaxis] * m2 + n12
# 입실론
eps = np.finfo(float).eps * 10
X = project([[1-eps,0,0], [0,1-eps,0], [0,0,1-eps]])
import matplotlib.tri as mtri
triang = mtri.Triangulation(X[:,0], X[:,1], [[0, 1, 2]])
refiner = mtri.UniformTriRefiner(triang)
triang2 = refiner.refine_triangulation(subdiv=6)
XYZ = project_reverse(np.dstack([triang2.x, triang2.y, 1-triang2.x-triang2.y])[0])
pdf = sp.stats.dirichlet(alpha).pdf(XYZ.T)
plt.tricontourf(triang2, pdf)
plt.axis("equal")
plt.show()
In [4]:
np.dstack([[1,0,0],[0,1,0],[0,0,1]])[0]
Out[4]:
In [5]:
theta0 = np.array([0.2, 0.6, 0.2])
np.random.seed(0)
# 0,1,2 중 theta의 확률로 20번 뽑아라.
x1 = np.random.choice(3, 20, p=theta0)
N1 = np.bincount(x1, minlength=3)
x2 = np.random.choice(3, 100, p=theta0)
N2 = np.bincount(x2, minlength=3)
x3 = np.random.choice(3, 1000, p=theta0)
N3 = np.bincount(x3, minlength=3)
In [6]:
a0 = np.ones(3) / 3
plot_dirichlet(a0)
In [7]:
a1 = a0 + N1
plot_dirichlet(a1)
print((a1 - 1)/(a1.sum() - 3))
In [8]:
a2 = a1 + N2
plot_dirichlet(a2)
print((a2 - 1)/(a2.sum() - 3))
In [9]:
a3 = a2 + N3
plot_dirichlet(a3)
print((a3 - 1)/(a3.sum() - 3))
사후분포의 파라미터가 업데이트됨에 따라 다항분포 파라미터의 구간 추정이 정확해지는 것을 알 수 있다.¶
Gaussian 정규 분포를 활용한 베이지안 모수 추정 ③¶
$$ -\infty < \mu < \infty $$정규 분포의 기대값 모수(mu)를 베이지안 방법으로 추정하자.
- 분산은 알고 있다고 가정하자.
$$ P(\mu) \sim N(\mu_0, \sigma^2_0) = \dfrac{1}{\sqrt{2\pi\sigma_0^2}} \exp \left(-\dfrac{(\mu-\mu_0)^2}{2\sigma_0^2}\right)$$(1) Prior
(2) Likelihood
- 데이터는 서로 독립이므로, $$ P(x_{1},\ldots,x_{N} \mid \mu) = \prod_{i=1}^N N(x_i \mid \mu ) = \prod_{i=1}^N \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp \left(-\dfrac{(x_i-\mu)^2}{2\sigma^2}\right) $$
(3) Posterior $$ \begin{eqnarray} P(\mu \mid x_{1},\ldots,x_{N}) &\propto & P(x_{1},\ldots,x_{N} \mid \mu) P(\mu) \\ &\propto & \prod_{i=1}^N \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp \left(-\dfrac{(x_i-\mu)^2}{2\sigma^2}\right)\cdot\dfrac{1}{\sqrt{2\pi\sigma_0^2}} \exp \left(-\dfrac{(\mu-\mu_0)^2}{2\sigma_0^2}\right) \\ \\&\propto & \exp \left(-\dfrac{(\mu-\mu'_0)^2}{2\sigma_0^{'2}}\right) \\ \end{eqnarray} $$
$$ \begin{eqnarray} \mu'_0 &=& \dfrac{\sigma^2}{N\sigma_0^2 + \sigma^2}\mu_0 + \dfrac{N\sigma_0^2}{N\sigma_0^2 + \sigma^2} \dfrac{\sum x_i}{N} \\ &=& A\cdot\mu_0 + B\cdot\bar{X} \\ \\ \dfrac{1}{\sigma_0^{'2}} &=& \dfrac{1}{\sigma_0^{2}} + \dfrac{N}{\sigma^{'2}} \end{eqnarray} $$
- 베이지안 규칙을 사용하여 사후 분포를 구하면 다음과 같이 갱신된 하이퍼 모수를 가지는 정규 분포가 된다.
분산의 역수가 증가하므로 파라미터가 갱신될수록 분산은 감소한다. $$ \dfrac{N}{\sigma^{'2}} > 0 $$
In [10]:
mu, sigma2 = 2, 4
mu0, sigma20 = 0, 1
xx = np.linspace(1, 3, 1000)
np.random.seed(0)
N = 10
x = sp.stats.norm(mu).rvs(N)
mu0 = sigma2/(N*sigma20 + sigma2) * mu0 + (N*sigma20)/(N*sigma20 + sigma2)*x.mean()
sigma20 = 1/(1/sigma20 + N/sigma2)
plt.plot(xx, sp.stats.norm(mu0, sigma20).pdf(xx), label="1st");
print(mu0)
N = 20
x = sp.stats.norm(mu).rvs(N)
mu0 = sigma2/(N*sigma20 + sigma2) * mu0 + (N*sigma20)/(N*sigma20 + sigma2)*x.mean()
sigma20 = 1/(1/sigma20 + N/sigma2)
plt.plot(xx, sp.stats.norm(mu0, sigma20).pdf(xx), label="2nd");
print(mu0)
N = 50
x = sp.stats.norm(mu).rvs(N)
mu0 = sigma2/(N*sigma20 + sigma2) * mu0 + (N*sigma20)/(N*sigma20 + sigma2)*x.mean()
sigma20 = 1/(1/sigma20 + N/sigma2)
plt.plot(xx, sp.stats.norm(mu0, sigma20).pdf(xx), label="3rd");
print(mu0)
N = 100
x = sp.stats.norm(mu).rvs(N)
mu0 = sigma2/(N*sigma20 + sigma2) * mu0 + (N*sigma20)/(N*sigma20 + sigma2)*x.mean()
sigma20 = 1/(1/sigma20 + N/sigma2)
plt.plot(xx, sp.stats.norm(mu0, sigma20).pdf(xx), label="4th");
print(mu0)
plt.axis([1, 3, 0, 20])
plt.legend()
plt.show()
Inverse Gamma 분포를 통해 가우시안 정규 분포의 분산 추정¶
Gaussian Normal 분포의 parameter $$ \sigma^2 > 0$$ $$ {(\sigma^2})^{-1}는\; Gamma\; 분포\; \text{Gam}(a, b)로\; 나타낼\; 수\; 있다. $$
Posterior Predictive¶
$$ p(x_{N+1}\;|\;x_{1:N}) = \int p(x\;|\;\theta)\cdot p(\theta\;|\;x_{1:N}) d\theta $$pf)
$$ With \;Law \;of \;total \;probability \\ P(x) = \int P(x\;|\;\theta)\cdot P(\theta)\;d\theta $$$$ \begin{eqnarray} ∴ \;\;\;\; p(x_{N+1}\;|\;x_{1:N}) &=& \int P(x_{N+1}\;|\;\theta,\;x_{1:N})\cdot P(\theta\;|\;x_{1:N})\;d\theta \\ &=& \int P(x_{N+1}\;|\;\theta)\cdot P(\theta\;|\;x_{1:N})\;d\theta \;\;\; (∵\; x_{N+1}과\; x_{1:N}은\; 서로\; 독립이므로)\\ \end{eqnarray} $$MLE
- plugin approximation $$ p(x_{N+1}\;|\;x_{1:N}) = p(x\;|\;\theta^{\ast}) $$
$$ p(\theta\;|\;\alpha^{\ast}) \;\;→\;\; 대부분\; 시뮬레이션으로\; 값을\; 얻어낸다. $$MAP $$ p(x_{N+1}\;|\;x_{1:N}) = \int p(x\;|\;\theta) \cdot p(\theta\;|\;\alpha^{\ast}) d\theta $$
MLE와 MAP의 차이점¶
MLE
- Data만으로 파라미터를 추정한다.
- 편향된 데이터로 MLE를 하게 되면 올바르지 않은 값을 구하게 된다.
- 극단적인 잘못된 데이터 수집을 했을 경우 정확성이 낮아진다.
예를 들어, 동전을 100번 던졌는데, 앞면이 70번 나왔다.
MLE의 경우 앞면이 나올 이론적 확률 = 0.5를 전혀 반영하지 않기 때문에 잘못된 값을 추정할 수 있다.
또다른 예는, 강남대성 재수생들의 평균 모의고사 점수와 전국 꼴지 고등학교의 모의고사 점수 평균은 실제 모의고사 평균값과 차이가 있을 것이다.
MAP
- Prior라는 강력한 가설을 추가함. (data 이외의 강력한 가설)
- 사전분포를 통하면 MLE의 단점을 커버할 수 있다.