R로 구현하는 Gibbs Sampling

R로 구현하는 Gibbs Sampling

 

Bayesian Linear Regression

\[y_i = \beta_0 + \beta_1 x_i + \epsilon, \epsilon \sim N(0,1/\tau) \] \[ y_i \sim N(\beta_0 + \beta_1 x_i, 1/\tau)\] 이 모델에서 아래 데이터를 관찰했을 때, likelihood는 \[ (y_i, x_i), i=1,2,3...,N \] \[ L(yi,xi | \beta_0, \beta_1, \tau) = \prod_{i=1}^{N}N(\beta_0, \beta_1 x_i, 1/\tau) \] Bayesian의 특징은 여기서 beta0과 beta1, tau의 prior distribution을 정의한다는 것이다. 예를 들어, conjugate prior를 가정한다.

 

\[\beta_0 = N(\mu_0, 1/\tau_0) \] \[\beta_1 = N(\mu_1, 1/\tau_1) \] \[\tau \sim Gamma(\alpha, \beta) \] 이 파라미터들의 posterior distribution을 구하는 것이 bayesian estimation의 목적이다. 현재 데이터를 관찰함으로써, prior distribution이 변한 것이 posterior이다.

 

Gibbs sampling

우리의 목적은 theta0과 theta1의 joint posterior distribution을 구하는 것이다. 즉, 아래 pdf를 구하는 것이다.

\[ p(\theta_1, \theta_2 || x) \] 이걸 구하기 위해서 Gibbs sampling 방법에서는 각각의 파라미터에 대해full conditional distribution을 구한다. (이것이 가장 힘든 파트이다.) 데이터와, 모든 파라미터를 상수로 가정하고, 관심 있는 파라미터의 조건부 pdf를 구하는 것이다.

 

1. 초기값을 정한다. \[ \theta_2^{(i)} \]

2. 표본을 아래 분포에서 뽑는다. \[\theta_1^{(i+1)} \sim p(\theta_1 | \theta_2^{(i)}, x) \]

3. 표본을 아래 분포에서 뽑는다. \[ \theta_2^{(i+1)} \sim p(\theta_2 | \theta_1^{(i)}, x) \]

그러면 위에 기술했던 Linear 모델의 posterior를 우선 구해보자.

우선 x,y data를 관찰했을 때의 beta0의 posterior는 아래와 같다.

\[ \beta_0|\beta_1,\tau,\tau_0,\mu_0,x,y \] 이것은 아래의 정규분포를 따른다. 왜 이렇게 되는지 자세한 설명은 이곳에 있다. 

\[ \beta_0|\beta_1,\tau,\tau_0,\mu_0,x,y \sim N(\frac{\tau_0\mu_0 + \tau \sum_i(y_i-\beta_1 x_i)}{\tau_0 + \tau N}, 1/(\tau_0 + \tau N)) \] 마찬가지로 나머지 추정하고 싶은 파라미터에 대한 posterior는 아래와 같다.

\[ \beta_1|\beta_0,\tau,\tau_0,\mu_0,x,y \sim N(\frac{\tau_1\mu_1 + \tau \sum_i(y_i-\beta_0 x_i)}{\tau_1 + \tau \sum_ix_i{^2}}, 1/(\tau_1 + \tau \sum_ix_i{^2}))\]

\[ \tau|\beta_0,\beta_1, \alpha, \beta ,x,y \sim Gamma(\alpha+\frac{N}{2}, \beta + \sum_i( \frac{(y_i-\beta_0-\beta_1x_i)^2}{2}) \] 이제 posterior를 알았으니, Gibbs sampling을 이용해 이를 추정해보자.

 

# Hyperparameters (Prior) 
k.mu0 <- 0
k.mu1 <- 0
k.tau0 <- 1
k.tau1 <- 1
k.alpha <- 2
k.beta <- 1
# Initial parameter (임의의 값)
k.beta0 <- 0
k.beta1 <- 0
k.tau <- 2

# Simulation을 위한 데이터 생성. 
beta.true <- 2
int.true <- -1
tau.true <- 1
x <- runif(min = 0, max = 4, 100)
y <- beta.true*x + int.true + rnorm(100, mean=0, sd=1/tau.true)
plot(x,y)

 

관찰한 데이터

 

full condition pdf에서 sample 추출하는 함수 정의

get.beta0.sample <- function(beta1, tau, tau0, mu0, x, y) { 
  N <- length(y) 
  precision <- tau0 + tau*N 
  mean <- tau0 * mu0 + tau * sum(y - beta1 * x) 
  mean <- mean/precision 
  result <- rnorm(1, mean=mean,sd=1/sqrt(precision)) 
  return(result) 
}
get.beta1.sample <- function(beta0, tau, tau1, mu1, x, y) {
  N <- length(y)
  precision <- tau1 + tau*sum(x*x)
  mean <- tau1 * mu1 + tau * sum((y - beta0) * x)
  mean <- mean/precision
  result <- rnorm(1, mean=mean,sd=1/sqrt(precision))
  return(result)
}
get.tau.sample <- function(x, y, beta0, beta1, alpha, beta){
  N <- length(y)
  alpha_new <- alpha + N / 2
  resid <- y - beta0 - beta1 * x
  beta_new <- beta + sum(resid * resid)/2
  return(rgamma(1, shape=alpha_new, scale=1/beta_new))
}

# Burn-in period를 제거하고 result를 저장
burnin <- 9000
result <- trace[-burnin,]

# Prior와 Posterior 비교하기
temp_x <- seq(from=-3,to=3,by = 0.01)
beta0.prior <- dnorm(x=temp_x, mean=k.mu0, sd=1/k.tau0)
beta0.post <- result[,'beta0']
plot(temp_x, beta0.prior)

 

hist(beta0.post)

 

temp_x <- seq(from=-3,to=3,by = 0.01)
beta1.prior <- dnorm(x=temp_x, mean=k.mu1, sd=1/k.tau1)
beta1.post <- result[,'beta1']
plot(temp_x, beta1.prior)

 

hist(beta1.post)

 

temp_x <- seq(from=0,to=10,by = 0.01)
tau.prior <- dgamma(x=temp_x, shape=k.alpha, scale = k.beta)
tau.post <- result[,'tau']
plot(temp_x, tau.prior)

hist(tau.post)

 

 

참고 

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