#  This script should be in a series with the other spline figures.
##Tried to edit/format this to (hopefully) be more like the 
##rest (poly, s5vs15, loctestf) -Spencer

library(splines)

n.observations <- 101

x.obs <- seq(from = -2, to = 2, length.out = n.observations)
f.obs <- sin(x.obs)+2*exp(-30*(x.obs^2))
f.sd <- sd(f.obs)

added.noise.sd <- f.sd/50
added.noise <- rnorm(n.observations, 0, added.noise.sd)
observations <- f.obs+added.noise

x.knots <- c(-1.8, -0.4, -0.2, 0, 0.2, 0.4, 1.8)
x.mat <- matrix(0, nrow = length(x.obs), ncol = length(x.knots)+3)
x.mat[,1:3] <- cbind(x.obs, x.obs^2, x.obs^3)
for(i in 1:length(x.knots)){
	x.mat[,(i+3)] <- mapply(max, (x.obs-x.knots[i])^3, 0)
}

x.df <- data.frame(x.mat)

x.values = seq(-2, 2, by = 0.001)
x.matpred <- matrix(0, nrow = length(x.values), ncol = length(x.knots)+3)
x.matpred[,1:3] <- cbind(x.values, x.values^2, x.values^3)
for(i in 1:length(x.knots)){
	x.matpred[,(i+3)] <- mapply(max, (x.values-x.knots[i])^3, 0)
}

x.df.predictions <- data.frame(x.matpred)

splines.model <- lm(observations~., cbind.data.frame(observations, x.df))
splines.predictions <- predict(splines.model, x.df.predictions)

#  Prepare some graphical parameters.
label.size <- 1.5
axis.size <- 1.4
line.color <- 153

#  Plot the observations and spline model.
plot(x.obs, observations, pch = 16,
     xlim = c(-2, 2), ylim = c(-1.2, 2.2), 
     xaxt = "n", yaxt = "n", 
     main = "", xlab = "x", ylab = "y", cex.lab = label.size)
lines(x.values, splines.predictions, col = line.color, lwd = 2)

x.axis <- seq(from = -2, to = 2, by = 1)
y.axis <- seq(from = -1, to = 2, by = 1)
axis(1, at = x.axis, labels = x.axis, cex.axis = axis.size)
axis(2, at = y.axis, labels = y.axis, cex.axis = axis.size)

dev.print(device = postscript, "15.4.eps", horizontal = TRUE)