#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter2: Section2.8.1
#
# Subset variable selection
#
# Author : Tomohiro Ando
#
################################################################################

# 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)

# Picking up predicotrs

#k <- 1 #Model M1
k <- c(1,2) #Model M2
#k <- c(1,2,3,4,5) #Model M3
#k <- c(1,3) #Model M4
#k <- c(2,4,5) #Model M5
#k <- c(3,4,5) #Model M6

X <- as.matrix(x[,k])
p <- ncol(X)

# Setting hyperparameters

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

# Bayes factor calculation
 
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
BayesFactor <- sqrt(det(Ahat))*sqrt(det(A))*(lambda0/2)^(nu0/2)*gamma(nuhat/2)/ ( (pi^(0.5*n))*gamma(nu0/2)*(lambdahat/2)^(nuhat/2) )

B  <- solve(t(X)%*%X+A)%*%t(X)%*%y
TE <- sum((y-x[,k]%*%matrix(B,ncol=1))^2)/n
PE <- sum((bt[1]*x[,1]+bt[2]*x[,2]-x[,k]%*%matrix(B,ncol=1))^2)/n

print(c(TE,PE,log(BayesFactor)))

# Plot surface

library(akima)

XX <- matrix(runif(1000*5,-2,2),1000,5)
z <- XX[,k]%*%matrix(B,ncol=1)
aa <- interp(XX[,1],XX[,2],z)
persp(aa$x,aa$y,aa$z, theta=30, phi=30, expand=1, r=6, col="lightblue", shade = 0.5, ticktype = "detailed", xlab = "X1",ylab = "X2", zlab ="f(X)",zlim=c(-2,2))