R : Copyright 2002, The R Development Core Team
Version 1.5.1  (2002-06-17)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> invisible(options(echo = TRUE))
> library(rpart)
> 
> lrn <- read.table("bv_lrn.dat",
+          header=F,colClasses="numeric",col.names=c("x","y"))
> val <- read.table("bv_val.dat",
+          header=F,colClasses="numeric",col.names=c("x","y"))
> tru <- read.table("bv_tru.dat",
+          header=F,colClasses="numeric",col.names=c("x","f"))
> 
> cx.parm <- 0.00075
> tree <- rpart(y ~ x, data=lrn, method="anova",cp=cx.parm,minbucket=1);
> yhat <- tree$frame[tree$where,5] 
> 
> mse.lrn <- mean((lrn$y-yhat)^2)
> mse.val <- mean((val$y-yhat)^2)
> 
> #print(tree$frame)
> #summary.rpart(tree,file="bv_tr05.out")
> 
> print(paste("cx.parm =", cx.parm))
[1] "cx.parm = 0.00075"
> print(paste("mse.lrn =", mse.lrn))
[1] "mse.lrn = 0.749993828236913"
> print(paste("mse.val =", mse.val))
[1] "mse.val = 1.33091181560658"
> 
> source("psopts.r");
> postscript(file="bv_tr05_tree.eps");
> par.default <- par(no.readonly=TRUE)
> par(mar=c(0,0,0,0))
> plot(tree)
> text(tree) 
> par(par.default)
> dev.off()
null device 
          1 
> 
> postscript(file="bv_tr05_lrn.eps")
> plot(x=lrn$x,y=lrn$y,type="n",xlab="feature",ylab="target")
> points(x=lrn$x,y=lrn$y,pch="o",cex=0.5)
> lines(x=lrn$x,y=yhat,type="l",lty="solid",col="red",lwd=2)
> dev.off()
null device 
          1 
> 
> postscript(file="bv_tr05_tru.eps")
> plot(x=lrn$x,y=lrn$y,type="n",xlab="feature",ylab="target")
> points(x=lrn$x,y=lrn$y,pch="o",cex=0.5)
> lines(x=lrn$x,y=yhat,type="l",lty="solid",col="red",lwd=2)
> lines(x=tru$x,y=tru$f,type="l",lty="solid",col="blue",lwd=2)
> dev.off()
null device 
          1 
> 
> 
> rm(val,tru)
> set.seed(100542)
> nobs <- 500
> sdv <- 1.0
> loop <-100
> 
> x <- seq(-1,1,length=nobs)
> f <- 6*x^3 + 1.5*sin(11*x)*exp(-x^2*6)
> f <- 1.4*f
> 
> yave <- yhat
> postscript(file="bv_tr05_ave.eps")
> plot(x=lrn$x,y=lrn$y,type="n",xlab="feature",ylab="target")
> for (i in 1:loop) {
+   e <- rnorm(n=nobs,mean=0.0,sd=sdv)
+   y <- f + e
+   lrn$y <- y
+   tree <- rpart(y ~ x, data=lrn, method="anova",cp=cx.parm,minbucket=1);
+   yhat <- tree$frame[tree$where,5]   
+   if (i<26) lines(x=lrn$x,y=yhat,type="l",lty="solid",col="black",lwd=0.25)
+   yave <- yave + yhat
+ }
> yave <- yave/loop
> lines(x=lrn$x,y=f,type="l",lty="solid",col="blue",lwd=2)
> lines(x=lrn$x,y=yave,type="l",lty="solid",col="red",lwd=2)
> dev.off()
null device 
          1 
> 
> proc.time()
[1] 12.07  0.20 12.38  0.00  0.00
>