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.001
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_tr04.out")

print(paste("cx.parm =", cx.parm))    
print(paste("mse.lrn =", mse.lrn))
print(paste("mse.val =", mse.val))

source("psopts.r");
postscript(file="bv_tr04_tree.eps");
par.default <- par(no.readonly=TRUE)
par(mar=c(0,0,0,0))
plot(tree)
text(tree) 
par(par.default)
dev.off()

postscript(file="bv_tr04_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()

postscript(file="bv_tr04_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()

postscript(file="bv_tr04_val.eps")
plot(x=lrn$x,y=lrn$y,type="n",xlab="feature",ylab="target")
points(x=lrn$x,y=val$y,pch="o",cex=0.5)
lines(x=lrn$x,y=yhat,type="l",lty="solid",col="red",lwd=1)
#lines(x=tru$x,y=tru$f,type="l",lty="solid",col="blue",lwd=2)
dev.off()

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_tr04_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()