### 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")