library(PolynomF)

f <- function(x) sin(2 * pi * x)
plot(f, 0, 1)

x1 <- seq(from=0, to=1, length=10)
d1 <- f(x1)
d1 <- d1 + rnorm(length(d1), sd=0.2)
training.data <- data.frame(x=x1, d=d1)

x2 <- seq(from=0, to=1, length=100)
d2 <- f(x2)
d2 <- d2 + rnorm(length(d2), sd=0.2)
test.data <- data.frame(x=x2, d=d2)

RMS.training <- RMS.test <- numeric()
tmp <- list()

#	pdf(file="overfitting.pdf")
cairo_pdf(file="overfitting.pdf", onefile=TRUE)

old <-par(no.readonly=TRUE)
par(mfrow=c(3, 3))

for (i in 1:9) {

	m <- lm (d ~ poly(x, degree=i, raw=TRUE), data=training.data)

	plot(f, ylim=c(-1.5, 1.5), col="green")
	points(training.data, col="blue")
	curve(polynom(coef(m))(x), col="red", add=TRUE)
	title(main=substitute(M == x, list(x = i)))

	RMS.training[i] <- sqrt(sum((predict(m) - training.data$d) ^ 2) / 10)
	RMS.test[i] <- sqrt(sum((predict(m, test.data[1]) - test.data$d) ^ 2) / 100)
	tmp[[i]] <- m
}

par(old)

plot(NULL,
	xlim=c(1, 9),
	ylim=range(RMS.training, RMS.test),
	xlab=expression(M),
	ylab=expression(E[RMS])
)
points(RMS.training, col="blue", pch=20, type="o")
points(RMS.test, col="red", pch=20, type="o")
legend(
	"bottomleft",
	c("Training", "Test"),
	col=c("blue", "red"),
	pch=20,
	lty="solid"
)

dev.off()

lapply(tmp, function(z) { z <- coef(z); attributes(z) <- NULL ; z})
