zvcv {ZVCV} | R Documentation |
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 |
samples |
An |
derivatives |
An |
log_weights |
(optional) A vector of length |
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 |
est_inds |
(optional) A vector of indices for the estimation-only samples. The default when |
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
|
folds |
The number of folds used in k-fold cross-validation for selecting the optimal control variate. Depending on the |
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
: Theintegrand_logged
input stored for reference. -
est_inds
: Theest_inds
input stored for reference. -
polyorder
: Thepolyorder
value used in the final estimate. -
regul_reg
: Theregul_reg
flag used in the final estimate. -
alpha_elnet
: Thealpha_elnet
value used in the final estimate. -
nfolds
: Thenfolds
value used in the final estimate. -
apriori
: Theapriori
vector used in the final estimate. -
intercept
: Theintercept
flag used in the final estimate. -
polyorder_max
: Thepolyorder_max
flag used in the final estimate, if multipleoptions
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