data(ato)
if (FALSE) {
##
## the code below was used to create the random partition
##
## recover the data in its original form
X <- X*19+1
Z <- Z*sqrt(Zv) + Zm
## code the inputs and outputs; i.e., undo the transformation
## above
X <- (X-1)/19
Zm <- mean(Z)
Zv <- var(as.vector(Z))
Z <- (Z - Zm)/sqrt(Zv)
## random training and testing partition
train <- sample(1:nrow(X), 1000)
Xtrain <- X[train,]
Xtest <- X[-train,]
Ztest <- as.list(as.data.frame(t(Z[-train,])))
Ztrain <- Ztrain.out <- list()
mult <- rep(NA, nrow(Xtrain))
kill <- rep(FALSE, nrow(Xtrain))
for(i in 1:length(train)) {
reps <- sample(1:ncol(Z), 1)
w <- sample(1:ncol(Z), reps)
Ztrain[[i]] <- Z[train[i],w]
if(reps < 10) Ztrain.out[[i]] <- Z[train[i],-w]
else kill[i] <- TRUE
mult[i] <- reps
}
## calculate training locations and outputs for replicates not
## included in Ztrain
Xtrain.out <- Xtrain[!kill,]
Ztrain.out <- Ztrain[which(!kill)]
## fit hetGP model
out <- mleHetGP(X=list(X0=Xtrain, Z0=sapply(Ztrain, mean), mult=mult),
Z=unlist(Ztrain), lower=rep(0.01, ncol(X)), upper=rep(30, ncol(X)),
covtype="Matern5_2", noiseControl=nc, known=list(beta0=0),
maxit=100000, settings=list(return.matrices=FALSE))
##
## the adaptive lookahead design is read in and fit as
## follows
##
Xa <- (ato.a[,1:8]-1)/19
Za <- ato.a[,9]
Za <- (Za - Zm)/sqrt(Zv)
## uses nc defined above
out.a <- mleHetGP(Xa, Za, lower=rep(0.01, ncol(X)),
upper=rep(30, ncol(X)), covtype="Matern5_2", known=list(beta0=0),
noiseControl=nc, maxit=100000, settings=list(return.matrices=FALSE))
}
##
## the following code duplicates a predictive comparison in
## the package vignette
##
## first using the model fit to the train partition (out)
out <- rebuild(out)
## predicting out-of-sample at the test sights
phet <- predict(out, Xtest)
phets2 <- phet$sd2 + phet$nugs
mhet <- as.numeric(t(matrix(rep(phet$mean, 10), ncol=10)))
s2het <- as.numeric(t(matrix(rep(phets2, 10), ncol=10)))
sehet <- (unlist(t(Ztest)) - mhet)^2
sc <- - sehet/s2het - log(s2het)
mean(sc)
## predicting at the held-out training replicates
phet.out <- predict(out, Xtrain.out)
phets2.out <- phet.out$sd2 + phet.out$nugs
s2het.out <- mhet.out <- Ztrain.out
for(i in 1:length(mhet.out)) {
mhet.out[[i]] <- rep(phet.out$mean[i], length(mhet.out[[i]]))
s2het.out[[i]] <- rep(phets2.out[i], length(s2het.out[[i]]))
}
mhet.out <- unlist(t(mhet.out))
s2het.out <- unlist(t(s2het.out))
sehet.out <- (unlist(t(Ztrain.out)) - mhet.out)^2
sc.out <- - sehet.out/s2het.out - log(s2het.out)
mean(sc.out)
if (FALSE) {
## then using the model trained from the "adaptive"
## sequential design, with comparison from the "batch"
## one above, using the scores function
out.a <- rebuild(out.a)
sc.a <- scores(out.a, Xtest = Xtest, Ztest = Ztest)
c(batch=mean(sc), adaptive=sc.a)
## an example of one iteration of sequential design
Wijs <- Wij(out.a$X0, theta=out.a$theta, type=out.a$covtype)
h <- horizon(out.a, Wijs=Wijs)
control = list(tol_dist=1e-4, tol_diff=1e-4, multi.start=30, maxit=100)
opt <- IMSPE_optim(out.a, h, Wijs=Wijs, control=control)
opt$par
}
Run the code above in your browser using DataLab