# R/integrand_edm_crit.R In pGPx: Pseudo-Realizations for Gaussian Process Excursions

```#' @title Integrand of the distance in measure criterion
#' @description Computes the integrand of the distance in measure criterion.
# Input:
#' @param x vector of dimension \eqn{d} representing the \eqn{ith} point where to compute the criterion
#' @param E matrix of dimension \eqn{d*(i-1)} containing the previously optimized simulation points
#' @param model km model
#' @param Thresh threshold selected for excursion set
#' @param batchsize number of simulation points
#' @param alpha value of Vorob'ev threshold
#' @param predE list containing the posterior mean and standard deviation at E
#' @param predx list containing the posterior mean and standard deviation at x
#' @param precalc.data list result of \link[KrigInv]{precomputeUpdateData} with \code{model} and \code{x}.
#' @return the value of the integrand at \eqn{x}
#' @references Azzimonti D. F., Bect J., Chevalier C. and Ginsbourger D. (2016). Quantifying uncertainties on excursion sets under a Gaussian random field prior. SIAM/ASA Journal on Uncertainty Quantification, 4(1):850–874.
#'
#' Azzimonti, D. (2016). Contributions to Bayesian set estimation relying on random field priors. PhD thesis, University of Bern.
#' @export
integrand_edm_crit = function (x, E, model, Thresh, batchsize, alpha,predE,
predx = NULL,precalc.data=NULL)
{

x<-t(x)

if(is.null(predx)){
predx <- KrigInv::predict_nobias_km(object=model,newdata=x,type="UK",se.compute=TRUE)
}
evalpoint.oldmean <- predx\$mean ;
evalpoint.oldsd <- ifelse(predx\$sd<.Machine\$double.eps,0,predx\$sd);

if(is.null(precalc.data))
precalc.data <- KrigInv::precomputeUpdateData(model,x)

d <- model@d

X.new <- matrix(t(E), nrow = d)

if (!identical(colnames(x), colnames(model@X)))
colnames(x) <- colnames(model@X)

X.new <- t(X.new)
mk <- predE\$mean
sk <- predE\$sd
newXvar <- sk * sk
F.newdata <- predE\$F.newdata
c.newdata <- predE\$c
Sigma.r <- predE\$cov
kn = KrigInv::computeQuickKrigcov(model, x, X.new,
precalc.data, F.newdata, c.newdata)
krig2 <- KrigInv::predict_update_km_parallel(newXmean = mk, newXvar = newXvar,
newXvalue = mk, Sigma.r = Sigma.r, newdata.oldmean = evalpoint.oldmean,
newdata.oldsd = evalpoint.oldsd, kn = kn)
if (!is.null(krig2\$error)){
current.vorob <- pnorm((evalpoint.oldmean-Thresh)/evalpoint.oldsd)
return(current.vorob)
}

gamma_n <- rowSums(krig2\$lambda*kn)
gamma_n <- ifelse(gamma_n <.Machine\$double.eps,0, gamma_n)

arg1 <- as.numeric((evalpoint.oldmean - Thresh)/evalpoint.oldsd)
arg2 <- as.numeric((Thresh-evalpoint.oldmean)/sqrt(gamma_n))
arg3 <- -sqrt(gamma_n)/evalpoint.oldsd

arg1[arg1 == -Inf] <- -1000
arg1[arg1 == Inf] <- 1000
arg1[is.nan(arg1)] <- 1000

arg2[arg2 == -Inf] <- -1000
arg2[arg2 == Inf] <- 1000
arg2[is.nan(arg2)] <- -1000*sign(arg1)

arg3[arg3 <= -1] <- -1
arg3[arg3 >= 1] <- 1
arg3[is.nan(arg3)] <- 0

term1 <- pbivnorm::pbivnorm(arg1, arg2, arg3)
term2 <- pbivnorm::pbivnorm(-arg1, -arg2, arg3)
crit <- term1 + term2

return(crit)
}
```

## Try the pGPx package in your browser

Any scripts or data that you put into this service are public.

pGPx documentation built on May 2, 2019, 3:28 a.m.