Bayesian 가우시안 기저 모형을 활용한 Linear 예측(R)

예제 가우시안 기저 모형

이성령

---

Prior

$p(w|\alpha) = N(w|0 , \alpha^{-1}I)$

---

Likelihood

$p(t|X,w,\beta) = \prod_{i=1}^N N(t_i | w^{T}\phi(x_i), \beta^{-1})$
$p(t|w,\beta) = p(t|x,w,\beta)$

---

Posterior

$p(w|t,\alpha , \beta) = \frac{p(t|w,\beta)p(w|\alpha)}{p(t|\alpha,\beta)}$ $= \frac{p(t|x,w,\beta)p(w|\alpha)}{p(t|\alpha,\beta)}$
$p(w|t,\alpha , \beta) = N(w | m_N , S_N)$

$m_{N} = \beta S_N \phi^{T}t$
$S_N^{-1} = \alpha I + \beta \phi^{T}\phi$

---

Posterior Predictive

$p(t|w,\beta) = p(t|x,w,\beta)$

$p(t^{*}|x,t,\alpha,\beta) = \int p(t^{*}|x,w,\beta)p(w|t,\alpha,\beta)dw$

$\quad \quad = N(t|m_N^{T}\phi(x) , \sigma^2_{N}(x))$

$\sigma^2_{N}(x) = \frac{1}{\beta} + \phi(x)^T S_N\phi(x)$

$p(t^{*}|X,w,\beta) = N(t^{*} | w^{T}\phi(x), \beta^{-1})$
$p(w|t,\alpha , \beta) = N(w | m_N , S_N)$

---

예제 가우시안 기저 모형

• 실제 DGP : $y= sum(2\pi x) + \epsilon$ , $x \sim U(0,1)$ , $\epsilon \sim N(0, \frac{1}{\beta})$
• 가우시안 기저 모형 : $sin(2\pi x) \approx \theta_0 + \theta_1\phi_1(x)+ \theta_2\phi_2(x) + ... + + \theta_9\phi_9(x)$
  • 단 $\phi_k(x) \sim N(\frac{k}{9} , \frac{1}{8}^2) , (k = 1,2,...,9)$
• Prior : $\theta = (\theta_0, \theta_1 , ..., \theta_9)^T \sim N(0, \alpha^{-1} \boldsymbol{I})$
• $\alpha$ = 2 , $\beta$ = 25
• N = 1, 2, 20

DGP <- function(X , noise_variance){
  return( 0.5 + sin(2 * pi * X ) + 
            matrix( rnorm(dim(X)[1]*dim(X)[2],mean=0,sd=noise_variance), dim(X)[1] , dim(X)[2])  )
}


Gaussian_Basis <- function(x , mu , sigma = 1/8){
  return( exp(-0.5*(x-mu)**2 / sigma**2))
}


expand <- function( x, bf , bf_args){
  if(is.null(dim(x))){
    a <- matrix(1, nrow =  1 , ncol = 1)
  }else{
    a = matrix(1, nrow =  dim(x)[1] , ncol = dim(x)[2])  
  }
  
  for( i in bf_args ){
    a <- cbind(a , bf(x , i))
  }
  return(a)
}

시각화

N_list = c(1,2,20 , 50)

beta = 25.0
alpha = 2.0

X = matrix( runif(1*N_list[length(N_list)],min = 0 , max = 1), 1 , N_list[length(N_list)])
t = DGP(X , noise_variance =  1/ beta)
X_test = matrix(seq(from = 0.0, to = 1.0, length.out =  100), 100, 1)
Phi_test = expand(X_test, bf=Gaussian_Basis,
                  bf_args=seq(1/9,1,length.out = 9))
y_true = DGP(X_test, noise_variance=0)

Visualization Loop

i = 1 
p <- list()
for( N in N_list ){
  X_N = matrix(X[0:N] , ncol = 1)
  dim(X_N)
  t_N = matrix(t[0:N] , ncol =  1 )
  dim(t_N)
  Phi_N = expand( X_N, bf=Gaussian_Basis, bf_args=seq(1/9 , 1 ,length.out = 9))
  dim(Phi_N)
  post = posterior(Phi_N, t_N, alpha, beta)
  m_N <- post$mu 
  dim(m_N)
  m_N <- matrix(m_N , ncol = 1)
  S_N <- post$sigma
  m_N <- matrix(m_N , ncol = 1, nrow = length(m_N))
  dim(S_N)
  OUTPUT = posterior_predictive(Phi_test, m_N, S_N, beta)
  pred <- OUTPUT$y
  std <- sqrt(OUTPUT$y_var)
  d = data.frame(X_test , pred , std , y_true)
  out <- ggplot(d) + 
    geom_ribbon(aes(X_test ,ymin=pred-std,ymax=pred+std),
                color="yellow",alpha=0.1) +
    geom_line(aes(X_test, pred, color = "darkblue"), size= 2) +
    geom_line(aes(X_test, y_true, color ="red") , size= 2) +
    scale_color_discrete(name = " ", labels = c("Prediction", "TRUE")) +
    theme_bw() +
    ggtitle(TeX(paste0("$sin(2\\pi x)+$", "$\\epsilon$  ", "(N : ",N, ")" ) )) +
    labs(x = "Data Range" , y = "") + 
    theme(plot.title = element_text(size = 25 , face = "bold" , hjust = 0.5) , 
          legend.position="top",
          axis.line = element_line(colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank()) +
    geom_point(data = data.frame(X_N , t_N) ,aes(X_N , t_N))
    
  p[[i]] <- out
  i = i + 1
}

Visualization

do.call(grid.arrange,p)

파이썬 코드를 R 코드로 변환

• http://krasserm.github.io/2019/02/23/bayesian-linear-regression/