| Title: | Quantile Kriging for Stochastic Simulations with Replication |
|---|---|
| Description: | A re-implementation of quantile kriging. Quantile kriging was described by Plumlee and Tuo (2014) <doi:10.1080/00401706.2013.860919>. With computational savings when dealing with replication from the recent paper by Binois, Gramacy, and Ludovski (2018) <doi:10.1080/10618600.2018.1458625> it is now possible to apply quantile kriging to a wider class of problems. In addition to fitting the model, other useful tools are provided such as the ability to automatically perform leave-one-out cross validation. |
| Authors: | Kevin R. Quinlan [aut, cre], Jim R. Leek [aut], Lawrence Livermore National Security [cph] |
| Maintainer: | Kevin R. Quinlan <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2025-02-21 05:00:33 UTC |
| Source: | https://github.com/cran/quantkriging |
Generates Leave-One-Out predictions for each location and quantile.
LOOquants(QKResults)LOOquants(QKResults)
QKResults |
Output from the quantKrig function. |
Returns the estimated quantiles and a plot of the leave-one-out predictions against the observed values.
Leave-one-out predictions at the input locations
X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) LOOquants(Qout)X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) LOOquants(Qout)
Create Quantile Kriging Results class from list
new_QKResults(qkList)new_QKResults(qkList)
qkList |
list(quants, yquants, g, l, ll <- -optparb$value, beta0, nu, xstar, ystar, Ki, quantv, mult, ylisto, type) |
New class QKResults
Generates new quantiles from a quantile object
newQuants(QKResults, quantv)newQuants(QKResults, quantv)
QKResults |
Output from the quantKrig function. |
quantv |
Vector of quantile values alpha between 0 and 1, |
The same quantile object with new estimated quantiles.
X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) Qout2 <- newQuants(Qout, c(0.025, 0.5, 0.975)) QuantPlot(Qout2)X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) Qout2 <- newQuants(Qout, c(0.025, 0.5, 0.975)) QuantPlot(Qout2)
Quantile Predictions using Quantile Kriging model (class QKResults)
## S3 method for class 'QKResults' predict(object, xnew, quantnew = NULL, ...)## S3 method for class 'QKResults' predict(object, xnew, quantnew = NULL, ...)
object |
Output from the quantKrig function. |
xnew |
Locations for prediction |
quantnew |
Quantiles for prediction, default is to keep the same as the quantile object |
... |
Ignore. No other arguments for this method |
Quantile predictions at the specified input locations
Kevin Quinlan [email protected]
X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) predict(Qout, xnew = c(0.4, 0.5, 0.6)) quantpreds <- predict(Qout, xnew = seq(0,1,length.out = 100), quantnew = seq(0.01,0.99,by = 0.01)) matplot(seq(0,1,length.out = 100), quantpreds, type = 'l')X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) predict(Qout, xnew = c(0.4, 0.5, 0.6)) quantpreds <- predict(Qout, xnew = seq(0,1,length.out = 100), quantnew = seq(0.01,0.99,by = 0.01)) matplot(seq(0,1,length.out = 100), quantpreds, type = 'l')
Implements Quantile Kriging from Plumlee and Tuo (2014).
quantKrig(x, y, quantv, lower, upper, method = "loo", type = "Gaussian", rs = TRUE, nm = TRUE, known = NULL, optstart = NULL, control = list())quantKrig(x, y, quantv, lower, upper, method = "loo", type = "Gaussian", rs = TRUE, nm = TRUE, known = NULL, optstart = NULL, control = list())
x |
Inputs |
y |
Univariate Response |
quantv |
Vector of Quantile values to estimate (ex: c(0.025, 0.975)) |
lower |
Lower bound of hyperparameters, if isotropic set lengthscale then nugget, if anisotropic set k lengthscales and then nugget |
upper |
Upper bound of hyperparameters, if isotropic set lengthscale then nugget, if anisotropic set k lengthscales and then nugget |
method |
Either maximum likelihood ('mle') or leave-one-out cross validation ('loo') optimization of hyperparameters |
type |
Covariance type, either 'Gaussian', 'Matern3_2', or 'Matern5_2' |
rs |
If TRUE, rescales inputs to [0,1] |
nm |
If TRUE, normalizes output to mean 0, variance 1 |
known |
Fixes all hyperparamters to a known value |
optstart |
Sets the starting value for the optimization |
control |
Control from optim function |
Fits quantile kriging using a double exponential or Matern covariance function. This emulator is for a stochastic simulation and models the distribution of the results (through the quantiles), not just the mean. The hyperparameters can be trained using maximum likelihood estimation or leave-one-out cross validation as recommended in Plumlee and Tuo (2014). The GP is trained using the Woodbury formula to improve computation speed with replication as shown in Binois et al. (2018). To get meaningful results, there should be sufficient replication at each input. The quantiles at a location are found using:
) where is the kernel of the design matrix (with nugget effect), the ordered sample closest to that quantile at each input, and the mean at each input.
The estimated quantile values in matrix form
The actual quantile values from the data in matrix form
The scaling parameter for the kernel
The lengthscale parameter(s)
The log likelihood
Estimated linear trend
Estimator of the variance
Matrix of unique input values
Average value at each unique input value
Inverted covariance matrix
Vector of alpha values between 0 and 1 for estimated quantiles, it is recommended that only a small number of quantiles are used for fitting and more quantiles can be found later using newQuants
Number of replicates at each input
Matthew Plumlee & Rui Tuo (2014) Building Accurate Emulators for Stochastic Simulations via Quantile Kriging, Technometrics, 56:4, 466-473, DOI: 10.1080/00401706.2013.860919
Mickael Binois, Robert B. Gramacy & Mike Ludkovski (2018) Practical Heteroscedastic Gaussian Process Modeling for Large Simulation Experiments, Journal of Computational and Graphical Statistics, 27:4, 808-821, DOI: 10.1080/10618600.2018.1458625
# Simple example X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) Ystar <- rnorm(length(Ystar),Ystar,1) lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, quantv = seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar) #fit for non-normal errors Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e Qout <- quantKrig(Xstar,Ystar, quantv = seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar)# Simple example X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) Ystar <- rnorm(length(Ystar),Ystar,1) lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, quantv = seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar) #fit for non-normal errors Ystar <- rep(Y,each = 100) e <- rchisq(length(Ystar),5)/5 - 1 Ystar <- Ystar + e Qout <- quantKrig(Xstar,Ystar, quantv = seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar)
Quantile Kriging is a method to model the uncertainty of a stochastic simulation by modelling both the overall simulation response and the output distribution at each sample point. The output distribution is characterized by dividing it into quantiles, where the division of each quantile is determined by kriging.
This library is our re-implementation of Quantile Kriging as described by Matthew Plumlee & Rui Tuo in their 2014 paper "Building Accurate Emulators for Stochastic Simulations via Quantile Kriging." With computational savings when dealing with replication from the recent paper "Practical heteroskedastic Gaussian process modeling for large simulation experiments " by Binois, M., Gramacy, R., and Ludovski, M. it is now possible to apply Quantile Kriging to a wider class of problems. In addition to fitting the model, other useful tools are provided such as the ability to automatically perform leave-one-out cross validation.
quantKrig : Implements Quantile Kriging
quantPlot : Plots the Quantile output from quantKrig if there is only one input.
newQuants : Generates new quantiles from a quantile object
LOOquants : Generates Leave-One-Out predictions for each location and quantile.
Plots the Quantile output from quantKrig if there is only one input.
QuantPlot(QKResults, X1 = NULL, Y1 = NULL, main = NULL, xlab = NULL, ylab = NULL, colors = NULL)QuantPlot(QKResults, X1 = NULL, Y1 = NULL, main = NULL, xlab = NULL, ylab = NULL, colors = NULL)
QKResults |
Output from the quantKrig function. |
X1 |
X values if ploting the original data in the background |
Y1 |
Y values if ploting the original data in the background |
main |
Plot Title defaults to Fitted Quantiles |
xlab |
Label for x-axis defaults to X |
ylab |
Label for y-axis defaults to Y |
colors |
Customize colors associated with the quantiles |
A ggplot object
X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) Ystar <- rnorm(length(Ystar),Ystar,1) Ystar <- (Ystar - mean(Ystar)) / sd(Ystar) Xstar <- (Xstar - min(Xstar)/ max(Xstar) - min(Xstar)) lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar)X <- seq(0,1,length.out = 20) Y <- cos(5*X) + cos(X) Xstar <- rep(X,each = 100) Ystar <- rep(Y,each = 100) Ystar <- rnorm(length(Ystar),Ystar,1) Ystar <- (Ystar - mean(Ystar)) / sd(Ystar) Xstar <- (Xstar - min(Xstar)/ max(Xstar) - min(Xstar)) lb <- c(0.0001,0.0001) ub <- c(10,10) Qout <- quantKrig(Xstar,Ystar, seq(0.05,0.95, length.out = 7), lower = lb, upper = ub) QuantPlot(Qout, Xstar, Ystar)