takousiki2.R

# original in http://d.hatena.ne.jp/syou6162/20080801/1217528370

url <- "http://research.microsoft.com/en-us/um/people/cmbishop/PRML/webdatasets/curvefitting.txt"
curvefitting <- read.csv(file=url, sep=" ",head=F)

colnames(curvefitting)<-c("x","t")



takousiki <- function(n){
  plot(function(x){sin(2*pi*x)},0,1,ylab="",ylim=c(-1.5,1.5))
  par(new=T)
  ##### plot raw data (before curvefitting) 
  # x <- seq(0,1,by=0.1)
  # y <- sin(2 * pi * seq(0,1,by=0.1)) + rnorm(11,sd=0.5)
  # plot(x,y,ylab="",ylim=c(-1.5,1.5))
  #
  plot(curvefitting,ylim=c(-1.5,1.5),xlim=c(0,1),axes=F,col="blue")
  x <- curvefitting$x
  y <- curvefitting$t
  par(new=T)
  ##### (xn <<- x^n) in 1:n ###########################################################################################
  # > x
  #  [1] 0.000000 0.111111 0.222222 0.333333 0.444444 0.555556 0.666667 0.777778 0.888889 1.000000
  # > x^2
  #  [1] 0.00000000 0.01234565 0.04938262 0.11111089 0.19753047 0.30864247 0.44444489 0.60493862 0.79012365 1.00000000
  #
  # (function(s){s^3})(seq(2)) => 1 8
  # (function(s){paste("x",s," <<- ","x^",s,sep="")})(seq(2)) => "x1 <<- x^1" "x2 <<- x^2"
  #
  invisible(mapply(function(v){eval(parse(text=v))},(function(x){paste("x",x," <<- ","x^",x,sep="")})(seq(n))))
  ##### sigma(x1+x2+...+xn) in 1:n ####################################################################################
  # ?Reduce => Higher-Order Function
  #
  # > lm(y ~ x1 + x2)$coefficients
  # (Intercept)           x1           x2  
  #      0.9143      -1.9028       0.6350  
  #
  # eval(parse(text=paste("lm <- lm(y ~ ",Reduce(function(x,init){paste(x,init,sep=" + ")},paste("x",seq(n),sep="")),")",sep="")))
  text.xsum = Reduce(function(x,init){paste(x,init,sep=" + ")},paste("x",seq(n),sep=""))
  eval(parse(text=paste("lm <- lm(y ~ ",text.xsum,")",sep="")))
  ##### (xn <<- newdata^1) in 1:n #####################################################################################
  newdata <- seq(0,1,by=0.01)
  text.ndf <- Reduce(function(x,init){paste(x,init,sep=",")},(function(x){paste("x",x,"=newdata^",x,sep="")})(seq(n)))
  eval(parse(text=paste("newdata <- data.frame(",text.ndf,")",sep="")))
  par(new=T)
  ##### plot fitting curve #################################################################
  plot(seq(0,1,by=0.01),predict(lm,newdata),type="l",col="red",ann=F,axes=F)
}