tgp2.R
In tgp: Bayesian Treed Gaussian Process Models

### R code from vignette source 'tgp2.Rnw'

###################################################
### code chunk number 1: tgp2.Rnw:33-35
###################################################
library(tgp)
options(width=65)


###################################################
### code chunk number 2: tgp2.Rnw:148-149
###################################################
seed <- 1; set.seed(seed)  ## seed zero problematic with null btlm map tree below


###################################################
### code chunk number 3: tgp2.Rnw:214-218
###################################################
fb.train <- fried.bool(500)
X <- fb.train[,1:13]; Z <- fb.train$Y
fb.test <- fried.bool(1000)
XX <- fb.test[,1:13]; ZZ <- fb.test$Ytrue


###################################################
### code chunk number 4: tgp2.Rnw:224-225
###################################################
names(X)


###################################################
### code chunk number 5: tgp2.Rnw:232-235
###################################################
fit1 <- bcart(X=X, Z=Z, XX=XX, verb=0)
rmse1 <- sqrt(mean((fit1$ZZ.mean - ZZ)^2))
rmse1


###################################################
### code chunk number 6: cat-fbcart-mapt
###################################################
tgp.trees(fit1, "map")


###################################################
### code chunk number 7: tgp2.Rnw:244-245
###################################################
graphics.off()


###################################################
### code chunk number 8: tgp2.Rnw:259-262
###################################################
fit2 <- btlm(X=X, Z=Z, XX=XX, verb=0)
rmse2 <- sqrt(mean((fit2$ZZ.mean - ZZ)^2))
rmse2


###################################################
### code chunk number 9: cat-fbtlm-mapt
###################################################
tgp.trees(fit2, "map")


###################################################
### code chunk number 10: tgp2.Rnw:270-271
###################################################
graphics.off()


###################################################
### code chunk number 11: tgp2.Rnw:305-308
###################################################
fit3 <- btlm(X=X, Z=Z, XX=XX, basemax=10, verb=0)
rmse3 <- sqrt(mean((fit3$ZZ.mean - ZZ)^2))
rmse3


###################################################
### code chunk number 12: cat-fbtlm-mapt
###################################################
tgp.trees(fit3, "map")


###################################################
### code chunk number 13: tgp2.Rnw:314-315
###################################################
graphics.off()


###################################################
### code chunk number 14: tgp2.Rnw:337-340
###################################################
fit4 <- btgpllm(X=X, Z=Z, XX=XX, verb=0)
rmse4 <- sqrt(mean((fit4$ZZ.mean - ZZ)^2))
rmse4


###################################################
### code chunk number 15: tgp2.Rnw:345-346
###################################################
fit4$gpcs


###################################################
### code chunk number 16: tgp2.Rnw:357-360
###################################################
fit5 <-  btgpllm(X=X, Z=Z, XX=XX, basemax=10, verb=0)
rmse5 <- sqrt(mean((fit5$ZZ.mean - ZZ)^2))
rmse5 


###################################################
### code chunk number 17: cat-fb-mapt
###################################################
h <- fit1$post$height[which.max(fit1$posts$lpost)]
tgp.trees(fit5, "map")


###################################################
### code chunk number 18: tgp2.Rnw:377-378
###################################################
graphics.off()


###################################################
### code chunk number 19: tgp2.Rnw:409-412
###################################################
fit6 <-  btgpllm(X=X, Z=Z, XX=XX, basemax=10, splitmin=11, verb=0)
rmse6 <- sqrt(mean((fit6$ZZ.mean - ZZ)^2))
rmse6


###################################################
### code chunk number 20: tgp2.Rnw:448-449
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 21: tgp2.Rnw:778-779
###################################################
f <- friedman.1.data(250) 


###################################################
### code chunk number 22: tgp2.Rnw:785-788
###################################################
Xf <- f[, 1:6] 
Zf <- f$Y 
sf <- sens(X=Xf, Z=Zf, nn.lhs=600, model=bgpllm, verb=0)


###################################################
### code chunk number 23: tgp2.Rnw:799-800
###################################################
names(sf$sens)


###################################################
### code chunk number 24: sens-full
###################################################
plot(sf, layout="sens", legendloc="topleft")


###################################################
### code chunk number 25: tgp2.Rnw:818-819
###################################################
graphics.off()


###################################################
### code chunk number 26: sens-mains
###################################################
par(mar=c(4,2,4,2), mfrow=c(2,3))
plot(sf, layout="sens", maineff=t(1:6))


###################################################
### code chunk number 27: tgp2.Rnw:846-847
###################################################
graphics.off()


