3. Bayesian Linear Regression

1. Frequentist 방식

 
   다음의 예측 식이 있다고 하자.
$$
f(\mathbf{X}) = \beta_0 + \sum_{j=1}^{p}\mathbf{X}_j\beta_j + \epsilon \\
= \beta^T \mathbf{X} + \epsilon
$$
   여기서 $\epsilon \sim N(0, \sigma^2)$인 에러텀이다. 여기서 학습을 위한 데이터가 $(x_1, y_1), ...(x_N, y_N)$과 같이 있다고 하면, 우리는 이 데이터를 가장 비슷하게 맞추는 $\beta$를 찾아야 한다. 이를 위해 일반적으로 OLS(ordinary least squares)방식을 사용한다. OLS 방식을 사용하기 위해 우선 측정 값과 $\beta$를 통한 예측치간의 차이를 나타내는 다음의 식을 정의한다.
$$
RSS(\beta) = \sum^{N}_{i=1}(y_i - f(x_i))^2
$$
   이 식(RSS: residual sum of squares)는 일종의 비용(에러)함수로 볼 수 있다. OLS의 목적은 이 함수의 값을 최소화하는 $\beta$를 찾는 것이고, 이를 위해 Maximum Likelihood Estimate 방식을 사용한다. 이 식은 다음과 같다.
$$
\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{y}
$$
   이렇게 찾아진 $\hat{\beta}$를 이용하여 $x_{N+1}$ 값에 대한 결과치 $y_{N+1}$을 추정한다.

2. Bayesian 방식

   Bayesian 방식으로 Linear Regression을 수행하기 위해 다음과 같이 $\mathbf{y}$에 대해 정의한다.
$$
\mathbf{y} \sim N(\beta^T \mathbf{X}, \sigma^2 \mathbf{I})
$$
   즉, 앞서 소개한 방식과 다르게 $y$를 어떤 점으로 보는 것이 아니라, 일종의 확률변수로 보는 것이다. 이를 수행하기 위해 다음의 두 가지를 고려해야 한다.

• Prior Distribution - $\beta$에 대한 선행 정보가 있을 경우, 이를 반영하기 위한 분포를 지정한다. 만약 정보가 없을 경우 'non-informative prior'를 선택한다. 
• Posterior Distribution - Bayesian 방식은 $\beta$도 분포를 갖는 확률변수로 취급하므로 데이터 추가에 따른 새로운 $\beta$의 분포를 도출한다. 

   위에서 등장하는 non-informative prior는 추가 정보를 찾아봐야 이해할 것 같다. 지금 여기선 잠시 넘어간다.
   앞선 포스트에서 설명하였듯이, $\beta$의 분포가 conjugate prior인 경우 닫힌 해를 이용하여 다음 $\beta$의 분포를 추정할 수 있다. 그렇지 않을 경우 일정 분포를 가정하고 Markov Chain Monte Carlo 방법을 이용하여 추정작업을 수행한다.

3. PyMC3를 이용한 Bayesian Linear Regression 예제

   너무나도 간단한 설명을 마치고 PyMC3를 이용하는 예제를 소개한다.

 3.1 Generalized Linear Model(GLM)

   GLM 모형은 일반적인 Linear Regression 모형의 확장형태라고 볼 수 있다. 결과 값으로 나오는 $y$의 오류값이 정규분포가 아닌 다른 분포를 갖을 수 있도록 허용한 방식이다. 자세한 내용은 위키페이지를 참고하기 바란다.
   GLM 모형은 결과값 $y$가 지수형태를 갖는 특정 분포를 통해 만들어진다고 가정한다. $y$는 다음과 같은 평균값을 갖는다.
$$
\mathbf{E}(\mathbf{y}) = \mu = g^{-1}(\mathbf{X}\beta)
$$
   여기서 $\mathbf{X}\beta$가 linear predictor 역학을 하는데, $\beta$가 미지수이다. $g()$는 연결함수(link function)이다. 그리고 분산은 평균의 함수 $V$로 표현한다.
