[통계] 베이지안 회귀분석

 회귀는 데이터로부터 모델을 추정하는 한 방법이다. 최소자승법이 잔차를 최소화 시키는 방법이라면, 베이지안 회귀는 가능도 최대화가 목적이다. 

 이 글의 최종 목표는 베이지안 회귀의 원리를 이해하고, python package인 pymc3을 활용까지 다뤄보는 것이다.

 이번 포스팅에서는 우선 원리를 이해하는 것을 목표로 한다.

---

모델

 가장 기본적인 선형 모델은 다음과 같이 쓸 수 있다.

 $y = \theta_{1}x + \theta_{2}$ + $\epsilon$

 이 때 $\theta_{1}$, $\theta_{2}$는 모델의 파라미터, $x$, $y$는 관측 값, $\epsilon$은 오차를 나타낸다. 

 $x$, $y$라는 관측값으로 부터 $\theta_{1}$, $\theta_{2}$를 추정해 나가는 것이 목표다.

---

가능도

 확률이 모델에서 관측값이 나올 확률이라면, 가능도는 관측값에서 모델이 지닌 가능도를 평가 한다. 이 모델이 그럴싸 한지 보기 때문에 영어로 likelihood라고 부른다고 이해(암기?)하고 있다.

 선형 모델에서 오차 $\epsilon$은 평균 0, 표준편차 $\sigma$의 정규분포를 가정한다. 이를 식으로 표현하면 다음과 같다

$\epsilon \sim N(0, \sigma^{2})$

 그리고 이를 통해 $y$의 확률 분포를 나타내면 다음과 같다. 

$y \sim N(\theta_{1}x + \theta_{2}, \sigma^{2})$

 가능도는 정규분포의 확률밀도 함수를 통해 구할 수 있다. 

 평균 $\mu$, 표준편차 $\sigma$인 정규분포를 생각해보자. 

$x \sim N(\mu, \sigma^{2})$

 이 정규분포의 확률밀도는 다음과 같다.

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

 이 때, 관측된 값 $x_{i}$의 가능도는 $x_{i}$가 확률밀도에서 갖는 값을 나타낸다. 그리고 가능도는 모든 관측된 값의 가능도의 곱으로 나타낸다. 우리의 관측값은 모델로부터 독립적으로 추출된 값이므로, 각 관측값에 대한 확률밀도의 곱이 가능도가 된다.

 이를 이용해, 위의 $y$식의 가능도, 즉, 관측값 $x$, $y$에서 현재 모델의 모수의 가능도를 구하면 다음과 같다.

$L(\theta_{1}, \theta_{2}; x, y) = \prod_{i=1}^N \frac{1}{{\sigma}\sqrt{2\pi}}\exp \left( -\frac{(y_{i} - (\theta_{1}x_{i} + \theta_{2}))^{2}}{2{\sigma}^{2}} \right)$

---

베이즈 추론

 베이즈 추론은 관측값을 통해, 모델의 확률분포를 업데이트 하는 추론이다. 자세한 것은 여기  

 관측값 $x$, $y$를 줄여서 $\mathbf{x}$, 파라미터 $\theta_{1}$, $\theta_{2}$를 줄여서 $\boldsymbol{\theta}$로 두자. 그럼 베이즈 식은 다음과 같다.

$P(\boldsymbol{\theta}|\mathbf{x}) = \frac{P(\mathbf{x}|\boldsymbol{\theta})P(\boldsymbol{\theta})}{P(\mathbf{x})}$

 

 목적은 $\mathbf{x}$를 통해 $\boldsymbol{\theta}$의 확률분포를 구하는 것이다.

 여기서 $P(\boldsymbol{\theta})$를 사전분포, $P(\boldsymbol{\theta}|\mathbf{x})$를 사후분포라고 한다. 이를 이용하면, 적당한 $\boldsymbol{\theta}$의 범위에 대해서 관측값에 대한 사후분포를 구할 수 있다. 하지만, $\boldsymbol{\theta}$의 개수가 늘어난다면, 우리가 사후분포를 탐색해야 하는 공간은 지수 함수 꼴로 늘어나기 때문에 파라미터 공간을 탐색을 하기 위한 전략을 짜야 한다.

---

