# GP_deviance: Computes the Deviance of a GP model In GPfit: Gaussian Processes Modeling

## Description

Evaluates the deviance (negative 2*log-likelihood), as defined in Ranjan et al. (2011), however the correlation is reparametrized and can be either power exponential or Matern as discussed in `corr_matrix`.

## Usage

 ```1 2``` ```GP_deviance(beta, X, Y, nug_thres = 20, corr = list(type = "exponential", power = 1.95)) ```

## Arguments

 `beta` a (d x 1) vector of correlation hyper-parameters, as described in `corr_matrix` `X` the (n x d) design matrix `Y` the (n x 1) vector of simulator outputs `nug_thres` a parameter used in computing the nugget. See `GP_fit`. `corr` a list of parameters for the specifing the correlation to be used. See `corr_matrix`.

## Value

the deviance (negative 2 * log-likelihood)

## Author(s)

Blake MacDonald, Hugh Chipman, Pritam Ranjan

## 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.

`corr_matrix` for computing the correlation matrix;
`GP_fit` for estimating the parameters of the GP model.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52``` ```## 1D Example 1 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) beta = rnorm(1) GP_deviance(beta,x,y) ## 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) beta = rnorm(1) GP_deviance(beta, x, y, corr = list(type = "matern", nu = 5/2)) ## 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) beta = rnorm(2) GP_deviance(beta, x, y) ```