tgp.R
In tgp: Bayesian Treed Gaussian Process Models

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

###################################################
### code chunk number 1: tgp.Rnw:26-28
###################################################
library(tgp)
options(width=65)


###################################################
### code chunk number 2: tgp.Rnw:185-186
###################################################
bgp


###################################################
### code chunk number 3: gpllm
###################################################
hist(c(rgamma(100000,1,20), rgamma(100000,10,10)), 
     breaks=50, xlim=c(0,2), freq=FALSE, ylim=c(0,3),
     main = "p(d) = G(1,20) + G(10,10)", xlab="d")
d <- seq(0,2,length=1000)
lines(d,0.2+0.7/(1+exp(-10*(d-0.5))))
abline(h=1, lty=2)
legend(x=1.25, y=2.5, c("p(b) = 1", "p(b|d)"), lty=c(1,2))


###################################################
### code chunk number 4: tgp.Rnw:668-669
###################################################
graphics.off()


###################################################
### code chunk number 5: tgp.Rnw:967-968
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 6: tgp.Rnw:984-988
###################################################
# 1-d linear data input and predictive data
X <- seq(0,1,length=50)  # inputs
XX <- seq(0,1,length=99) # predictive locations
Z <- 1 + 2*X + rnorm(length(X),sd=0.25) # responses


###################################################
### code chunk number 7: tgp.Rnw:993-994
###################################################
lin.blm <- blm(X=X, XX=XX, Z=Z)


###################################################
### code chunk number 8: linear-blm
###################################################
plot(lin.blm, main='Linear Model,', layout='surf')
abline(1,2,lty=3,col='blue')


###################################################
### code chunk number 9: tgp.Rnw:1002-1003
###################################################
graphics.off()


###################################################
### code chunk number 10: tgp.Rnw:1036-1037
###################################################
lin.gpllm <- bgpllm(X=X, XX=XX, Z=Z)


###################################################
### code chunk number 11: linear-gplm
###################################################
plot(lin.gpllm, main='GP LLM,', layout='surf')
abline(1,2,lty=4,col='blue')


###################################################
### code chunk number 12: tgp.Rnw:1045-1046
###################################################
graphics.off()


###################################################
### code chunk number 13: tgp.Rnw:1066-1070
###################################################
lin.gpllm.tr <- bgpllm(X=X, XX=0.5, Z=Z, pred.n=FALSE, trace=TRUE,
                       verb=0)
mla <- mean(lin.gpllm.tr$trace$linarea$la)
mla


###################################################
### code chunk number 14: tgp.Rnw:1076-1077
###################################################
1-mean(lin.gpllm.tr$trace$XX[[1]]$b1)


###################################################
### code chunk number 15: tgp.Rnw:1085-1086
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 16: tgp.Rnw:1104-1110
###################################################
X <- seq(0,20,length=100)
XX <- seq(0,20,length=99)
Ztrue <- (sin(pi*X/5) + 0.2*cos(4*pi*X/5)) * (X <= 9.6)
lin <- X>9.6; 
Ztrue[lin] <- -1 + X[lin]/10
Z <- Ztrue + rnorm(length(Ztrue), sd=0.1)


###################################################
### code chunk number 17: tgp.Rnw:1115-1116
###################################################
sin.bgp <- bgp(X=X, Z=Z, XX=XX, verb=0)


###################################################
### code chunk number 18: sin-bgp
###################################################
plot(sin.bgp, main='GP,', layout='surf')
lines(X, Ztrue, col=4, lty=2, lwd=2)


###################################################
### code chunk number 19: tgp.Rnw:1124-1125
###################################################
graphics.off()


###################################################
### code chunk number 20: tgp.Rnw:1141-1142
###################################################
sin.btlm <- btlm(X=X, Z=Z, XX=XX)


###################################################
### code chunk number 21: sin-btlm
###################################################
plot(sin.btlm, main='treed LM,', layout='surf')
lines(X, Ztrue, col=4, lty=2, lwd=2)


###################################################
### code chunk number 22: tgp.Rnw:1155-1156
###################################################
graphics.off()


###################################################
### code chunk number 23: sin-btlmtrees
###################################################
tgp.trees(sin.btlm)


###################################################
### code chunk number 24: tgp.Rnw:1163-1164
###################################################
graphics.off()


###################################################
### code chunk number 25: tgp.Rnw:1185-1186
###################################################
sin.btgp <- btgp(X=X, Z=Z, XX=XX, verb=0)


