# DTV: Compute Distance Transform Variability In pGPx: Pseudo-Realizations for Gaussian Process Excursions

## Description

Compute the expected L^2 distance between the average distance transform and the set realizations. If the input is the actual values of the gaussian process, compute also the random sets.

## Usage

 `1` ```DTV(rand.set, threshold, nsim, n.int.points) ```

## Arguments

 `rand.set` a matrix of size `n.int.points`x`nsim` containing the excursion set realizations stored as long vectors. For example the excursion set obtained from the result of `simulate_and_interpolate`. `threshold` threshold value `nsim` number of simulations of the excursion set `n.int.points` total length of the excursion set discretization. The size of the image is `sqrt(n.int.points)`.

## Value

A list containing

• `variance:`Value of the distance transform variability. The integral of `dvar` over the spatial domain.

• `dbar:`empirical distance average transform 1/N ∑_{i=1}^N d(x,Γ_i), a matrix of size `n.int.points` x `n.int.points`

• `dvar:`empirical variance of distance transform 1/N ∑_{i=1}^N (d(x,Γ_i) - dbar)^2, a matrix of size `n.int.points` x `n.int.points`

• `alldt:`distance transforms for all realizations, a matrix of size `n.int.points` x `nsim`

• `naTot:`Total number of infinite distance transform values. These are returned in realizations where there is no excursion.

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

Felzenszwalb, P. F. and Huttenlocher, D. P. (2012). Distance Transforms of Sampled Functions. Theory of Computing, 8(19):415-428.

## Examples

 ``` 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``` ```### Simulate and interpolate for a 2d example if (!requireNamespace("DiceKriging", quietly = TRUE)) { stop("DiceKriging needed for this example to work. Please install it.", call. = FALSE) } if (!requireNamespace("DiceDesign", quietly = TRUE)) { stop("DiceDesign needed for this example to work. Please install it.", call. = FALSE) } # Define the function g=function(x){ return(-DiceKriging::branin(x)) } d=2 # Fit OK km model design<-DiceDesign::maximinESE_LHS(design = DiceDesign::lhsDesign(n=50, dimension = 2, seed=42)\$design)\$design colnames(design)<-c("x1","x2") observations<-apply(X = design,MARGIN = 1,FUN = g) kmModel<-DiceKriging::km(formula = ~1,design = design,response = observations, covtype = "matern3_2",control=list(trace=FALSE)) # Get simulation points # Here they are not optimized, you can use optim_dist_measure to find optimized points simu_points <- DiceDesign::maximinSA_LHS(DiceDesign::lhsDesign(n=100, dimension = d, seed=1)\$design)\$design # obtain nsims posterior realization at simu_points nsims <- 30 nn_data<-expand.grid(seq(0,1,,50),seq(0,1,,50)) nn_data<-data.frame(nn_data) colnames(nn_data)<-colnames(kmModel@X) approx.simu <- simulate_and_interpolate(object=kmModel, nsim = nsims, simupoints = simu_points, interpolatepoints = as.matrix(nn_data), nugget.sim = 0, type = "UK") Dvar<- DTV(rand.set = approx.simu,threshold = -10, nsim = nsims,n.int.points = 50^2) image(matrix(Dvar\$dbar,ncol=50),col=grey.colors(20),main="average distance transform") image(matrix(Dvar\$dvar,ncol=50),col=grey.colors(20),main="variance of distance transform") points(design,pch=17) ```

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