blm

#------------------------------------------------------
# bayesian linear model
#------------------------------------------------------
blm <- function(formula, data, lambda){
  #model
  mm <- model.matrix(formula, data) #design matrix
  y <- model.response(model.frame(formula, data)) #response variable
  
  #hyper parameter
  Lambda <- diag(ncol(mm))
  m <- rep(0, ncol(mm))
  
  #learning
  Lambda_hat <- lambda*(t(mm) %*% mm) + Lambda
  Lambda_hat_inv <- solve(Lambda_hat)
  m_hat <- Lambda_hat_inv %*% t((lambda * (y %*% mm) + t(Lambda %*% m)))
  
  #model evidence
  mde = -0.5*(lambda * sum(y^2) - log(lambda) + log(2*pi) + as.vector(m %*% Lambda %*% m) - 
                log(det(Lambda)) - as.vector(t(m_hat) %*% Lambda_hat %*% m_hat) + log(det(Lambda_hat)))
  
  return(list(formula = formula, lambda = lambda, Lambda_hat = Lambda_hat, Lambda_hat_inv = Lambda_hat_inv, m_hat = m_hat,
              model_evidence = mde))
}
#predict
predict.blm <- function(model, data, type = "response"){
  mm <- model.matrix(model$formula, data)
  
  if(type == "response"){
    m_ast <- t(model$m_hat) %*% t(mm)
    ret <- as.vector(m_ast)
  }else if(type == "sd"){
    variance_ast <- model$lambda^-1 + apply(mm %*% model$Lambda_hat_inv * mm, 1, sum)
    ret <- as.vector(sqrt(variance_ast))
  }
  return(ret)
}

#----------------------------------------------------
#demo 
#----------------------------------------------------
X <- runif(10, 0, 10)
y <- sin(X)
df <- data.frame(X, y)

(model <- blm(y ~ ., df, 10))
(model <- blm(y ~ X + I(X^2), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) + I(X^6), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) + I(X^6) + I(X^7), df, 10))

df_test <- data.frame(X = runif(100, 0, 10), y = 0)
df_test <- df_test[order(df_test$X), ]

yhat <- predict.blm(model, df_test, type = "response")
sq_hat <- predict.blm(model, df_test, type = "sd")

#plot
plot(df$X, df$y)
lines(df_test$X, yhat)
lines(df_test$X, yhat + 2*sq_hat)
lines(df_test$X, yhat - 2*sq_hat)

blm2.r

#------------------------------------------------------
# bayesian linear model
#------------------------------------------------------
blm <- function(formula, data, lambda){
  #model
  mm <- model.matrix(formula, data) #design matrix
  y <- model.response(model.frame(formula, data)) #response variable
  
  #hyper parameter
  Lambda <- diag(ncol(mm))
  m <- rep(0, ncol(mm))
  
  #learning
  Lambda_hat <- lambda*(t(mm) %*% mm) + Lambda
  Lambda_hat_inv <- solve(Lambda_hat)
  m_hat <- Lambda_hat_inv %*% t((lambda * (y %*% mm) + t(Lambda %*% m)))
  
  #model evidence
  mde = -0.5*(lambda * sum(y^2) - log(lambda) + log(2*pi) + as.vector(m %*% Lambda %*% m) - 
                log(det(Lambda)) - as.vector(t(m_hat) %*% Lambda_hat %*% m_hat) + log(det(Lambda_hat)))
  
  return(list(formula = formula, lambda = lambda, Lambda_hat = Lambda_hat, Lambda_hat_inv = Lambda_hat_inv, m_hat = m_hat,
              model_evidence = mde))
}
#predict
predict.blm <- function(model, data, type = "response"){
  mm <- model.matrix(model$formula, data)
  
  if(type == "response"){
    m_ast <- t(model$m_hat) %*% t(mm)
    ret <- as.vector(m_ast)
  }else if(type == "sd"){
    variance_ast <- model$lambda^-1 + apply(mm %*% model$Lambda_hat_inv * mm, 1, sum)
    ret <- as.vector(sqrt(variance_ast))
  }
  return(ret)
}

#----------------------------------------------------
#demo 
#----------------------------------------------------
X <- runif(10, 0, 10)
y <- sin(X)
df <- data.frame(X, y)

(model <- blm(y ~ ., df, 10))
(model <- blm(y ~ X + I(X^2), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) + I(X^6), df, 10))
(model <- blm(y ~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) + I(X^6) + I(X^7), df, 10))

df_test <- data.frame(X = runif(100, 0, 10), y = 0)
df_test <- df_test[order(df_test$X), ]

yhat <- predict.blm(model, df_test, type = "response")
sq_hat <- predict.blm(model, df_test, type = "sd")

#plot
plot(df$X, df$y)
lines(df_test$X, yhat)
lines(df_test$X, yhat + 2*sq_hat)
lines(df_test$X, yhat - 2*sq_hat)