verification of PRML fig_1.5

lm_test.R

sr <- function() source('lm_test.R')

make_data <- function(n,sd=0.4) {
	x <- seq(0,2*pi,length=n)
	y <- sin(x) + rnorm(n,0,sd)
	return(data.frame(x=x,y=y))
}

my_test <- function(d,m) {
	z <- lm(y~poly(x,m),data=d)
	return(z)
}

main <- function(d_learn=10,d_test=100) {
	d <- make_data(d_learn)
	layout(t(matrix(1:9,3,3)))
	nd <- make_data(d_test)
	rl <- c()
	rt <- c()
	for(i in 2:9) {
		z <- my_test(d,i)
		r1 <- sqrt(sum(residuals(z)^2)/d_learn)
		aic1 <- AIC(z)
		rl <- cbind(rl,c(i,r1))

		y <- predict(z,nd)
		r2 <- sqrt(sum((y-nd$y)^2)/d_test)
		rt <- cbind(rt,c(i,r2))

		s<-sprintf('m=%d,r1=%f,aic1=%f,r2=%f',i,r1,aic1,r2)
		print(s)
		s<-sprintf('m=%d,r1=%.2f,r2=%.2f',i,r1,r2)
		plot(d$x,d$y,pch=19,main=s)
		lines(nd$x,predict(z,nd))
	}
	plot(rl[1,],rl[2,],ylim=c(0,1),col=1)
	lines(rl[1,],rl[2,],ylim=c(0,1),col=1)
	par(new=T)
	plot(rt[1,],rt[2,],ylim=c(0,1),col=6)
	lines(rt[1,],rt[2,],ylim=c(0,1),col=6)
}

main()