베이지안 회귀 분석 예제

import warnings
warnings.simplefilter('ignore')

import pymc3 as pm
import numpy as np
np.random.seed(1000)
import matplotlib as mpl
import matplotlib.pyplot as plt

단순 선형 회귀

• np.polyfit 명령으로 단순 선형 회귀 가능

x = np.linspace(0, 10, 500)
y = 4 + 2 * x + np.random.standard_normal(len(x)) * 2

plt.scatter(x, y)
plt.show()

reg = np.polyfit(x, y, 1)
reg

array([2.03384161, 3.77649234])

plt.scatter(x, y, c=y, marker='v', cmap=mpl.cm.jet)
plt.plot(x, reg[1] + reg[0] * x, lw=2.0)
plt.colorbar()
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

베이지안 회귀

with pm.Model() as model: 
    # PyMC3의 모형은 with 문 안에서 사용된다.
    # 사전 확률 정의
    alpha = pm.Normal('alpha', mu=0, sd=20)
    beta = pm.Normal('beta', mu=0, sd=20)
    sigma = pm.Uniform('sigma', lower=0, upper=10)
    
    # 선형 회귀 모형 정의
    y_est = alpha + beta * x
    
    # 가능도 분포 정의 
    likelihood = pm.Normal('y', mu=y_est, sd=sigma, observed=y)
    
    # 최적화를 사용하여 시작값 추정
    start = pm.find_MAP()
 
    # NUTS MCMC 샘플링 알고리즘 인스턴스 생성
    step = pm.NUTS(scaling=start)
    
    # 샘플링을 사용하여 100개의 사후 샘플 생성
    trace = pm.sampling.sample(100, step=step, start=start, progressbar=False)
    

logp = -1,068.5, ||grad|| = 60.625: 100%|██████████| 28/28 [00:00<00:00, 1040.77it/s]  
Only 100 samples in chain.
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta, alpha]

trace[0]

{'alpha': 3.7809674403284026,
 'beta': 2.0348240447962915,
 'sigma_interval__': -1.3284519594128563,
 'sigma': 2.0941554456646196}

pm.traceplot(trace, lines={'alpha': 4, 'beta': 2, 'sigma': 2}, figsize=(8, 8))
plt.show()

plt.scatter(x, y, c=y, marker='v', cmap=mpl.cm.jet)
plt.colorbar()
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
for i in range(len(trace)):
    plt.plot(x, trace['alpha'][i] + trace['beta'][i] * x)