## load needed libraries
library(BART)
library(nonlinvarsel)
library(foreach)


## read in cars data, rescale price and mileage, make sure categorical variables are factors
cdat = read.csv("http://www.rob-mcculloch.org/data/susedcars.csv")
cdat$price = cdat$price/1000
cdat$mileage = cdat$mileage/1000
facv = c(2,3,6,7)
for(i in 1:length(facv)) cdat[[facv[i]]] = as.factor(cdat[[facv[i]]])


## pull off y and x
y = cdat$price
# dummy up categorical variables
x = bartModMat(cdat[,-1])

## optional train/test split
set.seed(99)
n = length(y)
ii = sample(1:n,floor(.75*n))
xp = x[-ii,] ; yp = y[-ii]
x = x[ii,]; y = y[ii]

## run bart
set.seed(14)
bf = wbart(x,y,xp,nskip=500,ndpost = 3000) #BART fit

## forward variable selection
vsfr = vsf(xp,bf$yhat.test.mean)
plot(vsfr,cex.axis=.6)
print(vsfr)

## assess uncertainty
sp = sumpost(bf$yhat.test,vsfr$fit,bf$yhat.test.mean,distrmse)
plot(sp,diff=FALSE)
plot(sp)