다운로드
작성자: admin 작성일시: 2016-05-30 12:05:26 조회수: 2970 다운로드: 219
카테고리: 기초 수학 태그목록:

베이지안 모수 추정

베이지안 모수 추정(Bayesian parameter estimation) 방법은 모수의 값에 해당하는 특정한 하나의 숫자를 계산하는 것이 아니라 모수의 값이 가질 수 있는 모든 가능성, 즉 모수의 분포를 계산하는 작업이다.

이때 계산된 모수의 분포를 표현 방법은 두 가지가 있다.

  1. 비모수적(non-parametric) 방법

    • 샘플을 제시한 후 히스토그램와 같은 방법으로 임의의 분포를 표현한다. MCMC(Markov chain Monte Carlo)와 같은 몬테카를로 방법에서 사용한다.
  2. 모수적(parametric) 방법

    • 모수의 분포를 잘 알려진 확률 분포 모형을 사용하여 나타낸다. 이렇게 하면 모수를 나타내는 확률 분포 수식이 다시 모수(parameter)를 가지게 되는데 이를 hyper-parameter라고도 부른다. 모수적 방법은 결국 hypter-parameter의 값을 숫자로 계산하는 작업이 된다.

여기에서는 모수적 방법의 몇 가지 간단한 예를 보인다.

베이지안 모수 추정의 기본 원리

베이지안 모수 추정 방법은 다음 공식을 사용하여 모수의 분포 $p(\theta)$를 $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})} \propto p(x_{1},\ldots,x_{N} \mid \theta ) \cdot p(\theta) $$

이 식에서

$p(\theta)$ 는 사전(Prior) 분포라고 한다. 사전 분포는 베이지안 추정 작업을 하기 전에 이미 알고 있던 모수 $\theta$의 분포를 뜻한다. 아무런 지식이 없는 경우에는 보통 uniform 분포 $\text{Beta}(1,1)$나 0 을 중심으로하는 정규 분포 $\mathcal{N}(0, 1)$를 사용한다

$p(\theta \mid x_{1},\ldots,x_{N})$ 는 사후(Posterior) 분포라고 한다. 수학적으로는 데이터 $x_{1},\ldots,x_{N}$가 알려진 상태에서의 $\theta$에 대한 조건부 확률 분포이다. 우리가 베이지안 모수 추정 작업을 통해 구하고자 하는 것이 바로 이 사후 분포이다.

$p(x_{1},\ldots,x_{N} \mid \theta)$ 분포는 Likelihood 분포라고 한다. 조건부 확률 분포 $p(\theta \mid x_{1},\ldots,x_{N})$ 에서 데이터 $x_{1},\ldots,x_{N}$을 이미 알고 있고 $\theta$가 미지수이지만 Likelihood에서는 이와 반대로 $\theta$를 이미 알고 있고 데이터 $x_{1},\ldots,x_{N}$ 가 나올 확률을 구한다.

베르누이 분포의 모수 추정

가장 단순한 이산 확률 분포인 베르누이 분포의 모수 $\theta$를 베이지안 추정법으로 추정해 본다.

베르누이 분포의 모수는 0부터 1사이의 값을 가지므로 사전 분포는 하이퍼 모수 $a=b=1$인 베타 분포로 한다.

$$ P(\theta) \propto \theta^{a−1}(1−\theta)^{b−1} \;\;\; (a=1, b=1)$$

데이터는 모두 독립적인 베르누이 분포의 곱이므로 Likelihood는 다음과 같이 이항 분포가 된다.

$$ P(x_{1},\ldots,x_{N} \mid \theta) = \prod_{i=1}^N \theta^{x_i} (1 - \theta)^{1-x_i} $$

베이지안 규칙을 사용하여 사후 분포를 구하면 다음과 같이 갱신된 하이퍼 모수 $a'$, $b'$를 가지는 베타 분포가 된다.

$$ \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} (1 - \theta)^{\sum_{i=1}^N (1-x_i) + b−1 } \\ &=& \theta^{N_1 + a−1} (1 - \theta)^{N_0 + b−1 } \\ &=& \theta^{a'−1} (1 - \theta)^{b'−1 } \\ \end{eqnarray} $$

이렇게 사전 분포와 사후 분포가 같은 확률 분포 모형을 가지게 하는 사전 분포를 conjugate prior 라고 한다.

갱신된 하이퍼 모수의 값은 다음과 같다.

$$ a' = N_1 + a $$$$ b' = N_0 + b $$
In:
theta0 = 0.6
a0, b0 = 1, 1
print("step 0: mode = unknown")

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()
step 0: mode = unknown
step 1: mode = 0.52
step 2: mode = 0.62
step 3: mode = 0.58
step 4: mode = 0.59

카테고리 분포의 모수 추정

다음으로 클래스 갯수가 $K$인 카테고리 분포의 모수 $\theta$ 벡터를 베이지안 추정법으로 추정해 본다.

