#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter5: Section5.6.1
#
# Nonlinear regression models using basis expansion predictors with GBIC
#
# Author : Tomohiro Ando
#
################################################################################

# library

library(mva)
library(stats)
library(akima)

# Data generation

n <- 400
d <- 2
x <- matrix(runif(d*n,0,2),ncol=2)
z <- sin(2*pi*x[,1])+cos(2*pi*x[,2])
y <- z+rnorm(n,0,1)

# Setting hyperparameters

m  <- 30       # number of basis fucntions
nu <- 5        # hyperparameters
beta <- 0.0001 # smoothing parameter

# Posterior inference

Dk <- diff(diag(1,m+1),differences=2)
K <- t(Dk)%*%Dk
Km <- kmeans(x,m)
U <- Km$centers
W <- Km$withinss
Size <- Km$size
IN <- matrix(1,nrow=n,ncol=1)
P <- matrix(0,nrow=n,ncol=m)
for(i in 1:d){P <- P + (x[,i]%*%matrix(1,nrow=1,ncol=m)-IN%*%t(U[,i]))^2}
SS    <- nu*W/Size
Sigma <- diag(1/(2*SS))
E     <- cbind(1,exp(-P%*%Sigma))            #Design matrix
Beta  <- solve(t(E)%*%E+n*beta*K)%*%t(E)%*%y #Posterior mode
O     <- E%*%Beta                            # Predicted value 
S     <- sqrt((sum((O-y)^2))/n)              #Sample standard devitations

lambda <- beta/(S^2)


# GBIC score

L <- diag(as.vector(y-O))
J.BB <- t(E)%*%E+n*beta*K
J.SB <- t(IN)%*%L%*%E/(S^2)
J.BS <- t(J.SB)
J.SS <- n/(2*S^2)
Second <- rbind(cbind(J.BB,J.BS),cbind(J.SB,J.SS))/(n*S^2)

psi1 <- eigen(Second)$values
psi2 <- eigen(K)$values
det1 <- prod(psi1)
det2 <- prod(psi2[1:(m-1)])

BIC <- (n+m-1)*log(S^2)+n*lambda*t(Beta)%*%K%*%Beta+n+(n-3)*log(2*pi)+3*log(n)+log(det1)-log(det2)-(m-1)*log(beta)


# Data plot

D <- interp(x[,1],x[,2],y)
persp(D,ticktype="detailed",phi=40,theta=45)


# True function

D <- interp(x[,1],x[,2],z)
persp(D,ticktype="detailed",phi=40,theta=45)


# Estimated function

D <- interp(x[,1],x[,2],O)
persp(D,ticktype="detailed",phi=40,theta=45)