$$
Var(\mathbf{y}) = V(\mathbf{E}(\mathbf{y})) = V(g^{-1}(\mathbf{X}\beta))
$$
   이 모형을 가지고 $\beta$를 다양한 방식으로 추정하는데, 여기서는 Markov Chain Mote Carlo를 이용한다.

  3.2 수행 예제

   다음의 예제는 위의 GLM모형을 임의의 데이터에 적용하는 예제이다. 난수로 만들어진 일정 데이터에 맞게 $\beta$를 찾는 코드이다. 우선 다음의 그림은 frequentist 방식으로 fitting한 결과이다.

   이 작업을 MCMC를 이용한 결과는 다음과 같다. 다음의 그림은 5000번의 수행 중, 100개의결과를 임의 추출한 것이다.

   위 그림의 녹색선이 목표지점이다. 이 MCMC수행의 추적 차트는 다음과 같다.

   위 그림에서 intercept가 $\beta_0$, $x$가 $\beta_1$, sd가 error이다.
   이를 수행하는 코드를 다음에 첨부한다.

 
# simul_fit.py

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import pymc3 as pm

sns.set(style="darkgrid", palette="muted")

# simulate function
def simulate_linear_data(N, beta_0, beta_1, eps_sigma_sq):
    """
    simulate a random dataset using a noisy linear proces

    :param N: Number of data points
    :param beta_0: Intercept
    :param beta_1: slope of univariate predictor, X
    :param eps_sigma_sq:
    :return: pandas dataframe
    """

    # create pandas DataFrame
    # column 'x' : N uniformly sampled values b/w 0.0 and 1.0
    df = pd.DataFrame(
        { "x":
            np.random.RandomState(42).choice(
              list(map(
                  lambda x: float(x) / 100.0,
                  np.arange(N)
              )), N, replace=False
            )
        }
    )

    # Use linear model (y ~ beta_0 + beta_1*x + epsilon) to
    # generate a column 'y' of responses based on 'x'
    eps_mean = 0.0
    df["y"] = beta_0 + beta_1 * df["x"] + np.random.RandomState(42).normal(
        eps_mean, eps_sigma_sq, N)

    return df


def glm_mcmc_inference(df, iterations=5000):
    """
    Calculate MCMC trace of GLM
    :param df:
    :param iterations:
    :return:
    """

    # use PyMC3 to construct model context
    basic_model = pm.Model()
    with basic_model:
        # create the glm using the Patsy model syntax
        # Normal distribution for likelihood
        pm.glm.GLM.from_formula("y ~ x", df, family=pm.glm.families.Normal())

        # use Maximum A Posterior(MAP) optimization as initial value
        start = pm.find_MAP()

        # use the No-U-Turn Sampler
        step = pm.NUTS()

        # Calculate trace
        trace = pm.sample(iterations, step, start, random_seed=42, progressbar=True)

    return trace


# main
if __name__ == "__main__":
    # true parameters
    beta_0 = 1.0 # Intercept
    beta_1 = 2.0 # Slope

    # simulate 100 points
    N = 100
    eps_sigma_sq = 0.5

    # Simulate the linear data
    df = simulate_linear_data(N, beta_0, beta_1, eps_sigma_sq)

    # plot data, and Frequentist linear regression fit
    # using seaborn package
    sns.lmplot(x="x", y="y", data=df, size=10)
    plt.xlim(0.0, 1.0)

    trace = glm_mcmc_inference(df, iterations=5000)
    pm.traceplot(trace[500:])
    plt.show()

    # plot a sample of posterior regression lines
    sns.lmplot(x="x", y="y", data=df, size=10, fit_reg=False)
    plt.xlim(0.0, 1.0)
    plt.ylim(0.0, 4.0)
    pm.plot_posterior_predictive_glm(trace, samples=100)
    x = np.linspace(0, 1, N)
    y = beta_0 + beta_1 * x
    plt.plot(x, y, label="True Regression Line", lw=3., c="green")
    plt.legend(loc=0)
    plt.show()