#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter7: Section7.4.1
#
# P-spline regression model with Gaussian noise
#
# Author : Tomohiro Ando
#
################################################################################
library(splines)

# B-Splines design matrix

bspline <- function(x,xleft,xright,ndx,bdeg){
	dx <- (xright - xleft)/ndx
	knots <- seq(xleft - bdeg * dx , xright + bdeg * dx , by = dx)
	B <- spline.des(knots , x , bdeg +1 , 0 * x)$design
B
}


# Data

n <- 50
x <- seq(0,1,len=n)
z <- exp(-x*sin(2*pi*x))+1
y <- z+rnorm(n,0,0.3)

#Settings

p <- 10 # Number of basis functions

#lambda <- 100
#lambda <- 0.1
#lambda <- 0.001
lambda <- 0.00001

xleft  <- min(x)-0.0001
xright <- max(x)+0.0001
ndx    <- p-3
bdeg   <- 3

E <- bspline(x,xleft,xright,ndx,bdeg)
Dk <- diff(diag(1,p),differences=2); 
K <- t(Dk)%*%Dk; 
B <- solve(t(E)%*%E+n*lambda*K)%*%t(E)%*%y
yhat <- E%*%B


#Plot

plot(x,y,xlab="",ylab="",ylim=c(0,4))
lines(x,z)
lines(x,yhat,lty=2)