###################################################
### code chunk number 28: sens-indices
###################################################
plot(sf, layout="sens", maineff=FALSE)


###################################################
### code chunk number 29: tgp2.Rnw:859-860
###################################################
graphics.off()


###################################################
### code chunk number 30: tgp2.Rnw:910-915
###################################################
X <- airquality[,2:4]
Z <- airquality$Ozone
rect <- t(apply(X, 2, range, na.rm=TRUE))
mode <- apply(X , 2, mean, na.rm=TRUE)
shape <- rep(2,3)


###################################################
### code chunk number 31: sens-udraw
###################################################
Udraw <- lhs(300, rect=rect, mode=mode, shape=shape)
par(mfrow=c(1,3), mar=c(4,2,4,2))
for(i in 1:3){
  hist(Udraw[,i], breaks=10,xlab=names(X)[i], 
       main="",ylab="", border=grey(.9), col=8) 
}  


###################################################
### code chunk number 32: tgp2.Rnw:928-929
###################################################
graphics.off()


###################################################
### code chunk number 33: tgp2.Rnw:941-943
###################################################
s.air <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, rect=rect, 
                               shape=shape, mode=mode, verb=0))


###################################################
### code chunk number 34: sens-air1
###################################################
plot(s.air, layout="sens")


###################################################
### code chunk number 35: tgp2.Rnw:950-951
###################################################
graphics.off()


###################################################
### code chunk number 36: tgp2.Rnw:967-970
###################################################
rect[2,] <- c(0,5)
mode[2] <- 2
shape[2] <- 2


###################################################
### code chunk number 37: tgp2.Rnw:974-975
###################################################
sens.p <- suppressWarnings(sens(X=X,Z=Z,nn.lhs=300, model=NULL, rect=rect, shape=shape, mode=mode))


###################################################
### code chunk number 38: sens-air2
###################################################
s.air2 <- predict(s.air, BTE=c(1,1000,1), sens.p=sens.p, verb=0) 
plot(s.air2, layout="sens")


###################################################
### code chunk number 39: tgp2.Rnw:982-983
###################################################
graphics.off()


###################################################
### code chunk number 40: tgp2.Rnw:1006-1013
###################################################
X$Temp[X$Temp >70] <- 1
X$Temp[X$Temp >1] <- 0
rect <- t(apply(X, 2, range, na.rm=TRUE))
mode <- apply(X , 2, mean, na.rm=TRUE)
shape <- c(2,2,0)
s.air <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, rect=rect, 
                               shape=shape, mode=mode, verb=0, basemax=2))


###################################################
### code chunk number 41: sens-air3
###################################################
plot(s.air, layout="sens")


###################################################
### code chunk number 42: tgp2.Rnw:1019-1020
###################################################
graphics.off()


###################################################
### code chunk number 43: tgp2.Rnw:1033-1034
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 44: tgp2.Rnw:1207-1211
###################################################
rosenbrock <- function(x){ 
  x <- matrix(x, ncol=2)
  100*(x[,1]^2 - x[,2])^2 + (x[,1] - 1)^2 
}


###################################################
### code chunk number 45: tgp2.Rnw:1216-1217
###################################################
rosenbrock(c(1,1))


###################################################
### code chunk number 46: tgp2.Rnw:1225-1228
###################################################
rect <- cbind(c(-1,-1),c(5,5))
X <- lhs(40, rect)
Z <- rosenbrock(X)


###################################################
### code chunk number 47: tgp2.Rnw:1244-1246
###################################################
XX <- lhs(200, rect)
rfit <- bgp(X,Z,XX,improv=c(1,10), verb=0)


###################################################
### code chunk number 48: tgp2.Rnw:1254-1255
###################################################
cbind(rfit$improv,XX)[rfit$improv$rank <= 10,]


###################################################
### code chunk number 49: optim-fit1
###################################################
plot(rfit, as="improv")


###################################################
### code chunk number 50: tgp2.Rnw:1269-1270
###################################################
graphics.off()


###################################################
### code chunk number 51: optim-fit2
###################################################
rfit2 <- predict(rfit, XX=XX, BTE=c(1,1000,1), improv=c(5,20), verb=0) 
plot(rfit2, layout="as", as="improv")


###################################################
### code chunk number 52: tgp2.Rnw:1298-1299
###################################################
graphics.off()


###################################################
### code chunk number 53: tgp2.Rnw:1438-1439
###################################################
f <- function(x) { exp2d.Z(x)$Z }


###################################################
### code chunk number 54: tgp2.Rnw:1453-1456
###################################################
rect <- rbind(c(-2,6), c(-2,6))
X <- lhs(20, rect)
Z <- f(X)


