[통계] 베이지안 회귀분석 - 선형모델 with PyMC3

 python 패키지인 PyMC3을 이용, 지난 글에서 소개한 베이지안 회귀분석을 해보는 예제를 정리해본다.

 PyMC3은 베이지안 회귀를 이용한 다양한 모델링을 가능케 해주는 package다. 사용법도 쉽고, 무엇보다 홈페이지에 다양한 예제가 자세히 소개 되어 있다.

 Install 시, 주의사항은 python3.6 보다 높은 버전에서는 에러가 난다는 점. python3.6 버전에서 설치하는 것이 안정적이다. (홈페이지에서도 그렇게 언급하고 있다.) PyMC3는 theano를 기반으로, PyMC4는 tensorflow를 기반으로 동작한다고 한다.

---

예제 데이터 생성

 $y = 2x + 1 + \epsilon$

 위와 같은 모델을 통해 예제 데이터를 생성한다.

 이 때, 기울기 $\theta_{1}$는 2, 절편 $\theta_{2}$는 1, $\epsilon$에는 표준편차 $\sigma$ 0.5의 정규분포의 노이즈를 부여하고, 200개의 샘플을 작성한다.

코드)

import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm


""" 파라미터 """
size = 200
theta1 = 2
theta2 = 1
sigma = 0.5

""" 데이터 생성 """
x = np.linspace(start=0, stop=10, num=size)
y = theta1 * x + theta2
y = y + np.random.normal(loc=0, scale=sigma, size=size)

""" plot """
fig = plt.figure(figsize=(4, 4))
plt.scatter(x, y, s=2)
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()

 생성된 데이터는 다음과 같다.

---

회귀 모델

 pymc에서는 모델을 다음과 같이 생성한다.

코드)

""" model 생성 """
with pm.Model() as linear_model:

    # Prior probability (proposal distribution)
    a = pm.Normal('theta1', mu=0, sd=10)
    b = pm.Normal('theta2', mu=0, sd=10)
    eps = pm.HalfCauchy('sigma', 1)

    # Model
    y_est = a * x + b

    # Likelihood distribution
    likelihood = pm.Normal('y', mu=y_est, sd=eps, observed=y)

 먼저, with로 모델을 정의한다.

 그리고,  사전확률 분포 (prior probability)를 설계한다. 파라미터 값으로는 보통 정규분포나 균등분포와 같이 정보가 없는 분포를 주로 사용하고, 표준편차는 HalfCauchy라는 양의 값으로만 이루어진 heavy-tail의 분포를 주로 사용한다고 한다. (관련 정보는 여기를 참조)

 다음으로는 모델을 구성하고, 이로부터 Likelihood를 정의한다. 이 때 정의가 위의 사전확률 분포와 같은 pm.Normal이라는 랜덤변수 함수를 이용해 생성하는데, 차이점은 observed가 있다는 점이다. 문서에 따르면 랜덤변수 함수는 observed를 지정하면 likelihood distribution으로, 지정하지 않으면, 랜덤변수 함수로 정의된다고 한다. (여기의 2.Probability Distributions 항목 참조)

---

MCMC 샘플링

 다양한 MCMC샘플링이 가능하고 파라미터 조정이 가능하지만, pymc3에서는 NUTS라는 sampler를 권장하고 있다. NUTS sampler는 많은 부분이 자동화 되어 있고, 간편히 쓸 수 있다.

코드)

""" MCMC 샘플링 """
with linear_model:
    trace = pm.sample(3000, progressbar=True,  chains=2, cores=1)

 여기서는 with로 모델을 따로 열었는데, 사실 위에 이어서 pm.sample 함수를 동작 시켜도 된다.

 NUTS sampler를 이용, 3000번의 샘플링을 지정했다. progressbar는 True로 하면 현재 샘플링 진행과정을 보여준다.

 chains는 마르코프 체인의 개수를 나타내는데, 2개를 지정하면, 두 번의 샘플링을 진행한다. 한 개의 마르코프 체인으로는 현재 샘플링이 정상 분포 상태를 나타내는 것인지 확인할 방법이 없기 때문에, 2 이상의 chain을 설정하여 통계적으로 샘플링을 평가할 필요가 있다.

 cores는 현재 CPU core 개수로 병렬로 chain을 돌리는 기법에 사용되는데, 나는 1보다 큰 값에서 에러가 났고, 이를 디버깅할 욕구가 없었기에 1로 해두었다.

 이 외에도, tune이라는 단계를 조정 한다거나 하는 세부적인 인자들이 있다.

---

결과

가장 쉽게 뽑을 수 있는 결과는 사후분포 확인이다. 다음 코드로 바로 가시화 할 수 있다.

코드)

""" 사후분포 plot """
pm.traceplot(trace)

 theta1, theta2, sigma의 샘플링이 이루고 있는 사후분포를 왼쪽 3개의 그래프에서 확인할 수 있다. 

 우리가 설정한 값이 각각 2, 1, 0.5 였으므로, 대략적으로 올바른 값 주변으로 MCMC 샘플링이 이루어진 것을 알 수 있다.

 오른쪽 3개의 그래프는 각 마르코프 체인 샘플링 과정을 보여준다. 여기서 나온 sample value를 histgram으로 나타낸 것이 왼쪽 그래프다.

 오른쪽 그래프를 보고 샘플링이 잘 이뤄 졌는가를 평가해 볼 수 있다. 특정 sample value에서 스텝이 지나도 상태전이가 일어나지 않는 구간이 발생 한다거나 하는 경우에는 제대로 된 정상 상태에 도달하지 않은 경우가 있으므로, sampling을 잘 할 수 있는 전략을 도입해야 한다.

 오렌지 색 파란색 결과는 두 개의 마르코프 체인 샘플링을 나타낸다. 위의 pm.sample에서 chains인자를 4개 설정했다면, 4개의 색이 보일 것이다.

 이 외에도, arviz라는 패키지를 함께 이용하면 다양한 결과를 볼 수 있는데, 다음과 같이 forest plot을 뽑는다거나, summary를 출력할 수 있다.

코드)

import arviz as az


az.plot_forest(trace, var_names=["theta1",  "theta2", "sigma"], figsize=(5, 2), r_hat=True, combined=True, textsize=8)

코드)

sm = az.summary(trace, var_names=["theta1", "theta2", "sigma"], round_to=2)

arviz의 summary 출력은 pandas data frame 형태로 return된다.

 포레스트 플롯과 summary에서도 나오듯이 중요한 것이, 사후확률분포의 신뢰구간과, r_hat이라는 결관데, 신뢰구간은 파라미터 값의 불확실성을 나타내는 단위, r_hat은 MCMC샘플링의 신뢰도를 평가하기 위한 통계 값이다.

 r_hat은 체인 내부의 분산값과 전 MCMC샘플의 분산값의 비로 구하는 값인데, 이것이 1.1 이내로 오게 하는 것, 즉 전체 분산이 내부의 분산 보다 너무 크지 않게 샘플링 되는 것은 정상상태로 수속되었는지의 경험적 기준으로 삼는다고 한다. 

 이 외에도 정말 다양한 가시화와 모델링의 예시가 홈페이지에 있다

홈페이지 주소는 여기 -> https://docs.pymc.io/