Flexible Correlation in Bayesian Gaussian Process Emulation

To cite: Chen, H., Loeppky, J.L., and Welch, W.J. (2017), Flexible Correlation Structure for Accurate Prediction and Uncertainty Quantification in Bayesian Gaussian Process Emulation of a Computer Model, SIAM/ASA Journal on Uncertainty Quantification, in press.

Coded by Hao Chen under the supervision of William J. Welch.
Last change was made on May 4, 2017.

Illustrated using the N-K example (Bastos and O’Hagan, 2009). The following demonstration code and all files called are in GPflex.zip.

1. Read the data

dp = 5
x1 = read.table(paste(getwd(), "/data/x100.txt", sep = ""), header = F)
colnames(x1) = paste('x', c(1:dp), sep = '')
y1 = read.table(paste(getwd(), "/data/logy100.txt", sep = ""), header = F)[,  2]

x2 = read.table(paste(getwd(), "/data/x150.txt", sep = ""),header = F)
colnames(x2) = paste('x', c(1:dp), sep = '')
y2 = read.table(paste(getwd(), "/data/logy150.txt", sep = ""), header = F)[, 2]

2. Get the training set and hold-out set.

set.seed(100 + 69)
tst = sample(1:150, 100, replace = F)

Training set, randomly sample 100 from 150, n = 100

xf = x2[tst, ]
rownames(xf) = c(1:100)
yf = data.frame(y = y2[tst]) 
rownames(yf) = c(1:100)

Hold-out set, m = 100 + 50 = 150

xp = rbind(x1, x2[-tst, ])
rownames(xp) = c(1:150)
yp = data.frame(y = c(y1, y2[-tst]))

3. Model Fitting

3.1 Hybrid with PowExp

source(paste(getwd(),"/Hybrid_PowExp/hybrid_powexp.R", sep = ""))
########smoothness parameters
alpha = c(1.0086779, 2, 1.7810691, 1.9886917, 1.8338945)
########Fit the model
fit_h<-runBMCMC(nmcmc = 10000, burn = 4000, thin = 10, x = xf, y = yf, xtest1 = xp, 
                lambda.ini = rep(0.5, dp), lambda.w.ini = rep(0.1, dp), alpha)

## [1] 0.1
## [1] 0.2
## [1] 0.3
## [1] 0.4
## [1] 0.5
## [1] 0.6
## [1] 0.7
## [1] 0.8
## [1] 0.9
## [1] 1

RMSE

print(sqrt(sum((fit_h$pred.y - yp$y)^2) / length(yp$y)))

## [1] 0.05000984

Normalized RMSE

print(sqrt(sum((fit_h$pred.y-yp$y)^2) / length(yp$y)) / sqrt(sum((mean(yf$y)-yp$y)^2) / length(yp$y)))

## [1] 0.09915959

3.2 Full with PowExp

source(paste(getwd(),"/Full_PowExp/full_powexp.R", sep = ""))
########Fit the model
fit_f<-runBMCMC(nmcmc = 10000, burn = 4000, thin = 10, x = xf, y = yf, xtest1 = xp, 
                lambda.ini = rep(0.5, dp), lambda.w.ini = rep(0.1, dp), gamma.ini = rep(5, dp), 
                gamma.w.ini = rep(20, dp))

## [1] 0.1
## [1] 0.2
## [1] 0.3
## [1] 0.4
## [1] 0.5
## [1] 0.6
## [1] 0.7
## [1] 0.8
## [1] 0.9
## [1] 1

RMSE

print(sqrt(sum((fit_f$pred.y-yp$y)^2)/length(yp$y)))

## [1] 0.04970327

Normalized RMSE

print(sqrt(sum((fit_f$pred.y - yp$y)^2) / length(yp$y)) / sqrt(sum((mean(yf$y) - yp$y)^2) / length(yp$y)))

## [1] 0.09855173

3.3 Hybrid with Matern

source(paste(getwd(), "/Hybrid_Matern/hybrid_matern.R", sep = ""))
########smoothness parameters
cp = c(0, 2, 1, 1, 1)
########Fit the model
fit_m<-runBMCMC(nmcmc = 10000, burn = 4000, thin = 10, x = xf, y = yf, xtest1 = xp, 
               lambda.ini = rep(0.5, dp), lambda.w.ini = rep(0.1, dp), cp)

## [1] 0.1
## [1] 0.2
## [1] 0.3
## [1] 0.4
## [1] 0.5
## [1] 0.6
## [1] 0.7
## [1] 0.8
## [1] 0.9
## [1] 1

RMSE

print(sqrt(sum((fit_m$pred.y-yp$y)^2) / length(yp$y)))

## [1] 0.04860369

Normalized RMSE

print(sqrt(sum((fit_m$pred.y - yp$y)^2) / length(yp$y)) / sqrt(sum((mean(yf$y) - yp$y)^2) / length(yp$y)))

## [1] 0.09637148