###################################################
### code chunk number 55: tgp2.Rnw:1460-1473
###################################################
out <- progress <- NULL
for(i in 1:20) {
  
  ## get recommendations for the next point to sample
  out <- optim.step.tgp(f, X=X, Z=Z, rect=rect, prev=out, verb=0)

  ## add in the inputs, and newly sampled outputs
  X <- rbind(X, out$X)
  Z <- c(Z, f(out$X))
  
  ## keep track of progress and best optimum
  progress <- rbind(progress, out$progress)
}


###################################################
### code chunk number 56: optim-progress
###################################################
par(mfrow=c(1,2))
matplot(progress[,1:2], main="x progress",
        xlab="rounds", ylab="x[,1:2]", type="l", lwd=2)
legend("topright", c("x1", "x2"), lwd=2, col=1:2, lty=1:2)
plot(log(progress$improv), type="l", main="max log improv",
     xlab="rounds", ylab="max log(improv)")


###################################################
### code chunk number 57: tgp2.Rnw:1490-1491
###################################################
graphics.off()


###################################################
### code chunk number 58: tgp2.Rnw:1506-1507
###################################################
out$progress[1:2]


###################################################
### code chunk number 59: tgp2.Rnw:1532-1533
###################################################
formals(optim)$method


###################################################
### code chunk number 60: tgp2.Rnw:1537-1538
###################################################
formals(optim.ptgpf)$method


###################################################
### code chunk number 61: tgp2.Rnw:1552-1553
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 62: tgp2.Rnw:1670-1673
###################################################
geo <- default.itemps(type="geometric")
har <- default.itemps(type="harmonic")
sig <- default.itemps(type="sigmoidal")


###################################################
### code chunk number 63: it-itemps
###################################################
par(mfrow=c(2,1))
all <- cbind(geo$k, har$k, sig$k)
matplot(all, pch=21:23,
        main="inv-temp ladders", xlab="indx", ylab="itemp")
legend("topright", pch=21:23, 
       c("geometric","harmonic","sigmoidal"), col=1:3)
matplot(log(all), pch=21:23,
        main="log(inv-temp) ladders", xlab="indx", ylab="itemp")


###################################################
### code chunk number 64: tgp2.Rnw:1688-1689
###################################################
graphics.off()


###################################################
### code chunk number 65: tgp2.Rnw:1755-1762
###################################################
ESS <- function(w)
{
  mw <- mean(w)
  cv2 <- sum((w-mw)^2)/((length(w)-1)*mw^2)
  ess <- length(w)/(1+cv2)
  return(ess)
}


###################################################
### code chunk number 66: tgp2.Rnw:1908-1911
###################################################
exp2d.data<-exp2d.rand() 
X<-exp2d.data$X 
Z<-exp2d.data$Z 


###################################################
### code chunk number 67: tgp2.Rnw:1917-1920
###################################################
its <- default.itemps(m=10)
exp.btlm <- btlm(X=X,Z=Z, bprior="b0", R=2, itemps=its, pred.n=FALSE,
                 BTE=c(1000,3000,2)) 


###################################################
### code chunk number 68: tgp2.Rnw:1945-1946
###################################################
exp.btlm$ess


###################################################
### code chunk number 69: tgp2.Rnw:1957-1960
###################################################
library(MASS)
moto.it <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
        bprior="b0", R=3, itemps=geo, trace=TRUE, pred.n=FALSE, verb=0)


###################################################
### code chunk number 70: tgp2.Rnw:1964-1965
###################################################
moto.it$ess$combined


###################################################
### code chunk number 71: tgp2.Rnw:1969-1971
###################################################
p <- moto.it$trace$post
ESS(p$wlambda)


###################################################
### code chunk number 72: tgp2.Rnw:1977-1978
###################################################
ESS(p$w)


###################################################
### code chunk number 73: tgp2.Rnw:1983-1984
###################################################
as.numeric(c(sum(p$itemp == 1), moto.it$ess$each[1,2:3]))


###################################################
### code chunk number 74: tgp2.Rnw:1995-1997
###################################################
moto.reg <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
        R=3, bprior="b0", trace=TRUE, pred.n=FALSE, verb=0)


###################################################
### code chunk number 75: tgp2.Rnw:2003-2006
###################################################
L <- length(p$height)
hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE))
b <- hist2bar(cbind(moto.reg$trace$post$height, p$height, hw))


###################################################
### code chunk number 76: it-moto-height
###################################################
barplot(b, beside=TRUE, col=1:3, xlab="tree height", ylab="counts", 
         main="tree heights encountered")
legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3)


