Objective :
Related material:
Readings :
- Rasmussen, Carl Edward. "Gaussian processes in machine learning." In Advanced lectures on machine learning, pp. 63-71. Springer Berlin Heidelberg, 2004.
- see http://www.GaussianProcess.org/gpml
- Chapters: 2, 4, 5.1, & 5.4.2
- Slides provided
Software :
- R-cran (https://cran.r-project.org/)
- R packages
- DiceKrigin (https://cran.r-project.org/web/packages/DiceKriging/index.html)
- lhs (https://cran.r-project.org/web/packages/lhs/index.html)
- Roustant, Olivier, David Ginsbourger, and Yves Deville. "DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization." (2012).
- Plus dependences ...The Piston Simulation function models the circular motion of a piston within a cylinder. It involves a chain of nonlinear functions.
The response C is cycle time (the time it takes to complete one cycle), in seconds.
$$ C(x) = 2\pi \sqrt{\frac{M}{k+S^2}\frac{P_0 V_0}{T_0}\frac{T_0}{V^2}} $$$$ V = \frac{S}{2k} ( \sqrt{A^2 +4k\frac{P_0 V_0}{T_0}T_a} -A) $$$$ A = P_0 S 19.62M-\frac{kV_0}{S} $$The input variables and their usual input ranges are:
-------------------------------------------------- M ∈ [30, 60] | piston weight (kg) S ∈ [0.005, 0.020] | piston surface area (m2) V0 ∈ [0.002, 0.010] | initial gas volume (m3) k ∈ [1000, 5000] | spring coefficient (N/m) P0 ∈ [90000, 110000] | atmospheric pressure (N/m2) Ta ∈ [290, 296] | ambient temperature (K) T0 ∈ [340, 360] | filling gas temperature (K) --------------------------------------------------
Here, we pretend that parameters M, k, P0, Ta, and T0 are fixed and equal to the mid-value of their ranges (M=45, k=3000, P0=1e+5000, Ta=293, and T0=350), and that the only free input parameters are S, and V0. Furthermore, we consider the output in the log scale.
=> As a result we wish to emulate the function f(S,V0) := log(C(S,V0)).
Reference: http://www.sfu.ca/~ssurjano/piston.html
Kenett, R., & Zacks, S. (1998). Modern industrial statistics: design and control of quality and reliability. Pacific Grove, CA: Duxbury press.
In [2]:
# DOWNLOAD THE R PACKAGES REQUIRED
install.packages('DiceKriging', repos = "http://cran.us.r-project.org")
install.packages('lhs', repos = "http://cran.us.r-project.org")
# install.packages('tcltk', repos = "http://cran.us.r-project.org")
# install.packages('aplpack', repos = "http://cran.us.r-project.org")
In [89]:
# LOAD THE R PACKAGES REQUIRED
library('lhs')
library('DiceKriging')
# library('tcltk')
# library('aplpack')
In [90]:
# Consider the simulator
piston <- function(xx) {
##########################################################################
#
# PISTON FUNCTION
#
# Authors: Sonja Surjanovic, Simon Fraser University
# Derek Bingham, Simon Fraser University
# Questions/Comments: Please email Derek Bingham at dbingham@stat.sfu.ca.
#
# Copyright 2013. Derek Bingham, Simon Fraser University.
#
# THERE IS NO WARRANTY, EXPRESS OR IMPLIED. WE DO NOT ASSUME ANY LIABILITY
# FOR THE USE OF THIS SOFTWARE. If software is modified to produce
# derivative works, such modified software should be clearly marked.
# Additionally, this program is free software; you can redistribute it
# and/or modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; version 2.0 of the License.
# Accordingly, this program is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the implied warranty
# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# For function details and reference information, see:
# http://www.sfu.ca/~ssurjano/
#
##########################################################################
#
# OUTPUT AND INPUT:
#
# C = cycle time
# xx = c(M, S, V0, k, P0, Ta, T0)
#
##########################################################################
M <- xx[1]
S <- xx[2]
V0 <- xx[3]
k <- xx[4]
P0 <- xx[5]
Ta <- xx[6]
T0 <- xx[7]
Aterm1 <- P0 * S
Aterm2 <- 19.62 * M
Aterm3 <- -k*V0 / S
A <- Aterm1 + Aterm2 + Aterm3
Vfact1 <- S / (2*k)
Vfact2 <- sqrt(A^2 + 4*k*(P0*V0/T0)*Ta)
V <- Vfact1 * (Vfact2 - A)
fact1 <- M
fact2 <- k + (S^2)*(P0*V0/T0)*(Ta/(V^2))
C <- 2 * pi * sqrt(fact1/fact2)
return(C)
}
input_min <- c(30, 0.005, 0.002, 1000, 90000, 290, 340 )
input_max <- c(60, 0.020, 0.020, 5000, 110000, 296, 360 )
myfun <- piston
# Here:
# - we consider only 2 inputs for the presentation purpose and fix the rest
# - wetransform the output in the log scale
# Namely :
myfun <- function(xx) {
zz <- c(NaN,NaN,NaN,NaN,NaN,NaN,NaN)
zz[1] <- M <- 0.5*(30 + 60)
zz[2] <- S <- 0.5*(0.005 + 0.020)
zz[3] <- V0 <- 0.5*(0.002 + 0.010)
zz[4] <- k <- 0.5*(1000 + 5000)
zz[5] <- P0 <- 0.5*(90000 + 110000)
zz[6] <- Ta <- 0.5*(290 + 296)
zz[7] <- T0 <- 0.5*(340 + 360)
zz[2] <- Ta <- xx[1]
zz[3] <- T0 <- xx[2]
C <- piston(zz)
C <- log(C)
return(C)
}
input_min <- c(0.005 , 0.002)
input_max <- c(0.020 , 0.010)
n.grid <- 100 ;
x1.grid <-seq(input_min[1],input_max[1],length.out=n.grid) ;
x2.grid <-seq(input_min[2],input_max[2],length.out=n.grid) ;
X.grid <- expand.grid( x1=x1.grid, x2=x2.grid )
y.grid <- apply(X.grid,1,myfun)
# Contour plot
contour(x1.grid, x2.grid, matrix(y.grid, n.grid, n.grid), 10,
main = "Real function",
xlab = "x1",
ylab = "x2",
xlim = c(input_min[1],input_max[1]),
ylim = c(input_min[2],input_max[2]))
# Surface plot
persp(x1.grid, x2.grid, matrix(y.grid, n.grid, n.grid ),
xlab = "x1", ylab = "x2", zlab = "f(x)",
theta = 30, phi = 30,
xlim = c(input_min[1],input_max[1]),
ylim = c(input_min[2],input_max[2]),
col='blue',
main="Real function")
In [104]:
# Generate a training data-set
n_data <- 7 ;
n_dim <- 2
X_data <- t(input_min + (input_max-input_min)*t(optimumLHS(n=n_data, k=n_dim))) ;
Y_data <- apply(X_data, 1, myfun) ;
plot(X_data, xlab = "x1", ylab = "x2", main="Design points")
In [88]:
# We use the function:
# km {DiceKriging}
?km
Out[88]:
In [136]:
# Diagnostics (Feel free to use others that you know)
R2 <-function (Y, Ypred) {
Ypred <- as.numeric(Ypred)
Y <- as.numeric(Y)
return(1 - mean((Y - Ypred)^2)/mean((Y - mean(Y))^2))
}
RMSE <- function (Y, Ypred) {
Ypred <- as.numeric(Ypred)
Y <- as.numeric(Y)
return(sqrt(mean((Y - Ypred)^2)))
}
MAE <- function (Y, Ypred) {
Ypred <- as.numeric(Ypred)
Y <- as.numeric(Y)
return(mean(abs(Y - Ypred)))
}
In [105]:
# Train different GP regressions differing on the covariance function
myfun_km_matern5_2 <- km( formula = ~ 1 , # linear trend formula; E.G.: "formula = ~ 1 +x1 +exp(x2) +I(x1*x2) +I(x1^2)"
design = data.frame(x1=X_data[,1],
x2=X_data[,2]), # a data frame representing the design of experiments.
response=Y_data, # the values of the 1-dimensional outpu
covtype="matern5_2", # covariance structure
coef.trend = NULL, # values for the trend parameters
coef.cov = NULL, # values for the covariance parameters
coef.var = NULL, # values for the variance parameters
nugget= 1e-7, # the homogeneous nugget effect
noise.var = NULL, # containing the noise variance at each observation
optim.method = "BFGS", # optimization method is chosen for the likelihood maximization.
lower = NULL, # lower bounds of the correlation parameters for optimization
upper = NULL, # upper bounds of the correlation parameters for optimizati
control =list(trace=FALSE), # =list(trace=FALSE) to supress optimization trace information
kernel=NULL) # a user's new covariance structure
myfun_km_matern3_2 <- km( formula = ~ 1 , # linear trend formula; E.G.: "formula = ~ 1 +x1 +exp(x2) +I(x1*x2) +I(x1^2)"
design = data.frame(x1=X_data[,1],
x2=X_data[,2]), # a data frame representing the design of experiments.
response=Y_data, # the values of the 1-dimensional outpu
covtype="matern3_2", # covariance structure
coef.trend = NULL, # values for the trend parameters
coef.cov = NULL, # values for the covariance parameters
coef.var = NULL, # values for the variance parameters
nugget= 1e-7, # the homogeneous nugget effect
noise.var = NULL, # containing the noise variance at each observation
optim.method = "BFGS", # optimization method is chosen for the likelihood maximization.
lower = NULL, # lower bounds of the correlation parameters for optimization
upper = NULL, # upper bounds of the correlation parameters for optimizati
control =list(trace=FALSE), # =list(trace=FALSE) to supress optimization trace information
kernel=NULL)
myfun_km_gauss <- km( formula = ~ 1 , # linear trend formula; E.G.: "formula = ~ 1 +x1 +exp(x2) +I(x1*x2) +I(x1^2)"
design = data.frame(x1=X_data[,1],
x2=X_data[,2]), # a data frame representing the design of experiments.
response=Y_data, # the values of the 1-dimensional outpu
covtype="gauss", # covariance structure
coef.trend = NULL, # values for the trend parameters
coef.cov = NULL, # values for the covariance parameters
coef.var = NULL, # values for the variance parameters
nugget= 1e-7, # the homogeneous nugget effect
noise.var = NULL, # containing the noise variance at each observation
optim.method = "BFGS", # optimization method is chosen for the likelihood maximization.
lower = NULL, # lower bounds of the correlation parameters for optimization
upper = NULL, # upper bounds of the correlation parameters for optimizati
control =list(trace=FALSE), # =list(trace=FALSE) to supress optimization trace information
kernel=NULL) # a user's new covariance structure
myfun_km_exp <- km( formula = ~ 1 , # linear trend formula; E.G.: "formula = ~ 1 +x1 +exp(x2) +I(x1*x2) +I(x1^2)"
design = data.frame(x1=X_data[,1],
x2=X_data[,2]), # a data frame representing the design of experiments.
response=Y_data, # the values of the 1-dimensional outpu
covtype="exp", # covariance structure
coef.trend = NULL, # values for the trend parameters
coef.cov = NULL, # values for the covariance parameters
coef.var = NULL, # values for the variance parameters
nugget= 1e-7, # the homogeneous nugget effect
noise.var = NULL, # containing the noise variance at each observation
optim.method = "BFGS", # optimization method is chosen for the likelihood maximization.
lower = NULL, # lower bounds of the correlation parameters for optimization
upper = NULL, # upper bounds of the correlation parameters for optimizati
control =list(trace=FALSE), # =list(trace=FALSE) to supress optimization trace information
kernel=NULL) # a user's new covariance structure
# CROSS VALIDATION CRITERIA
Y_pred <- leaveOneOut.km(myfun_km_matern5_2, "UK")$mean
CV_matern5_2 <-data.frame( R2 = R2(Y_data, Y_pred),
RMSE = RMSE(Y_data, Y_pred),
MAE = MAE(Y_data, Y_pred), row.names="matern5_2")
Y_pred <- leaveOneOut.km(myfun_km_matern3_2, "UK")$mean
CV_matern3_2 <-data.frame( R2 = R2(Y_data, Y_pred),
RMSE = RMSE(Y_data, Y_pred),
MAE = MAE(Y_data, Y_pred), row.names="matern3_2")
Y_pred <- leaveOneOut.km(myfun_km_gauss, "UK")$mean
CV_gauss <-data.frame( R2 = R2(Y_data, Y_pred),
RMSE = RMSE(Y_data, Y_pred),
MAE = MAE(Y_data, Y_pred), row.names="gauss~~~~")
Y_pred <- leaveOneOut.km(myfun_km_exp, "UK")$mean
CV_exp <-data.frame( R2 = R2(Y_data, Y_pred),
RMSE = RMSE(Y_data, Y_pred),
MAE = MAE(Y_data, Y_pred), row.names="exp~~~~~~")
CV_matern5_2
CV_matern3_2
CV_gauss
CV_exp
Out[105]:
Out[105]:
Out[105]:
Out[105]:
In [106]:
# Train the GP regression
myfun_km <- myfun_km_gauss
In [137]:
# CROSS VALIDATION CRITERIA
Y_pred <- leaveOneOut.km(myfun_km, "UK")$mean
data.frame(
R2 = R2(Y_data, Y_pred),
RMSE = RMSE(Y_data, Y_pred),
MAE = MAE(Y_data, Y_pred))
# PLOT DIAGNOSTICS
plot(myfun_km)
Out[137]:
In [107]:
# PRINT THE PARAMETERS
show(myfun_km)
In [140]:
# Plot the
n.grid <- 100 ;
x1.grid <-seq(input_min[1],input_max[1],length.out=n.grid) ;
x2.grid <-seq(input_min[2],input_max[2],length.out=n.grid) ;
X.grid <- expand.grid( x1=x1.grid, x2=x2.grid )
y.grid <- apply(X.grid, 1, myfun)
y.pred.grid <- predict(myfun_km, X.grid, "UK")
par(mfrow = c(3,3))
contour(x1.grid, x2.grid, matrix(y.grid, n.grid, n.grid), 11, main = "Real function")
points(X_data[ , 1], X_data[ , 2], pch = 19, cex = 1.5, col = "red")
contour(x1.grid, x2.grid, matrix(y.pred.grid$mean, n.grid, n.grid), 11, main = "Posterior GP regression mean")
points(X_data[ , 1], X_data[ , 2], pch = 19, cex = 1.5, col = "red")
contour(x1.grid, x2.grid, matrix(y.pred.grid$sd^2, n.grid, n.grid), 15, main = "Posterior GP regression variance")
points(X_data[ , 1], X_data[ , 2], pch = 19, cex = 1.5, col = "red")