#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter4: Section4.6.1
#
# A direct Monte Carlo algorithm for linear regression model
#
# Author : Tomohiro Ando
#
################################################################################

library(MCMCpack)

# Data generation

n <- 30
p <- 5

bt<- c(-0.25,0.5,0,0,0)
x <- matrix(runif(n*p,-2,2),n,p)
z <- bt[1]*x[,1]+bt[2]*x[,2]
y <- z+rnorm(n,0,0.3)


# NUmber of draws

L  <- 10000

# Setting hyperparameters

lambda0 <- 10^-10
nu0     <- 10^-10
A       <- 10^-5*diag(1,2)

# Posterior inference

k <- c(1,2) 
X <- as.matrix(x[,k])
Bmle      <- solve(t(X)%*%X)%*%t(X)%*%y
Ahat      <- solve(A+t(X)%*%X)
nuhat     <- n+nu0
lambdahat <- lambda0+sum((y-X%*%Bmle)^2)+t(Bmle)%*%solve(solve(t(X)%*%X)+solve(A))%*%Bmle
Bhat <- solve(t(X)%*%X+A)%*%t(X)%*%y

p <- ncol(X)

for(i in 1:L){
sigma <- rinvgamma(1,shape=nuhat/2,scale=lambdahat/2)
Beta  <- mvrnorm(n=1,Bhat,sigma*solve(t(X)%*%X+A))
print(c(Beta,sigma))
}