MCMC 샘플링

 여기서 등장하는 것이 마르코프 체인 몬테카를로법(MCMC)다. 

 마르코프 체인은 어떤 상태가 바로 전 상태에 의존하는 경우를 말한다. 몬테카를로는 무작위 샘플링을 하는 것이니 둘을 합친 MCMC는 어떤 상태가 끊임없이 움직이고 있는 체인 상태의 샘플링을 말한다.

 MCMC는 특정 조건에서 마르코프 체인이 정상 상태 분포(steady-state distribution)로 수속한다는 특징(ergodic)을 이용한다. 마르코프 연쇄가 ergodic하기 위해서는, 전이 상태가 하나의 값으로 정착 되지 않을 조건, 그리고 전이가 주기성을 띄지 않을 조건을 충족해야 한다. 이 조건을 만족하는 전이 알고리즘을 이용하면 어떤 파라미터 값에서 탐색을 시작하더라도 탐색 과정에서 발견된 값들이 이루는 분포가 정상 상태 분포를 이룰 것이고, 이 분포를 통해 모델 파라미터의 확률분포를 추정할 수 있다는 것이 MCMC 베이지안 추론이다. 

 MCMC의 정상 상태 분포는 파라미터의 확률분포를 타겟으로 한다. 이를 위해서는 파라미터 값을 탐색할 때 마르코프 체인의 상태 전이 확률이 더 그럴싸한 파라미터 값의 상태로 전이를 해줘야 한다. 즉, 샘플링이 더 그럴싸한 파라미터를 중심으로 이루어져야, 샘플링이 이루는 분포가 우리가 타겟으로 하는 파라미터 주변으로 이루어질 것이다.

 이를 위해 Metropolis-Hastings(MH) 알고리즘이 고안 되었다. MH 알고리즘에서는 사전분포를 제안하고, 현재 샘플링 파라미터 $\boldsymbol{\theta_{t}}$가 다음 샘플 $\boldsymbol{\theta_{new}}$로 전이될 확률을 사후확률의 비로 결정한다. 

 우선 사후확률의 비를 식으로 표현하면 다음과 같다.

$r = \frac{p(\boldsymbol{\theta_{new}} | \mathbf{x})}{p(\boldsymbol{\theta_{t}} | \mathbf{x})} = \frac{\frac{p(\mathbf{x}|\boldsymbol{\theta_{new}}) p(\boldsymbol{\theta_{new}})}  {p(\mathbf{x})}}{\frac{p(\mathbf{x}|\boldsymbol{\theta_{t}}) p(\boldsymbol{\theta_{t}})}  {p(\mathbf{x})}} = \frac{p(\mathbf{x}|\boldsymbol{\theta_{new}}) p(\boldsymbol{\theta_{new}})}{p(\mathbf{x}|\boldsymbol{\theta_{t}}) p(\boldsymbol{\theta_{t}})}$

 그리고 전이 여부를 다음 식으로 결정한다. 

$\alpha(\boldsymbol{\theta_{t}}, \boldsymbol{\theta_{new}}) = \min(1, r)$

 윗 식들을 보면 $p(\mathbf{x})$가 사라지고, 가능도 $p(\mathbf{x} | \boldsymbol{\theta})$와 사전확률분포 $p(\boldsymbol{\theta})$만으로 파라미터 탐색을 하고 있다는 것을 알 수 있다. 사전 확률 분포는 MH알고리즘에서는 proposal 분포라고도 하는데, 사전에 파라미터에 대한 정보가 없는 경우, 이 분포는 정보를 포함하지 않는 uniform 분포나 normal 분포를 이용한다. 위 파트에서 가능도를 구했고, 분포에 대한 정보는 설계자가 지정하므로, 우리는 r을 구할 수 있다!

 그럼, 전이 확률을 살펴보자. 만약, $\boldsymbol{\theta_{new}}$에서의 사후확률이 $\boldsymbol{\theta_{t}}$ 보다 크다면 r은 1보다 크게 되어 다음 상태 $\boldsymbol{\theta_{t+1}}$은 $\boldsymbol{\theta_{new}}$로 전이된다. 만약 새로 샘플링 된 파라미터가 기존 파라미터보다 작은 사후확률을 갖는다면 r이 작아지며, $\boldsymbol{\theta_{new}}$로 전이될 확률이 낮아지고, $\boldsymbol{\theta_{t}}$는 그대로 유지되어 $\boldsymbol{\theta_{t+1}}$이 되는 경우가 발생한다.

 이런 식으로 파라미터 공간을 탐색하면 우리가 찾고자 하는 파라미터 부근 공간에서 효율적인 탐색을 할 수 있게 된다.

 베이즈 추론에서 MH 알고리즘이 작동하는 원리를 간략하자면,  proposal 분포에서 추출한 파라미터를 그럴싸한 공간을 위주로 탐색하고, 그 탐색 경로가 이루는 분포를 통해 파라미터의 확률분포를 추정하는 알고리즘이다. 

---

 다음 글에서는 python 패키지인 pymc3을 활용한 베이지안 회귀분석을 정리 해봐야겠다.

참고 사이트

MH 알고리즘 :  http://www.secmem.org/blog/2019/01/11/mcmc/

선형 회귀모델 : https://coffeewhale.com/bayesian/linear/regression/2019/10/19/bayesian-lr/

(일본어)

http://bayesmax.sblo.jp/article/187238305.html

https://www.slideshare.net/yoshitaket/32-35647139

https://stats.biopapyrus.jp/bayesian-statistics/bayesian-inference-mcmc.html

http://web.sfc.keio.ac.jp/~maunz/BS14/BS14-10_GC.pdf