#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter3: Section3.1.2
#
# Comparison of the constructed density function of the posterior mode
#
# Author : Tomohiro Ando
#
################################################################################

#library

library(MASS)        
library(KernSmooth)  

# Settings

n <- 50
x <- seq(0,1,len=n)
z <- 0.5+0.35*x
TRUEP <- exp(z)/(1+exp(z))
E <- cbind(1,x)

# Setting hyperparameters

lambda <- 0.001


# Matrix stores the estimated posterior mode

BB <- matrix(0,200,2)

# Repeat

for(RR in 1:200){

y <- as.vector(rep(0,len=n))

for(i in 1:n){
a <- runif(1,0,1)
if(a<=TRUEP[i]){y[i] <- 1}
}

Bmode <- runif(2,-0.1,0.1)
for(k in 1:100){
Bmode.old <- Bmode; O <- E%*%Bmode; Pi <- as.vector(exp(O)/(exp(O)+1))
W <- diag(Pi*(1-Pi)); Eta <- (y-Pi)/(Pi*(1-Pi))+O
Bmode <- solve(t(E)%*%W%*%E+lambda*diag(1,2))%*%t(E)%*%W%*%Eta
if(sum(abs(Bmode-Bmode.old))<=10^(-20)){break}
}

BB[RR,] <- t(Bmode)
}

x <- BB[,1]           
y <- BB[,2]           
f1 <- bkde2D(cbind(x, y), bandwidth=c(width.SJ(x), width.SJ(y)))
persp(f1$fhat,phi = 30, theta = 20, d = 5, col="lightblue", ticktype = "detailed", xlab="B0",ylab = "B1", zlab ="",zlim=c(0,0.1))




# Plot

x <- BB[,1]         
y <- BB[,2]         
f1 <- bkde2D(cbind(x, y), bandwidth=c(width.SJ(x), width.SJ(y)))
persp(f1$fhat,phi = 30, theta = 20, d = 5, col="lightblue", ticktype = "detailed", xlab="B0",ylab = "B1", zlab ="",zlim=c(0,0.1))