[Statistics] Bayesian Linear Regression

저번 포스팅에서 Bayesian 통계학과 회귀분석에 대해서 알아보았죠? 이번 포스팅과 또 다음 포스팅을 통해서 Bayesian 접근법의 구체적인 케이스에 대해서 알아보겠습니다. (깁스샘플링의 컨셉을 간단하게 다시 복습해봅시다. 결합확률분포를 직접적으로 구하는것이 어렵거나 불가능하지만 조건부 분포를 알고 있을 때 이를 활용하여 결합확률분포에서 추출된 샘플을 표집할 수 있다는 것입니다. 그리고 이 과정에서 iteration을 돌 때 관심 대상 변수 외에 나머지 변수는 고정된 것으로 간주하고 업데이트를 진행합니다)

 

먼저 Bayesian Lienar Regression입니다. 우리는 N개의 독립변수/종속변수 쌍을 가지고 있으며 이것을 선형 회귀모델로 모형화하려고 합니다. 

 

그런데 추정하고자 하는 parameter β, τ등은 Bayesian 관점에서 분포입니다. Conjugate 조건을 만족시키기 위해(그래야 analytic form으로 풀이가 가능합니다) 이 분포들을 다음과 같이 Normal과 Gamma 분포로 가정합니다. 

 

Gibbs sampling 은 추정하고자 하는 변수를 업데이트할 때 나머지 변수는 고정한 상태에서 업데이트가 이루어집니다. 

 

1) 업데이트:  에 대해서 아래와 같은 조건부 분포를 구해야합니다. 

 

 

여기 오른쪽 항에서 두번째 그러니까 prior부분은 이미 윗 줄에서 Normal 분포로 정의를 했구요, 첫번째 likelihood 부분은 우리가 당초에 살펴봤던 simple linear regression의 likelihood로 보면 되겠습니다. (전 포스팅의 parameter a가 , b가 이라고 생각하시면 됩니다)

 

log가 붙어있는 꼴에서 두 항을 잘 교통정리하면

 

와 같고, 이것을 quadratic form으로 잘 묶어주면 결국에 우리가 이 단계에서 원하는 parameter의 조건부 사후분포(Conditional Posterior distribution)를 얻을 수가 있습니다. (Norml X Normal 이었으니까 그 결과도 Normal의 꼴입니다)

 

 

2) 업데이트: 마찬가지로 정리하면 똑같은 조건부 사후분포를 얻을 수 있습니다.  

 

 

3)  업데이트:

이번에는 조건부 사후분포가 약간 다른 꼴입니다.  

 

 

tau에 대한 사전분포(prior)의 hyper-parameter가 라고 하고 log-likelihood꼴에서 구하는 조건부 사후분포는 아래와 같아집니다. prior에서 온 부분 제외하면 위의 식과 비슷하죠?

 

 

이와 같은 업데이트를 계속해서 1)~3) 반복해 나가다보면 우리가 구하는 분포의 참값에 다가가게 된다는 것입니다. 

 

아래 사이트에서 Gibbs sampling의 파이썬 코드를 퍼왔습니다. (자세한 내용은 이 페이지를 방문해보세요)

https://kieranrcampbell.github.io/blog/2016/05/15/gibbs-sampling-bayesian-linear-regression.html

Gibbs sampling for Bayesian linear regression in Python | Kieran R Campbell - blog

beta_0_true = -1
beta_1_true = 2
tau_true = 1

N = 50
x = np.random.uniform(low = 0, high = 4, size = N)
y = np.random.normal(beta_0_true + beta_1_true * x, 1 / np.sqrt(tau_true))

synth_plot = plt.plot(x, y, "o")
plt.xlabel("x")
plt.ylabel("y")
# simulation을 위한 데이터 생성

# Gibbs sampling을 위한 초기값 설정/hyper-parameter설정
# Markov chain 부분에서도 말씀드렸지만 초기 스타트 위치 자체는 크게 중요하지 않습니다. 
# 결국 일정 시간이 지나면 어디서 출발했든 참값에 가까이 가게 됩니다

## specify initial values
init = {"beta_0": 0,
        "beta_1": 0,
        "tau": 2}

## specify hyper parameters
hypers = {"mu_0": 0,
         "tau_0": 1,
         "mu_1": 0,
         "tau_1": 1,
         "alpha": 2,
         "beta": 1}

def sample_beta_0(y, x, beta_1, tau, mu_0, tau_0):
    N = len(y)
    assert len(x) == N
    precision = tau_0 + tau * N
    mean = tau_0 * mu_0 + tau * np.sum(y - beta_1 * x)
    mean /= precision
    return np.random.normal(mean, 1 / np.sqrt(precision))
    
def sample_beta_1(y, x, beta_0, tau, mu_1, tau_1):
    N = len(y)
    assert len(x) == N
    precision = tau_1 + tau * np.sum(x * x)
    mean = tau_1 * mu_1 + tau * np.sum( (y - beta_0) * x)
    mean /= precision
    return np.random.normal(mean, 1 / np.sqrt(precision))

def sample_tau(y, x, beta_0, beta_1, alpha, beta):
    N = len(y)
    alpha_new = alpha + N / 2
    resid = y - beta_0 - beta_1 * x
    beta_new = beta + np.sum(resid * resid) / 2
    return np.random.gamma(alpha_new, 1 / beta_new)


# Gibbs sampling의 몸통부분
def gibbs(y, x, iters, init, hypers):
    assert len(y) == len(x)
    beta_0 = init["beta_0"]
    beta_1 = init["beta_1"]
    tau = init["tau"]
    
    trace = np.zeros((iters, 3)) ## trace to store values of beta_0, beta_1, tau
    
    for it in range(iters):
        beta_0 = sample_beta_0(y, x, beta_1, tau, hypers["mu_0"], hypers["tau_0"])
        # 자 여기가 중요합니다. 지금 보시면 샘플링 된 beta_0 를 그다음 샘플링에 넣어주고 있습니다
        beta_1 = sample_beta_1(y, x, beta_0, tau, hypers["mu_1"], hypers["tau_1"])
        # tau를 뽑을 때도 앞단계에서 뽑혀나온 beta_0, beta_1을 그대로 넣어줍니다.
        tau = sample_tau(y, x, beta_0, beta_1, hypers["alpha"], hypers["beta"])
        trace[it,:] = np.array((beta_0, beta_1, tau))
        
    trace = pd.DataFrame(trace)
    trace.columns = ['beta_0', 'beta_1', 'tau']
        
    return trace

# parameter를 추출하고 다음번에 넣어주고 이 작업을 반복해서 1000번을 수행하면
# 이제 어느정도 stationary 상태에 도달하고 구하고자 하는 parameter도 참값에 가까이 가게 됩니다.
iters = 1000
trace = gibbs(y, x, iters, init, hypers)

traceplot = trace.plot()
traceplot.set_xlabel("Iteration")
traceplot.set_ylabel("Parameter value")

그럼, MCMC를 활용해서 Bayesian Regression을 추정하는 작업에 대한 포스팅은 여기서 마치겠습니다.