# October 4, 2005 # Repeat Clive's matlab demo polynomial regression yr <- seq(1900,1980,10) pop<- c(76.2,92.2,106.0,123.2,132.2,151.3,179.3,203.3,226.5) myplt <- function(extra=F,xl=c(1900,2000),yl=c(0,300)) { plot(yr,pop,pch=16,xlim=xl,ylim=yl,ylab="u.s. pop") if(extra) points( c(1990,2000), c(248.7,281.4), pch=16, col=2 ) } # fit polynomial model uspop <- function() { if(length(dev.list())==1) { my.ps.color() } x <- seq(0,.8,.1); xx <- seq(0, 1,l=351) X <- rep(1,9); XX <- rep(1,length(xx)) for(k in 1:8) { X <- cbind(X,x^k); XX <- cbind(XX,xx^k) ab <- lsfit(X,pop,int=F)$coef; yh <- c( XX %*% ab ) myplt(T); lines(100*xx+1900,yh,lwd=2,col=3) title(paste("polynomial order ",k)) print(" use ps.copy and mv to poly1.eps ") browser() } } uspop()