# MCMC sampling for Heteroskedastic GP with variance changing in a subset of dimensions Conducts MCMC sampling of hyperparameters and latent noise process `llam` for a hetGP. Separate length scale parameters `theta_lam` and `theta_y` govern the correlation strength of the hidden layer and outer layer respectively. `lam` layer may have a non-zero nugget `g` which governs noise for the latent noise layer. `tau2_y` and `tau2_lam` control the amplitude of the mean and noise process respectively. In Matern covariance, `v` governs smoothness. ```r bhetGP_vdims( x = NULL, y = NULL, reps_list = NULL, nmcmc = 1000, sep = TRUE, inits = NULL, priors = NULL, reps = TRUE, vdims = NULL, cov = c("exp2", "matern", "ARD matern"), v = 2.5, stratergy = c("default", "flat"), vecchia = FALSE, m = min(25, length(y) - 1), ordering = NULL, verb = TRUE, omp_cores = 4, vecchia_var = FALSE ) ``` ## Arguments - `x`: vector or matrix of input locations - `y`: vector of response values - `reps_list`: list object from hetGP::find_reps - `nmcmc`: number of MCMC iterations - `sep`: logical indicating whether to fit isotropic GP (`sep = FALSE`) or seperable GP (`sep = TRUE`) - `inits`: set initial values for hyparameters: `llam`, `theta_y`, `theta_lam`, `g`, `mean_y`, `mean_lam`, `scale_y`, `scale_lam`. Additionally, set initial conditions for tuning: * `theta_check`: logical; if `theta_check = TRUE`, then ensures that theta_lam \> theta_y i.e., decay of correlation for noise process is slower than mean process. * `prof_ll_lam`: logical; if `prof_ll_lam = TRUE`, infers tau2_lam i.e., scale parameter for latent noise process * `noise`: logical; if `noise = TRUE`, infers nugget `g` throught M-H for latent noise process. - `priors`: hyperparameters for priors and proposals (see details) - `reps`: logical; if `reps = TRUE` uses Woodbury inference adjusting for replication of design points and `reps = FALSE` does not use Woodbury inference - `vdims`: vector of dimension index within which noise changes - `cov`: covariance kernel, either Matern, ARD Matern or squared exponential (`"exp2"`) - `v`: Matern smoothness parameter (only used if `cov = "matern"`) - `stratergy`: choose initialization stratergy; "default" uses hetGP for `vecchia = FALSE` settings and sVecchia for `vecchia = TRUE`. See details. - `vecchia`: logical indicating whether to use Vecchia approximation - `m`: size of Vecchia conditioning sets (only used if `vecchia = TRUE`) - `ordering`: optional ordering for Vecchia approximation, must correspond to rows of `x`, defaults to random, is applied to `x` (only used if `vecchia = TRUE`) - `verb`: logical indicating whether to print progress - `omp_cores`: logical; if `vecchia = TRUE`, user may specify the number of cores to use for OpenMP parallelization. Uses min(4, limit) where limit is max openMP cores available on the machine. - `vecchia_var`: logical indicating whether to use Vecchia approximation for noise process ## Returns a list of the S3 class `bhetgp` or `bhetgp_vec` with elements: * `x`: copy of input matrix * `y`: copy of response mean at inputs (x) * `Ylist`: list of all responses observed per location (x) * `A`: number of replicates at each location * `nmcmc`: number of MCMC iterations * `priors`: copy of proposal/priors * `settings`: setting for predictions using `bhetgp` or `bhetgp_vec` object * `theta_y`: vector of MCMC samples for `theta_y` (length scale of mean process) * `theta_lam`: matrix of MCMC samples for `theta_lam` (length scale of latent noise process) * `llam_samples`: matrix of ESS samples for `log lambda` (latent noise process samples) * `g`: vector of MCMC samples for `g` if infered * `tau2`: vector of MAP estimates for `tau2` (scale parameter of mean process) * `tau2_lam`: vector of MAP estimates for `tau2_lam` (scale parameter of latent noise process) * `llik_y`: vector of MVN log likelihood of Y for reach Gibbs iteration * `llik_lam`: vector of MVN log likelihood of `llam` i.e. the latent noise process for reach Gibbs iteration * `x_approx`: conditioning set, NN and ordering for `vecchia = TRUE` * `m`: copy of size of conditioning set for `vecchia = TRUE` * `time`: computation time in seconds ## Details Maps inputs `x` to mean response `y` and noise levels `llam`. Noise levels vary with subset `vdims` dimensions wihtin `x`. Conducts sampling of the latent noise process using Elliptical Slice sampling. Utilizes Metropolis Hastings sampling of the length scale and nugget parameters with proposals and priors controlled by `priors`. `g` for the noise process is set to a specific value, and by default, is not estimated. When `vecchia = TRUE`, all calculations leverage the Vecchia approximation with specified conditioning set size `m`. `tau2_y` is always inferred from likelihood; `tau2_lam` is inferred by default but may be pre-specified and fixed. NOTE on OpenMP: The Vecchia implementation relies on OpenMP parallelization for efficient computation. This function will produce a warning message if the package was installed without OpenMP (this is the default for CRAN packages installed on Apple machines). To set up OpenMP parallelization, download the package source code and install using the gcc/g++ compiler. Proposals for `g` and `theta` follow a uniform sliding window scheme, e.g. `theta_star <- runif(1, l * theta_t / u, u * theta_t / l)`, with defaults `l = 1` and `u = 2` provided in `priors`. To adjust these, set `priors = list(l = new_l, u = new_u)`. Priors on `g`, `theta_y`, and `theta_lam` follow Gamma distributions with shape parameters (`alpha`) and rate parameters (`beta`) controlled within the `priors` list object. Defaults are * `priors$alpha$theta_lam \<- 1.5` * `priors$beta$theta_lam \<- 4` * `priors$alpha$theta_y \<- 1.5` * `priors$beta$theta_y \<- 4` * `priors$alpha$g \<- 1.5` * `priors$beta$g \<- 4` `tau2_y` and `tau2_lam` are not sampled; rather directly inferred under conjugate Inverse Gamma prior with shape (`alpha`) and scale parameters (`beta`) within the `priors` list object * `priors$a$tau2_y \<- 10` * `priors$a$tau2_y \<- 4` * `priors$a$tau2_lam \<- 10` * `priors$a$tau2_lam \<- 4` These priors are designed for `x` scaled to [0, 1] and `y` having mean `mean_y`. These may be adjusted using the `priors` input. Initial values for `theta_y`, `theta_lam`, `llam` may be specified by the user. If no initial values are specified, `stratergy` will determine the initialization method. `stratergy = "default"` leverages mleHetGP for initial values of hyper-parameters if `vecchia = FALSE` and Scaled-Vecchia with Stochastic Kriging (Sk-Vec) hybrid approach if `vecchia = TRUE`. For SK-Vec hybrid approach, code appropriated from [https://github.com/katzfuss-group/scaledVecchia/blob/master/vecchia_scaled.R](https://github.com/katzfuss-group/scaledVecchia/blob/master/vecchia_scaled.R) is used to fit two GPs using the Vecchia approximation. The first for (x, y) pairs, which result in estimated residual sums of squares based on predicted y values. Another GP on (x[, vdims], s) to obtain latent noise estimates which are smoothed. A script is leveraged internally within this package that fits this method. Optionally, choose stratergy = "flat" which which will start at uninformative initial values; `llam` = log(var(y) * 0.1) or specify initial values. The output object of class `bhetgp` or `bhetgp_vec` is designed for use with `trim`, `predict`, and `plot`. ## References Binois, Mickael, Robert B. Gramacy, and Mike Ludkovski. "Practical heteroscedastic Gaussian process modeling for large simulation experiments." Journal of Computational and Graphical Statistics 27.4 (2018): 808-821. Katzfuss, Matthias, Joseph Guinness, and Earl Lawrence. "Scaled Vecchia approximation for fast computer-model emulation." SIAM/ASA Journal on Uncertainty Quantification 10.2 (2022): 537-554. Sauer, Annie Elizabeth. "Deep Gaussian process surrogates for computer experiments." (2023).

bhetGP_vdims function