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

data(lidar)
attach(lidar)
x = range
y = logratio
n = length(x)
namex = "Distance travelled before the light is reflected back to its source"
namey = "Logarithm of the ratio of received light from two laser sources"


par(mfrow=c(1,1))
plot(x,y,xlab=namex,ylab=namey,main="")

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

mae = rbind(
apply(abs(error[x<600,]),2,median),
apply(abs(error[x>=600,]),2,median),
apply(abs(error),2,median))
rmae = round(100*(mae/mae[,1]-1),1)[,2:5]

model = c("Spline","BNN","CART","Random forest","BART")


pdf(file="brnn-1.pdf",width=15,height=10)
par(mfrow=c(1,1))
plot(x,y,xlab=namex,ylab=namey,main="")
lines(x,pred.spline,col=1,lwd=2)
lines(x,pred.nn,col=2,lwd=2) 
lines(x,pred.tree,col=3,lwd=2)
lines(x,pred.rf,col=4,lwd=2) 
lines(x,pred.bart,col=6,lwd=2)
abline(v=600,lty=2)
legend("bottomleft",legend=model,col=c(1:4,6),lwd=2,bty="n")

par(mfrow=c(2,3))
for (i in 1:5){
  plot(x,error[,i],xlab=namex,ylab="Error",ylim=range(error),main=model[i])
  abline(v=600,lty=2)
  if (i==1){
    text(500,0.4,paste("MAE = ",round(mae[1,1],4),sep=""))
    text(650,0.4,paste("MAE = ",round(mae[2,1],4),sep=""))
  }else{
    text(500,0.4,paste("RMAE = ",rmae[1,i-1],"%",sep=""))
    text(650,0.4,paste("RMAE = ",rmae[2,i-1],"%",sep=""))
  	
  }
}
boxplot(error,names=model,ylim=c(min(error),1.1*max(error)))
for (i in 1:5)
  text(i,1.05*max(error),round(mae[3,i],3))
dev.off()