###################################################
### code chunk number 26: sin-btgp
###################################################
plot(sin.btgp, main='treed GP,', layout='surf')
lines(X, Ztrue, col=4, lty=2, lwd=2)


###################################################
### code chunk number 27: tgp.Rnw:1194-1195
###################################################
graphics.off()


###################################################
### code chunk number 28: tgp.Rnw:1221-1222
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 29: tgp.Rnw:1240-1243
###################################################
exp2d.data <- exp2d.rand()
X <- exp2d.data$X; Z <- exp2d.data$Z
XX <- exp2d.data$XX


###################################################
### code chunk number 30: tgp.Rnw:1253-1254
###################################################
exp.bgp <- bgp(X=X, Z=Z, XX=XX, corr="exp", verb=0)   


###################################################
### code chunk number 31: exp-bgp
###################################################
plot(exp.bgp, main='GP,')


###################################################
### code chunk number 32: tgp.Rnw:1261-1262
###################################################
graphics.off()


###################################################
### code chunk number 33: tgp.Rnw:1285-1286
###################################################
exp.btgp <- btgp(X=X, Z=Z, XX=XX, corr="exp", verb=0)


###################################################
### code chunk number 34: exp-btgp
###################################################
plot(exp.btgp, main='treed GP,')


###################################################
### code chunk number 35: tgp.Rnw:1293-1294
###################################################
graphics.off()


###################################################
### code chunk number 36: exp-btgptrees
###################################################
tgp.trees(exp.btgp)


###################################################
### code chunk number 37: tgp.Rnw:1301-1302
###################################################
graphics.off()


###################################################
### code chunk number 38: tgp.Rnw:1326-1327
###################################################
exp.btgpllm <- btgpllm(X=X, Z=Z, XX=XX, corr="exp", R=2)  


###################################################
### code chunk number 39: exp-btgpllm
###################################################
plot(exp.btgpllm, main='treed GP LLM,')


###################################################
### code chunk number 40: tgp.Rnw:1334-1335
###################################################
graphics.off()


###################################################
### code chunk number 41: exp-1dbtgpllm1
###################################################
plot(exp.btgpllm, main='treed GP LLM,', proj=c(1))


###################################################
### code chunk number 42: tgp.Rnw:1357-1358
###################################################
graphics.off()


###################################################
### code chunk number 43: exp-1dbtgpllm2
###################################################
plot(exp.btgpllm, main='treed GP LLM,', proj=c(2))


###################################################
### code chunk number 44: tgp.Rnw:1364-1365
###################################################
graphics.off()


###################################################
### code chunk number 45: tgp.Rnw:1389-1390
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 46: tgp.Rnw:1413-1416
###################################################
library(MASS)
X <- data.frame(times=mcycle[,1])
Z <- data.frame(accel=mcycle[,2])


###################################################
### code chunk number 47: tgp.Rnw:1422-1423
###################################################
moto.bgp <- bgp(X=X, Z=Z, verb=0)


###################################################
### code chunk number 48: moto-bgp
###################################################
plot(moto.bgp, main='GP,', layout='surf')


###################################################
### code chunk number 49: tgp.Rnw:1431-1432
###################################################
graphics.off()


###################################################
### code chunk number 50: tgp.Rnw:1445-1446
###################################################
moto.btlm <- btlm(X=X, Z=Z, verb=0)


###################################################
### code chunk number 51: moto-btlm
###################################################
plot(moto.btlm, main='Bayesian CART,', layout='surf')


###################################################
### code chunk number 52: tgp.Rnw:1455-1456
###################################################
graphics.off()


###################################################
### code chunk number 53: tgp.Rnw:1473-1475
###################################################
moto.btgpllm <- btgpllm(X=X, Z=Z, bprior="b0", verb=0)
moto.btgpllm.p <- predict(moto.btgpllm) ## using MAP


###################################################
### code chunk number 54: moto-btgp
###################################################
par(mfrow=c(1,2))
plot(moto.btgpllm, main='treed GP LLM,', layout='surf')
plot(moto.btgpllm.p, center='km', layout='surf')


###################################################
### code chunk number 55: tgp.Rnw:1486-1487
###################################################
graphics.off()


###################################################
### code chunk number 56: moto-btgpq
###################################################
par(mfrow=c(1,2))
plot(moto.btgpllm, main='treed GP LLM,', layout='as')
plot(moto.btgpllm.p, as='ks2', layout='as')


