# input space based on latin-hypercube sampling (not required)
# two dimensional example with N=216 sized sample
if(require(lhs)) { X <- randomLHS(216, 2)
} else { X <- matrix(runif(216*2), ncol=2) }
# pseudo responses, not important for visualizing design
Y <- runif(216)
## BLHS sample with m=6 divisions in each coordinate
sub <- blhs(y=Y, X=X, m=6)
Xsub <- sub$xs # the bootstrapped subsample
# visualization
plot(X, xaxt="n", yaxt="n", xlim=c(0,1), ylim=c(0,1), xlab="factor 1",
ylab="factor 2", col="cyan", main="BLHS")
b <- seq(0, 1, by=1/6)
abline(h=b, v=b, col="black", lty=2)
axis(1, at=seq (0, 1, by=1/6), cex.axis=0.8,
labels=expression(0, 1/6, 2/6, 3/6, 4/6, 5/6, 1))
axis(2, at=seq (0, 1, by=1/6), cex.axis=0.8,
labels=expression(0, 1/6, 2/6, 3/6, 4/6, 5/6, 1), las=1)
points(Xsub, col="red", pch=19, cex=1.25)
## Comparing global lengthscale MLE based on BLHS and random subsampling
if (FALSE) {
# famous borehole function
borehole <- function(x){
rw <- x[1] * (0.15 - 0.05) + 0.05
r <- x[2] * (50000 - 100) + 100
Tu <- x[3] * (115600 - 63070) + 63070
Tl <- x[5] * (116 - 63.1) + 63.1
Hu <- x[4] * (1110 - 990) + 990
Hl <- x[6] * (820 - 700) + 700
L <- x[7] * (1680 - 1120) + 1120
Kw <- x[8] * (12045 - 9855) + 9855
m1 <- 2 * pi * Tu * (Hu - Hl)
m2 <- log(r / rw)
m3 <- 1 + 2*L*Tu/(m2*rw^2*Kw) + Tu/Tl
return(m1/m2/m3)
}
N <- 100000 # number of observations
if(require(lhs)) { xt <- randomLHS(N, 8) # input space
} else { xt <- matrix(runif(N*8), ncol=8) }
yt <- apply(xt, 1, borehole) # response
colnames(xt) <- c("rw", "r", "Tu", "Tl", "Hu", "Hl", "L", "Kw")
## prior on the GP lengthscale parameter
da <- darg(list(mle=TRUE, max=100), xt)
## make space for two sets of boxplots
par(mfrow=c(1,2))
# BLHS calculating with visualization of the K MLE lengthscale estimates
K <- 10 # number of Bootstrap samples; Sun, et al (2017) uses K <- 31
sub_blhs <- blhs.loop(y=yt, X=xt, K=K, m=2, da=da, maxit=200, plot.it=TRUE)
# a random subsampling analog for comparison
sn <- sub_blhs$ly # extract a size that is consistent with the BLHS
that.rand <- matrix(NA, ncol=8, nrow=K)
for(i in 1:K){
sub <- sample(1:nrow(xt), sn)
gpsepi <- newGPsep(xt[sub,], yt[sub], d=da$start, g=1e-3, dK=TRUE)
mle <- mleGPsep(gpsepi, tmin=da$min, tmax=10*da$max, ab=da$ab, maxit=200)
deleteGPsep(gpsepi)
that.rand[i,] <- mle$d
}
## put random boxplots next to BLHS ones
boxplot(that.rand, xlab="input", ylab="theta-hat", col=2,
main="random subsampling")
}
Run the code above in your browser using DataLab