###################################################
### code chunk number 77: tgp2.Rnw:2014-2015
###################################################
graphics.off()


###################################################
### code chunk number 78: it-moto-ktrace
###################################################
plot(log(moto.it$trace$post$itemp), type="l", ylab="log(k)", xlab="samples",
     main="trace of log(k)")


###################################################
### code chunk number 79: tgp2.Rnw:2048-2049
###################################################
graphics.off()


###################################################
### code chunk number 80: it-moto-khist
###################################################
b <- itemps.barplot(moto.it, plot.it=FALSE)
barplot(t(cbind(moto.it$itemps$counts, b)), col=1:2,
        beside=TRUE, ylab="counts", xlab="itemps", 
        main="inv-temp observation counts")
legend("topleft", c("observation counts", "posterior samples"), fill=1:2)


###################################################
### code chunk number 81: tgp2.Rnw:2080-2081
###################################################
graphics.off()


###################################################
### code chunk number 82: tgp2.Rnw:2104-2106
###################################################
moto.it.sig <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
                      R=3, bprior="b0", krige=FALSE, itemps=sig, verb=0)


###################################################
### code chunk number 83: tgp2.Rnw:2110-2111
###################################################
moto.it.sig$ess$combined


###################################################
### code chunk number 84: it-moto-pred
###################################################
plot(moto.it.sig)


###################################################
### code chunk number 85: tgp2.Rnw:2117-2118
###################################################
graphics.off()


###################################################
### code chunk number 86: tgp2.Rnw:2144-2147
###################################################
Xcand <- lhs(10000, rbind(c(-6,6),c(-6,6)))
X <- dopt.gp(400, X=NULL, Xcand)$XX
Z <- exp2d.Z(X)$Z


###################################################
### code chunk number 87: tgp2.Rnw:2152-2154
###################################################
exp.reg <- btgpllm(X=X, Z=Z, BTE=c(2000,52000,10), bprior="b0", 
                   trace=TRUE, krige=FALSE, R=10, verb=0)


###################################################
### code chunk number 88: it-exp-pred
###################################################
plot(exp.reg)


###################################################
### code chunk number 89: tgp2.Rnw:2160-2161
###################################################
graphics.off()


###################################################
### code chunk number 90: tgp2.Rnw:2173-2175
###################################################
h <- exp.reg$post$height[which.max(exp.reg$posts$lpost)]
h


###################################################
### code chunk number 91: it-exp-mapt
###################################################
tgp.trees(exp.reg, "map")


###################################################
### code chunk number 92: tgp2.Rnw:2184-2185
###################################################
graphics.off()


###################################################
### code chunk number 93: tgp2.Rnw:2209-2212
###################################################
its <- default.itemps(k.min=0.02)
exp.it <- btgpllm(X=X, Z=Z, BTE=c(2000,52000,10), bprior="b0", 
               trace=TRUE, krige=FALSE, itemps=its, R=10, verb=0)


###################################################
### code chunk number 94: tgp2.Rnw:2217-2219
###################################################
exp.it$gpcs
exp.reg$gpcs


###################################################
### code chunk number 95: tgp2.Rnw:2227-2229
###################################################
p <- exp.it$trace$post
data.frame(ST=sum(p$itemp == 1), nIT=ESS(p$w), oIT=exp.it$ess$combined)


###################################################
### code chunk number 96: tgp2.Rnw:2241-2244
###################################################
L <- length(p$height)
hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE))
b <- hist2bar(cbind(exp.reg$trace$post$height, p$height, hw))


###################################################
### code chunk number 97: it-exp-height
###################################################
barplot(b, beside=TRUE, col=1:3, xlab="tree height", ylab="counts", 
         main="tree heights encountered")
legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3)


###################################################
### code chunk number 98: tgp2.Rnw:2252-2253
###################################################
graphics.off()


###################################################
### code chunk number 99: it-exp-trace-height
###################################################
ylim <- range(p$height, exp.reg$trace$post$height)
plot(p$height, type="l", main="trace of tree heights", 
     xlab="t", ylab="height", ylim=ylim)
lines(exp.reg$trace$post$height, col=2)
legend("topright", c("tempered", "reg MCMC"), lty=c(1,1), col=1:2)


###################################################
### code chunk number 100: tgp2.Rnw:2277-2278
###################################################
graphics.off()


###################################################
### code chunk number 101: it-expit-pred
###################################################
plot(exp.it)


###################################################
### code chunk number 102: it-expit-trees
###################################################
tgp.trees(exp.it, "map")


###################################################
### code chunk number 103: tgp2.Rnw:2305-2306
###################################################
graphics.off()