###################################################
### code chunk number 57: tgp.Rnw:1497-1498
###################################################
graphics.off()


###################################################
### code chunk number 58: tgp.Rnw:1545-1546
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 59: tgp.Rnw:1577-1581
###################################################
f <- friedman.1.data(200)
ff <- friedman.1.data(1000)
X <- f[,1:10]; Z <- f$Y
XX <- ff[,1:10]


###################################################
### code chunk number 60: tgp.Rnw:1588-1591
###################################################
fr.btlm <- btlm(X=X, Z=Z, XX=XX, tree=c(0.95,2), pred.n=FALSE, verb=0)
fr.btlm.mse <- sqrt(mean((fr.btlm$ZZ.mean - ff$Ytrue)^2))
fr.btlm.mse


###################################################
### code chunk number 61: tgp.Rnw:1594-1597
###################################################
fr.bgpllm <- bgpllm(X=X, Z=Z, XX=XX, pred.n=FALSE, verb=0)
fr.bgpllm.mse <- sqrt(mean((fr.bgpllm$ZZ.mean - ff$Ytrue)^2))
fr.bgpllm.mse


###################################################
### code chunk number 62: tgp.Rnw:1606-1609
###################################################
XX1 <- matrix(rep(0,10), nrow=1)
fr.bgpllm.tr <- bgpllm(X=X, Z=Z, XX=XX1, pred.n=FALSE, trace=TRUE, 
                       m0r1=FALSE, verb=0)


###################################################
### code chunk number 63: tgp.Rnw:1619-1621
###################################################
trace <- fr.bgpllm.tr$trace$XX[[1]]
apply(trace[,27:36], 2, mean)


###################################################
### code chunk number 64: tgp.Rnw:1627-1628
###################################################
mean(fr.bgpllm.tr$trace$linarea$ba)


###################################################
### code chunk number 65: tgp.Rnw:1634-1635
###################################################
summary(trace[,9:10])


###################################################
### code chunk number 66: tgp.Rnw:1638-1639
###################################################
apply(trace[,11:15], 2, mean)


###################################################
### code chunk number 67: tgp.Rnw:1645-1646
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 68: tgp.Rnw:1652-1656
###################################################
exp2d.data <- exp2d.rand(lh=0, dopt=10)
X <- exp2d.data$X
Z <- exp2d.data$Z
Xcand <- lhs(1000, rbind(c(-2,6),c(-2,6)))


###################################################
### code chunk number 69: tgp.Rnw:1669-1670
###################################################
exp1 <- btgpllm(X=X, Z=Z, pred.n=FALSE, corr="exp", R=5, verb=0)


###################################################
### code chunk number 70: as-mapt
###################################################
tgp.trees(exp1)


###################################################
### code chunk number 71: tgp.Rnw:1677-1678
###################################################
graphics.off()


###################################################
### code chunk number 72: tgp.Rnw:1693-1695
###################################################
XX <- tgp.design(200, Xcand, exp1)
XX <- rbind(XX, c(-sqrt(1/2),0))


###################################################
### code chunk number 73: as-cands
###################################################
plot(exp1$X, pch=19, cex=0.5)
points(XX)
mapT(exp1, add=TRUE)


###################################################
### code chunk number 74: tgp.Rnw:1712-1713
###################################################
graphics.off()


###################################################
### code chunk number 75: tgp.Rnw:1727-1729
###################################################
exp.as <- btgpllm(X=X, Z=Z, XX=XX, corr="exp", improv=TRUE, 
                        Ds2x=TRUE, R=5, verb=0)


###################################################
### code chunk number 76: as-expas
###################################################
par(mfrow=c(1,3), bty="n")
plot(exp.as, main="tgpllm,", layout="as", as="alm")
plot(exp.as, main="tgpllm,", layout='as', as='alc')
plot(exp.as, main="tgpllm,", layout='as', as='improv')


###################################################
### code chunk number 77: tgp.Rnw:1747-1748
###################################################
graphics.off()


###################################################
### code chunk number 78: tgp.Rnw:1822-1823
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 79: tgp.Rnw:1859-1863
###################################################
exp2d.data <- exp2d.rand(n2=150, lh=0, dopt=10)
X <- exp2d.data$X
Z <- exp2d.data$Z
XX <- rbind(c(0,0),c(2,2),c(4,4))


