GPfit (version 1.0-0)

predict.GP: Model Predictions from GPfit

Description

Computes the regularized predicted response \(\hat{y}_{\delta_{lb},M}(x)\) and the mean squared error \(s^2_{\delta_{lb},M}(x)\) for a new set of inputs using the fitted GP model.

Usage

# S3 method for GP
predict(object, xnew = object$X,M=1,…)

Arguments

object

a class GP object estimated by GP_fit

xnew

the (n_new x d) design matrix of test points where model predictions and MSEs are desired

M

the number of iterations. See `Details'

for compatibility with generic method predict

Value

Returns a list containing the predicted values (Y_hat), the mean squared errors of the predictions (MSE), and a matrix (complete_data) containing xnew, Y_hat, and MSE

Details

The value of M determines the number of iterations (or terms) in approximating \(R^{-1} \approx R^{-1}_{\delta_{lb},M}\). The iterative use of the nugget \(\delta_{lb}\), as outlined in Ranjan et al. (2011), is used in calculating \(\hat{y}_{\delta_{lb},M}(x)\) and \(s^2_{\delta_{lb},M}(x)\), where \(R_{\delta,M}^{-1} = \sum_{t = 1}^{M} \delta^{t - 1}(R+\delta I)^{-t}\).

References

Ranjan, P., Haynes, R., and Karsten, R. (2011). A Computationally Stable Approach to Gaussian Process Interpolation of Deterministic Computer Simulation Data, Technometrics, 53(4), 366 - 378.

See Also

GP_fit for estimating the parameters of the GP model; plot.GP for plotting the predicted and error surfaces.

Examples

Run this code
# NOT RUN {
## 1D Example
n = 5; d = 1;
computer_simulator <- function(x){
x = 2*x+0.5;
y = sin(10*pi*x)/(2*x) + (x-1)^4;
return(y)
}
set.seed(3);
library(lhs);
x = maximinLHS(n,d);
y = computer_simulator(x);
xvec <- seq(from=0,to=1,length.out=10);
GPmodel = GP_fit(x,y);
GPprediction = predict.GP(GPmodel,xvec);
yhat = GPprediction$Y_hat;
mse = GPprediction$MSE;
completedata = GPprediction$complete_data;
completedata;


## 1D Example 2
n = 7; d = 1;
computer_simulator <- function(x) {
y <- log(x+0.1)+sin(5*pi*x);
return(y)
}
set.seed(1)
library(lhs);
x = maximinLHS(n,d)
y = computer_simulator(x)
xvec <- seq(from=0,to=1,length.out=10);
GPmodel = GP_fit(x,y)
GPprediction = predict.GP(GPmodel,xvec);
yhat = GPprediction$Y_hat;
mse = GPprediction$MSE;
completedata = GPprediction$complete_data;


## 2D Example: GoldPrice Function
computer_simulator <- function(x) {
x1=4*x[,1] - 2; x2=4*x[,2] - 2;
t1 = 1 + (x1 + x2 + 1)^2*(19 - 14*x1 + 3*x1^2 - 14*x2 + 
6*x1*x2 + 3*x2^2);
t2 = 30 + (2*x1 -3*x2)^2*(18 - 32*x1 + 12*x1^2 + 48*x2 - 
36*x1*x2 + 27*x2^2);
y = t1*t2;
return(y)
}
n = 10; d = 2;
set.seed(1);
library(lhs);
x = maximinLHS(n,d); 
y = computer_simulator(x);
GPmodel = GP_fit(x,y);

xvector = seq(from=0,to=1,length.out=10);
xvec = expand.grid(x = xvector, y=xvector);
xvec = as.matrix(xvec);
GPprediction = predict.GP(GPmodel,xvec);
yhat = GPprediction$Y_hat;
mse = GPprediction$MSE;
completedata = GPprediction$complete_data;
# }

Run the code above in your browser using DataCamp Workspace