# zvcv: ZV-CV for general expectations In ZVCV: Zero-Variance Control Variates

## Description

The function `zvcv` is used to perform (regularised) ZV-CV given a set of samples, derivatives and function evaluations.

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```zvcv( integrand, samples, derivatives, log_weights, integrand_logged = FALSE, est_inds, options = list(polyorder = 2, regul_reg = TRUE, alpha_elnet = 1, nfolds = 10, apriori = (1:NCOL(samples)), intercept = TRUE, polyorder_max = Inf), folds = 5 ) ```

## Arguments

 `integrand` An N by k matrix of integrands (evaluations of the functions of interest) `samples` An N by d matrix of samples from the target `derivatives` An N by d matrix of derivatives of the log target with respect to the parameters `log_weights` (optional) A vector of length N containing log weights of the samples. The default is equal weights. `integrand_logged` (optional) Sometimes it is better to input the integrand on the logged scale for stability. If the actual integrand is the exponential of `integrand`, then `integrand_logged = TRUE`. Otherwise, the default of `integrand_logged = FALSE` should be used. `est_inds` (optional) A vector of indices for the estimation-only samples. The default when `est_inds` is missing or `NULL` is to perform both estimation of the control variates and evaluation of the integral using all samples. Otherwise, the samples from `est_inds` are used in estimating the control variates and the remainder are used in evaluating the integral. Splitting the indices in this way can be used to reduce bias from adaption and to make computation feasible for very large sample sizes (small `est_inds` is faster), but in general in will increase the variance of the estimator. `options` A list of control variate specifications. This can be a single list containing the elements below (the defaults are used for elements which are not specified). Alternatively, it can be a list of lists containing any or all of the elements below. Where the latter is used, the function `zvcv` automatically selects the best performing option based on cross-validation. `polyorder`: The order of the polynomial, with a default of 2. A value of `Inf` will get the cross-validation method to choose between orders. `regul_reg`: A flag for whether regularised regression is to be used. The default is TRUE, i.e. regularised regression is used. `alpha_elnet`: The alpha parameter for elastic net. The default is 1, which correponds to LASSO. A value of 0 would correspond to ridge regression. `nfolds`: The number of folds used in cross-validation to select lambda for LASSO or elastic net. The default is 10. `apriori`: A vector containing the subset of parameter indices to use in the polynomial. Typically this argument would only be used if the dimension of the problem is very large or if prior information about parameter dependencies is known. The default is to use all parameters 1:d where d is the dimension of the target. In `zvcv`, this is equivalent to using only the relevant columns in `samples` and `derivatives`). `intercept`: A flag for whether the intercept should be estimated or fixed to the empirical mean of the integrand in the estimation set. The default is to include an intercept (`intercept = TRUE`) as this tends to lead to better variance reductions. Note that an `intercept = TRUE` flag may be changed to `intercept = FALSE` within the function if `integrand_logged = TRUE` and a `NaN` is encountered. See South et al. (2018) for further details. `polyorder_max`: The maximum allowable polynomial order. This may be used to prevent memory issues in the case that the polynomial order is selected automatically. A default maximum polynomial order based on the regression design matrix having no more than ten million elements will be selected if the `polyorder` is infinite and in this case a warning will be given. Recall that setting your default R settings to `options(warn=1)` will ensure that you receive these warnings in real time. Optimal polynomial order selection may go to at most this maximum value, or it may stop earlier. `folds` The number of folds used in k-fold cross-validation for selecting the optimal control variate. Depending on the `options`, this may include selection of the optimal polynomial order, regression type and subset of parameters in the polynomial. The default is five.

## Value

A list is returned, containing the following components:

• `expectation`: The estimates of the expectations.

• `num_select`: The number of non-zero coefficients in the polynomial.

• `mse`: The mean square error for the evaluation set.

• `coefs`: The estimated coefficients for the regression (columns are for the different integrands).

• `integrand_logged`: The `integrand_logged` input stored for reference.

• `est_inds`: The `est_inds` input stored for reference.

• `polyorder`: The `polyorder` value used in the final estimate.

• `regul_reg`: The `regul_reg` flag used in the final estimate.

• `alpha_elnet`: The `alpha_elnet` value used in the final estimate.

• `nfolds`: The `nfolds` value used in the final estimate.

• `apriori`: The `apriori` vector used in the final estimate.

• `intercept`: The `intercept` flag used in the final estimate.

• `polyorder_max`: The `polyorder_max` flag used in the final estimate, if multiple `options` are specified.

Leah F. South

## References

Mira, A., Solgi, R., & Imparato, D. (2013). Zero variance Markov chain Monte Carlo for Bayesian estimators. Statistics and Computing, 23(5), 653-662.

South, L. F., Oates, C. J., Mira, A., & Drovandi, C. (2019). Regularised zero variance control variates for high-dimensional variance reduction. https://arxiv.org/abs/1811.05073

See ZVCV and `VDP` for additional examples. See `evidence` for functions which use `zvcv` to estimate the normalising constant of the posterior.
 ``` 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``` ```# An example where ZV-CV can result in zero-variance estimators # Estimating some expectations when theta is bivariate normally distributed with: mymean <- c(-1.5,1.5) mycov <- matrix(c(1,0.5,0.5,2),nrow=2) # Perfect draws from the target distribution (could be replaced with # approximate draws from e.g. MCMC or SMC) N <- 30 require(mvtnorm) set.seed(1) samples <- rmvnorm(N, mean = mymean, sigma = mycov) # derivatives of Gaussian wrt x derivatives <- t( apply(samples,1,function(x) -solve(mycov)%*%(x - mymean)) ) # The integrands are the marginal posterior means of theta, the variances and the # covariance (true values are c(-1.5,1.5,1,2,0.5)) integrand <- cbind(samples[,1],samples[,2],(samples[,1] - mymean)^2, (samples[,2] - mymean)^2, (samples[,1] - mymean)*(samples[,2] - mymean)) # Estimates without ZV-CV (i.e. vanilla Monte Carlo integration) # Vanilla Monte Carlo sprintf("%.15f",colMeans(integrand)) # ZV-CV with fixed specifications # For this example, polyorder = 1 with OLS is exact for the first two integrands and # polyorder = 2 with OLS is exact for the last three integrands # ZV-CV with 2nd order polynomial, OLS and a polynomial based on only x_1. # For diagonal mycov, this would be exact for the first and third expectations. sprintf("%.15f",zvcv(integrand, samples, derivatives, options = list(polyorder = 2, regul_reg = FALSE, apriori = 1))\$expectation) # ZV-CV with 1st order polynomial and OLS (exact for the first two integrands) sprintf("%.15f",zvcv(integrand, samples, derivatives, options = list(polyorder = 1, regul_reg = FALSE))\$expectation) # ZV-CV with 2nd order polynomial and OLS (exact for all) sprintf("%.15f",zvcv(integrand, samples, derivatives, options = list(polyorder = 2, regul_reg = FALSE))\$expectation) # ZV-CV with cross validation myopts <- list(list(polyorder = Inf, regul_reg = FALSE),list(polyorder = Inf, nfolds = 4)) temp <- zvcv(integrand,samples,derivatives,options = myopts, folds = 2) temp\$polyorder # The chosen control variate order temp\$regul_reg # Flag for if the chosen control variate uses regularisation # Cross-val ZV-CV to choose the polynomial order and whether to perform OLS or LASSO sprintf("%.15f",temp\$expectation) # Estimate based on the chosen control variate ```