install.packages("SemiPar")
install.packages("brnn")
library("MASS")
library("SemiPar")
library("brnn")
library("tree")
library("BART")
library("randomForest")

# Data set 1
x1 = seq(0,0.23,length.out=25)
y1 = 4*x1+rnorm(25,sd=0.1)
x2 = seq(0.25,0.75,length.out=50) 
y2 = 2-4*x2+rnorm(50,sd=0.1) 
x3 = seq(0.77,1,length.out=25) 
y3 = 4*x3-4+rnorm(25,sd=0.1) 
x1  = c(x1,x2,x3) 
y1  = c(y1,y2,y3) 

# Data set 2
n = nrow(mcycle)
x2 = mcycle$times
y2 = mcycle$accel

# Data set 3
data(lidar)
attach(lidar)
x3 = range
y3 = logratio



par(mfrow=c(2,2))
plot(x1,y1)
plot(x2,y2)
plot(x3,y3)


x = x1
y = y1
n = length(x)
namex = "X"
namey = "Y"

x = x2
y = y2
n = length(x)
x = x + rnorm(n,0,0.01)
namex = "Milliseconds after impact"
namey = "Acceleration in g"

x = x3
y = y3
n = length(x)
namex = "distance travelled"
namey = "Log of the ratio of received light"


error = matrix(0,n,5)

# Smoothed-spline 
spline.fit = smooth.spline(x,y,all.knots=TRUE)
pred.spline = spline.fit$y
error[,1] = y-pred.spline

# CART 
tree.fit = tree(y~x)
pred.tree  = predict(tree.fit)
error[,2] = y-pred.tree

# Random forest
rf         = randomForest(y~x) 
pred.rf    = rf$predicted
error[,3] = y-pred.rf

# BART for classification
bart       = wbart(x,y)
pred.bart  = apply(bart$yhat.train,2,mean)
error[,4] = y-pred.bart

# Bayesian NN
nn.fit = brnn(y~x,neurons=4)
pred.nn = predict(nn.fit)
error[,5] = y-pred.nn

par(mfrow=c(3,3),mai=c(0.8,0.8,0.2,0.2))
plot(x,y,xlab=namex,ylab=namey,main="")
lines(x,pred.spline,col=2,lwd=2)
lines(x,pred.nn,col=4,lwd=2) 
legend("bottomleft",legend=c("Smooth spline","BNN"),col=c(2,4),lwd=2,bty="n")

plot(x,y,xlab=namex,ylab=namey,main="")
lines(x,pred.tree,col=2,lwd=2)
lines(x,pred.rf,col=3,lwd=2) 
lines(x,pred.bart,col=4,lwd=2)
legend("bottomleft",legend=c("CART","Random forest","BART"),col=2:4,lwd=2,bty="n")

boxplot(error,names=c("Spline","CART","RF","BNN","BART"),ylim=c(min(error),1.1*max(error)))
for (i in 1:5)
  text(i,1.05*max(error),round(median(abs(error[,i])),4))