카테고리 분포의 모수의 각 원소는 모두 0부터 1사이의 값을 가지므로 사전 분포는 하이퍼 모수 $\alpha_k=1$인 디리클리 분포로 한다.

$$ P(\theta) \propto \prod_{k=1}^K \theta_k^{\alpha_k - 1} \;\;\; (\alpha_k = 1, \; \text{ for all } k) $$

데이터는 모두 독립적인 카테고리 분포의 곱이므로 Likelihood는 다음과 같이 다항 분포가 된다.

$$ P(x_{1},\ldots,x_{N} \mid \theta) = \prod_{i=1}^N \prod_{k=1}^K \theta_k^{x_{i,k}} $$

베이지안 규칙을 사용하여 사후 분포를 구하면 다음과 같이 갱신된 하이퍼 모수 $\alpha'_i$를 가지는 디리클리 분포가 된다.

$$ \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 \prod_{k=1}^K \theta_k^{x_{i,k}} \cdot \prod_{k=1}^K \theta_k^{\alpha_k - 1} \\ &=& \prod_{k=1}^K \theta^{\sum_{i=1}^N x_{i,k} + \alpha_k − 1} \\ &=& \prod_{k=1}^K \theta^{N_k + \alpha_k −1} \\ &=& \prod_{k=1}^K \theta^{\alpha'_k −1} \\ \end{eqnarray} $$

이 경우에도 conjugate prior 임을 알 수 있다.

갱신된 하이퍼 모수의 값은 다음과 같다. $$ \alpha'_k = N_k + \alpha_k $$

In:
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:
theta0 = np.array([0.2, 0.6, 0.2])

np.random.seed(0)
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:
a0 = np.ones(3) / 3
plot_dirichlet(a0)
In:
a1 = a0 + N1
plot_dirichlet(a1)
print((a1 - 1)/(a1.sum() - 3))
[ 0.12962963  0.62962963  0.24074074]
In:
a2 = a1 + N2
plot_dirichlet(a2)
print((a2 - 1)/(a2.sum() - 3))
[ 0.20621469  0.62146893  0.17231638]
In:
a3 = a2 + N3
plot_dirichlet(a3)
print((a3 - 1)/(a3.sum() - 3))
[ 0.2033393   0.59421586  0.20244484]

정규 분포의 기댓값 모수 추정

이번에는 정규 분포의 기댓값 모수를 베이지안 방법으로 추정한다. 분산 모수 $\sigma^2$은 알고 있다고 가정한다.

기댓값은 $-\infty$부터 $\infty$까지의 모든 수가 가능하기 때문에 모수의 사전 분포로는 정규 분포를 사용한다.

$$ P(\mu) = 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)$$

데이터는 모두 독립적인 정규 분포의 곱이므로 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) $$$$ \begin{eqnarray} P(\theta \mid x_{1},\ldots,x_{N}) &\propto & P(x_{1},\ldots,x_{N} \mid \theta) P(\theta) \\ &\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} \\ \dfrac{1}{\sigma_0^{'2}} &=& \dfrac{1}{\sigma_0^{2}} + \dfrac{N}{\sigma^{'2}} \end{eqnarray} $$
In:
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()
1.95573083623
2.15546157111
1.87567401511
2.04777201998

질문/덧글

베이지언 추정 결과에 대해 tara*** 2016년 10월 1일 4:59 오후

모수적 방법으로 계산된 베이지언 추정치의 유의성에 관해 실무, 혹은 학계에서 자주 사용하는 기준이 있나요? 피크의 폭이 어느 이상 좁아야 유의한 결과로 사용가능하다던가, 혹은 피크 높이가 어느 이상 높아야 한다던가 하는, 여러 사람들에게 통용되는 기준이 존재하는지 궁금합니다.

카테고리 분포의 모수 추정에서 poio*** 2017년 2월 21일 8:44 오후

사전분포는 αk = 1/K 에서 알파가 아니고 세타가 맞지 않나요?

답변: 카테고리 분포의 모수 추정에서 관리자 2017년 2월 22일 3:10 오후

하이퍼 모수 $\alpha$를 설정하는 것은 맞습니다. 다만 $1/K$가 아니라 $1$이 되어야 합니다. 이 부분은 수정하였습니다.

반면이라고 하셨는데 똑같아서 혼동이 왔습니다. kekm*** 2018년 1월 5일 3:22 오후

p(x1,…,xN∣θ) 분포는 Likelihood 분포라고 한다. 현재 우리가 알고 있는 값은 데이터 x1,…,xN 이고 θ가 미지수이다. 이와 반대로 θ를 알고 있는 상태에서의 데이터 x1,…,xN 가 나올 조건부 확률 분포를 Likelihood 라고 한다.

답변: 반면이라고 하셨는데 똑같아서 혼동이 왔습니다. 관리자 2018년 1월 10일 10:01 오전

문구를 조금 바꾸어보았습니다. 여전히 어렵다면 다시 지적 부탁드립니다.