zvcv {ZVCV}

ZV-CV for general expectations

Description

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

Usage

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.

Author(s)

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 Also

See ZVCV and VDP for additional examples. See evidence for functions which use zvcv to estimate the normalising constant of the posterior.

Examples

# 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[1])^2,
    (samples[,2] - mymean[2])^2, (samples[,1] - mymean[1])*(samples[,2] - mymean[2]))

# 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