###################################################
### code chunk number 80: tgp.Rnw:1868-1871
###################################################
out <- btgpllm(X=X, Z=Z, XX=XX, corr="exp", bprior="b0", 
               pred.n=FALSE, Ds2x=TRUE, R=10, 
               trace=TRUE, verb=0)


###################################################
### code chunk number 81: tgp.Rnw:1875-1876
###################################################
out$trace


###################################################
### code chunk number 82: traces-XXd
###################################################
trXX <- out$trace$XX; ltrXX <- length(trXX)
y <- trXX[[1]]$d
for(i in 2:ltrXX) y <- c(y, trXX[[i]]$d)
plot(log(trXX[[1]]$d), type="l", ylim=range(log(y)), ylab="log(d)",
     main="range (d) parameter traces")
names <- "XX[1,]"
for(i in 2:ltrXX) {
  lines(log(trXX[[i]]$d), col=i, lty=i)
  names <- c(names, paste("XX[", i, ",]", sep=""))
}
legend("bottomleft", names, col=1:ltrXX, lty=1:ltrXX)


###################################################
### code chunk number 83: tgp.Rnw:1908-1909
###################################################
graphics.off()


###################################################
### code chunk number 84: tgp.Rnw:1926-1928
###################################################
linarea <- mean(out$trace$linarea$la)
linarea


###################################################
### code chunk number 85: traces-la
###################################################
hist(out$trace$linarea$la)


###################################################
### code chunk number 86: tgp.Rnw:1935-1936
###################################################
graphics.off()


###################################################
### code chunk number 87: tgp.Rnw:1951-1957
###################################################
m <- matrix(0, nrow=length(trXX), ncol=3)#ncol=5)
for(i in 1:length(trXX))
  m[i,] <- as.double(c(out$XX[i,], mean(trXX[[i]]$b)))
m <- data.frame(cbind(m, 1-m[,3]))
names(m)=c("XX1","XX2","b","pllm")
m


###################################################
### code chunk number 88: traces-alc
###################################################
trALC <- out$trace$preds$Ds2x
y <- trALC[,1]
for(i in 2:ncol(trALC)) y <- c(y, trALC[,i])
plot(log(trALC[,1]), type="l", ylim=range(log(y)), ylab="Ds2x",
     main="ALC: samples from Ds2x")
names <- "XX[1,]"
for(i in 2:ncol(trALC)) {
  lines(log(trALC[,i]), col=i, lty=i)
  names <- c(names, paste("XX[", i, ",]", sep=""))
}
legend("bottomright", names, col=1:ltrXX, lty=1:ltrXX)


###################################################
### code chunk number 89: tgp.Rnw:1978-1979
###################################################
graphics.off()


###################################################
### code chunk number 90: tgp.Rnw:2065-2066
###################################################
seed <- 0; set.seed(seed)


###################################################
### code chunk number 91: tgp.Rnw:2076-2081
###################################################
library(MASS)
out <- btgpllm(X=mcycle[,1], Z=mcycle[,2], bprior="b0", 
         pred.n=FALSE, verb=0) 
save(out, file="out.Rsave")
out <- NULL


###################################################
### code chunk number 92: tgp.Rnw:2090-2093
###################################################
load("out.Rsave")
XX <- seq(2.4, 56.7, length=200)
out.kp <- predict(out, XX=XX, pred.n=FALSE)


###################################################
### code chunk number 93: tgp.Rnw:2098-2099
###################################################
out.p <- predict(out, XX=XX, pred.n=FALSE, BTE=c(0,1000,1))


###################################################
### code chunk number 94: tgp.Rnw:2108-2109
###################################################
out2 <- predict(out, XX, pred.n=FALSE, BTE=c(0,2000,2), MAP=FALSE)


###################################################
### code chunk number 95: pred-kp
###################################################
plot(out.kp, center="km", as="ks2")


###################################################
### code chunk number 96: tgp.Rnw:2120-2121
###################################################
graphics.off()


###################################################
### code chunk number 97: pred-p
###################################################
plot(out.p)


###################################################
### code chunk number 98: tgp.Rnw:2128-2129
###################################################
graphics.off()


###################################################
### code chunk number 99: pred-2
###################################################
plot(out2)


###################################################
### code chunk number 100: tgp.Rnw:2136-2137
###################################################
graphics.off()


###################################################
### code chunk number 101: tgp.Rnw:2160-2161
###################################################
unlink("out.Rsave")