##################################################################
# Smoothing spline

# define a smooth 'true' function

# a sine curve
fun = function(x) {1 + 3*sin(2*pi*x-pi)}

# generate random x, and y as a noisy version of fun(x);
# plot points, and the true curve on a fine grid

set.seed(200)
x = sort(runif(100))
y = fun(x)+rnorm(100)
xgrid = seq(0,1, by=.0025)
ygridt = fun(xgrid)
plot(x,y)
lines(xgrid,ygridt,lty=2)

postscript("scatterplot.ps", width=4, height=4, horizontal=F)
plot(x,y, col=4)
dev.off()

# fit a cubic smoothing spline: overfit
# Read R help to know the relationship between "spar" and "lambda"
fit1 = smooth.spline(x,y, spar=1/4)
postscript("smooth1.ps", width=4, height=4, horizontal=F)
plot(x,y, col=4)
lines(predict(fit1,xgrid))
dev.off()

# fit a cubic smoothing spline: underfit
fit2 = smooth.spline(x,y, spar=2)
postscript("smooth2.ps", width=4, height=4, horizontal=F)
plot(x,y, col=4)
lines(predict(fit2,xgrid))
dev.off()

# fit a cubic smoothing spline using the default (GCV) tuning parameter
spl.fit = smooth.spline(x,y)
postscript("smooth3.ps", width=4, height=4, horizontal=F)
plot(x,y, col=4)
lines(predict(spl.fit,xgrid))
dev.off()

# compare the fitted function with the true function
postscript("truevssmooth.ps", width=4, height=4, horizontal=F)
plot(x,y, type="n")
lines(xgrid,fun(xgrid),col=4)
lines(predict(spl.fit,xgrid))
dev.off()

######################################################################
# Interpolating spline

library(splines)

# generate points
set.seed(100)
n = 5
x = (0:n)/n
y = rnorm(n+1)

postscript("splinebend.ps", width=5, height=4, horizontal=F)

#plot(x, y, main = paste("Spline through", n+1, "points"),ylim=c(-1,1),cex=2)
#lines(spline(x, y, n = 201, method = "natural"), col = 2)

plot(x, y, main = paste("Spline through", n+1, "points"),
     xlim=c(-.1,1.1), ylim=c(-1,1),cex=2)
lines(spline(x, y, xmin=-.1, xmax=1.1, n = 201, method = "natural"), col = 2)
abline(v=(0:n)/n, lty=2, col=3)
dev.off()

ispl <- interpSpline(y ~ x)

postscript("splineg.ps", width=8, height=6, horizontal=F)
par(mfrow = c(2, 2), las = 1)

## plot over the range of the knots: [first knot, last knot]

#plot(predict(ispl, nseg = 201),    
#     main = "Interpolating spline", type = "l", xlab = "x", ylab = "y") 
#points(x, y, col = 2)
#abline(v=(1:(n-1))/n, lty=2, col=3)
#plot(predict(ispl, nseg = 201, deriv = 1),
#     main = "First derivative", type = "l", xlab = "x", ylab = "") 
#abline(v=(1:(n-1))/n, lty=2, col=3)
#plot(predict(ispl, nseg = 201, deriv = 2),
#     main = "Second derivative", type = "l", xlab = "x", ylab = "")
#abline(v=(1:(n-1))/n, lty=2, col=3) 
#plot(predict(ispl, nseg = 401, deriv = 3),
#     main = "Third derivative", type = "l", xlab = "x", ylab = "")
#abline(v=(1:(n-1))/n, lty=2, col=3) 

# include the boundaries

plot(predict(ispl, seq(-.1,1.1,length.out= 201)),    
     main = "Interpolating spline", type = "l", xlab = "x", ylab = "y") 
points(x, y, col = 2)
abline(v=(0:n)/n, lty=2, col=3)
plot(predict(ispl, seq(-.1,1.1,length.out= 201), deriv = 1),
     main = "First derivative", type = "l", xlab = "x", ylab = "") 
abline(v=(0:n)/n, lty=2, col=3)
plot(predict(ispl, seq(-.1,1.1,length.out= 201), deriv = 2),
     main = "Second derivative", type = "l", xlab = "x", ylab = "")
abline(v=(0:n)/n, lty=2, col=3) 
plot(predict(ispl, seq(-.1,1.1,length.out= 401), deriv = 3),
     main = "Third derivative", type = "l", xlab = "x", ylab = "")
abline(v=(0:n)/n, lty=2, col=3) 

dev.off()


########################################################################
#  B-splines  

cubicBs = bs(xgrid, knots=(1:9)/10, intercept=T)
quadBs = bs(xgrid, knots=(1:9)/10, degree=2, intercept=T)
linBs = bs(xgrid, knots=(1:9)/10, degree=1, intercept=T)


postscript("Bi2.ps", width=5, height=3, horizontal=F)
plot(xgrid,linBs[,4],col=1,type="l",lwd=2,xlab="x",ylab="", 
main=expression(B["2,2"]))
lines(xgrid,linBs[,3], col="gray30", lty=1)
lines(xgrid,linBs[,5], col="gray30", lty=1)
abline(v=(1:9)/10,col=3,lty=2)
dev.off()

postscript("Bsplines.ps", width=8, height=6, horizontal=F)
par(mfrow=c(2,1))
plot(xgrid,quadBs[,5],col=1,type="l",lwd=2,xlab="x",ylab="", ylim=c(0,1), main=expression(B["2,3"]))
lines(xgrid,linBs[,4], col="gray30", lty=1)
lines(xgrid,linBs[,5], col="gray30", lty=1)
abline(v=(1:9)/10,col=3,lty=2)

plot(xgrid,cubicBs[,6],col=1,type="l",lwd=2,xlab="x",ylab="", ylim=c(0,1), main=expression(B["2,4"]))
lines(xgrid,quadBs[,5], col="gray30", lty=1)
lines(xgrid,quadBs[,6], col="gray30", lty=1)
abline(v=(1:9)/10,col=3,lty=2)
dev.off()