#  This script fits a 20th order polynomial to noisy observations of 
#     f(x) = sin(x) + 2*exp( - 30 * x^2)
#
#
#  Figure caption:  Data simulated from function f(x) = sin(x) +
#  2exp(-30*x*x) together with twentieth order polynomial fit
#  (shown as line).  Note that the polynomial is over-fitting
#  (under-smoothing) in the relatively smooth regions of f(x), and
#  under-fitting (over-smoothing) in the peak.  (In the data shown
#  here, the noise standard deviation is 1/50 times the standard
#  deviation of the function values.)

n.observations <- 101
polynomial.order <- 20

#  Set up observations and function values.
x.at.obs <-seq(from = -2, to = 2, length.out = n.observations)
f.at.obs <- sin(x.at.obs) + 2*exp(-30 * (x.at.obs^2))

#  Obtain a reasonable noise level and add noise.  
f.sd <- sd(f.at.obs)
relative.noise <- 1/50
added.noise.sd <- f.sd * relative.noise
added.noise <- rnorm(n.observations, 0, added.noise.sd)
observations <- f.at.obs + added.noise

#  Obtain all polynomials.  
m.polynomials <- matrix(0, n.observations, polynomial.order)
for(i in 1:polynomial.order){
	m.polynomials[,i] <- x.at.obs^i
}


polynomial.model <- lm(observations ~ m.polynomials)
polynomial.predictions <- fitted(polynomial.model)

plot(x.at.obs, observations, 
     xlim = c(-2, 2), ylim = c(-1, 2), xaxt = "n", yaxt = "n", 
     main = "", xlab = "x", ylab = "y", cex.lab = 1.5)
lines(x.at.obs, polynomial.predictions, col = 153, lwd = 2)
points(x.at.obs, observations)
    
#  Clean up axes.
axis.size <- 1.4
x.axis <- seq(from = -2, to = 2, by = 1)
y.axis <- seq(from = -1, to = 2, by = 0.5)
axis(1, at = x.axis, labels = x.axis, cex.axis = axis.size)
axis(2, at = y.axis, labels = y.axis, cex.axis = axis.size)

#  Close the graphics device.  
dev.print(device = postscript, "15.2.added.noise", horizontal = TRUE)