Package 'ZVCV'

Approximate semi-exact control functionals (aSECF)

Description

This function performs approximate semi-exact control functionals as described in South et al (2020). It uses a nystrom approximation and conjugate gradient to speed up SECF. This is faster than SECF for large $N$. If you would like to choose between different kernels using cross-validation, then you can use aSECF_crossval.

Usage

aSECF(
  integrands,
  samples,
  derivatives,
  polyorder = NULL,
  steinOrder = NULL,
  kernel_function = NULL,
  sigma = NULL,
  K0 = NULL,
  nystrom_inds = NULL,
  est_inds = NULL,
  apriori = NULL,
  conjugate_gradient = TRUE,
  reltol = 0.01,
  diagnostics = FALSE
)

Arguments

integrands	An $N$ by $k$ matrix of integrands (evaluations of the function 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
polyorder	(optional) The order of the polynomial to be used in the parametric component, with a default of $1$. We recommend keeping this value low (e.g. only 1-2).
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma	(optional) The tuning parameters of the specified kernel. This involves a single length-scale parameter in "gaussian" and "RQ", a length-scale and a smoothness parameter in "matern" and two parameters in "product" and "prodsim". See below for further details.
K0	(optional) The kernel matrix. One can specify either this or all of sigma, steinOrder and kernel_function. The former involves pre-computing the kernel matrix using K0_fn and is more efficient when using multiple estimators out of CF, SECF and aSECF or when using the cross-validation functions.
nystrom_inds	(optional) The sample indices to be used in the Nystrom approximation.
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.
apriori	(optional) 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.
conjugate_gradient	(optional) A flag for whether to perform conjugate gradient to further speed up the nystrom approximation (the default is true).
reltol	(optional) The relative tolerance for choosing when the stop conjugate gradient iterations (the default is 1e-02). using squareNorm, as long as the nystrom_inds are NULL.
diagnostics	(optional) A flag for whether to return the necessary outputs for plotting or estimating using the fitted model. The default is false since this requires some additional computation when est_inds is NULL.

Value

A list with the following elements:

• expectation: The estimate(s) of the ($k$) expectations(s).

• cond_no: (Only if conjugate_gradient = TRUE) The condition number of the matrix being solved using conjugate gradient.

• iter: (Only if conjugate_gradient = TRUE) The number of conjugate gradient iterations

• f_true: (Only if est_inds is not NULL) The integrands for the evaluation set. This should be the same as integrands[setdiff(1:N,est_inds),].

• f_hat: (Only if est_inds is not NULL) The fitted values for the integrands in the evaluation set. This can be used to help assess the performance of the Gaussian process model.

• a: (Only if diagnostics = TRUE) The value of $a$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• b: (Only if diagnostics = TRUE) The value of $b$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• ny_inds: (Only if diagnostics = TRUE) The indices of the samples used in the nystrom approximation (this will match nystrom_inds if this argument was not NULL).

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

See Also

See ZVCV for examples and related functions. See aSECF_crossval for a function to choose between different kernels for this estimator.

---

Approximate semi-exact control functionals (aSECF) with cross-validation

Description

This function chooses between a list of kernel tuning parameters (sigma_list) or a list of K0 matrices (K0_list) for the approximate semi-exact control functionals method described in South et al (2020). The latter requires calculating and storing kernel matrices using K0_fn but it is more flexible because it can be used to choose the Stein operator order and the kernel function, in addition to its parameters. It is also faster to pre-specify K0_fn. For estimation with fixed kernel parameters, use aSECF.

Usage

aSECF_crossval(
  integrands,
  samples,
  derivatives,
  polyorder = NULL,
  steinOrder = NULL,
  kernel_function = NULL,
  sigma_list = NULL,
  est_inds = NULL,
  apriori = NULL,
  num_nystrom = NULL,
  conjugate_gradient = TRUE,
  reltol = 0.01,
  folds = NULL,
  diagnostics = FALSE
)

Arguments

integrands	An $N$ by $k$ matrix of integrands (evaluations of the function 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
polyorder	(optional) The order of the polynomial to be used in the parametric component, with a default of $1$. We recommend keeping this value low (e.g. only 1-2).
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma_list	(optional between this and K0_list) A list of tuning parameters for the specified kernel. This involves a list of single length-scale parameter in "gaussian" and "RQ", a list of vectors containing length-scale and smoothness parameters in "matern" and a list of vectors of the two parameters in "product" and "prodsim". See below for further details. When sigma_list is specified and not K0_list, the $K0$ matrix is computed twice for each selected tuning parameter.
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.
apriori	(optional) 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.
num_nystrom	(optional) The number of samples to use in the Nystrom approximation, with a default of ceiling(sqrt(N)). The nystrom indices cannot be passed in here because of the way the cross-validation has been set up.
conjugate_gradient	(optional) A flag for whether to perform conjugate gradient to further speed up the nystrom approximation (the default is true).
reltol	(optional) The relative tolerance for choosing when the stop conjugate gradient iterations (the default is 1e-02). using squareNorm, as long as the nystrom_inds are NULL.
folds	(optional) The number of folds for cross-validation. The default is five.
diagnostics	(optional) A flag for whether to return the necessary outputs for plotting or estimating using the fitted model. The default is false since this requires some additional computation when est_inds is NULL.

Value

A list with the following elements:

• expectation: The estimate(s) of the ($k$) expectations(s).

• mse: A matrix of the cross-validation mean square prediction errors. The number of columns is the number of tuning options given and the number of rows is $k$, the number of integrands of interest.

• optinds: The optimal indices from the list for each expectation.

• cond_no: (Only if conjugate_gradient = TRUE) The condition number of the matrix being solved using conjugate gradient.

• iter: (Only if conjugate_gradient = TRUE) The number of conjugate gradient iterations

• f_true: (Only if est_inds is not NULL) The integrands for the evaluation set. This should be the same as integrands[setdiff(1:N,est_inds),].

• f_hat: (Only if est_inds is not NULL) The fitted values for the integrands in the evaluation set. This can be used to help assess the performance of the Gaussian process model.

• a: (Only if diagnostics = TRUE) The value of $a$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• b: (Only if diagnostics = TRUE) The value of $b$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• ny_inds: (Only if diagnostics = TRUE) The indices of the samples used in the nystrom approximation (this will match nystrom_inds if this argument was not NULL).

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

See Also

See ZVCV for examples and related functions. See aSECF_crossval for a function to choose between different kernels for this estimator.

---

Control functionals (CF)

Description

This function performs control functionals as described in Oates et al (2017). To choose between different kernels using cross-validation, use CF_crossval.

Usage

CF(
  integrands,
  samples,
  derivatives,
  steinOrder = NULL,
  kernel_function = NULL,
  sigma = NULL,
  K0 = NULL,
  est_inds = NULL,
  one_in_denom = FALSE,
  diagnostics = FALSE
)

Arguments

integrands	An $N$ by $k$ matrix of integrands (evaluations of the function 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
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma	(optional) The tuning parameters of the specified kernel. This involves a single length-scale parameter in "gaussian" and "RQ", a length-scale and a smoothness parameter in "matern" and two parameters in "product" and "prodsim". See below for further details.
K0	(optional) The kernel matrix. One can specify either this or all of sigma, steinOrder and kernel_function. The former involves pre-computing the kernel matrix using K0_fn and is more efficient when using multiple estimators out of CF, SECF and aSECF or when using the cross-validation functions.
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.
one_in_denom	(optional) Whether or not to include a $1 +$ in the denominator of the control functionals estimator, as in equation 2 on p703 of Oates et al (2017). The $1 +$ in the denominator is an arbitrary choice so we set it to zero by default.
diagnostics	(optional) A flag for whether to return the necessary outputs for plotting or estimating using the fitted model. The default is false since this requires some additional computation when est_inds is NULL.

Value

A list with the following elements:

• expectation: The estimate(s) of the ($k$) expectation(s).

• f_true: (Only if est_inds is not NULL) The integrands for the evaluation set. This should be the same as integrands[setdiff(1:N,est_inds),].

• f_hat: (Only if est_inds is not NULL) The fitted values for the integrands in the evaluation set. This can be used to help assess the performance of the Gaussian process model.

• a: (Only if diagnostics = TRUE) The value of $a$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + 1*b$ for heldout K0 and estimators using heldout samples are of the form $mean(f - f_hat) + b$.

• b: (Only if diagnostics = TRUE) The value of $b$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + 1*b$ for heldout K0 and estimators using heldout samples are of the form $mean(f - f_hat) + b$.

• ksd: (Only if diagnostics = TRUE) An estimated kernel Stein discrepancy based on the fitted model that can be used for diagnostic purposes. See South et al (2020) for further details.

• bound_const: (Only if diagnostics = TRUE and est_inds=NULL) This is such that the absolute error for the estimator should be less than $ksd \times bound_const$.

Warning

Solving the linear system in CF has $O(N^3)$ complexity and is therefore not suited to large $N$. Using $est_inds$ will instead have an $O(N_0^3)$ cost in solving the linear system and an $O((N-N_0)^2)$ cost in handling the remaining samples, where $N_0$ is the length of $est_inds$. This can be much cheaper for large $N$.

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

Oates, C. J., Girolami, M. & Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3), 695-718.

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

See Also

See ZVCV for examples and related functions. See CF_crossval for a function to choose between different kernels for this estimator.

---

Control functionals (CF) with cross-validation

Description

This function chooses between a list of kernel tuning parameters (sigma_list) or a list of K0 matrices (K0_list) for the control functionals method described in Oates et al (2017). The latter requires calculating and storing kernel matrices using K0_fn but it is more flexible because it can be used to choose the Stein operator order and the kernel function, in addition to its parameters. It is also faster to pre-specify K0_fn. For estimation with fixed kernel parameters, use CF.

Usage

CF_crossval(
  integrands,
  samples,
  derivatives,
  steinOrder = NULL,
  kernel_function = NULL,
  sigma_list = NULL,
  K0_list = NULL,
  est_inds = NULL,
  log_weights = NULL,
  one_in_denom = FALSE,
  folds = NULL,
  diagnostics = FALSE
)

Arguments

integrands	An $N$ by $k$ matrix of integrands (evaluations of the function 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
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma_list	(optional between this and K0_list) A list of tuning parameters for the specified kernel. This involves a list of single length-scale parameter in "gaussian" and "RQ", a list of vectors containing length-scale and smoothness parameters in "matern" and a list of vectors of the two parameters in "product" and "prodsim". See below for further details. When sigma_list is specified and not K0_list, the $K0$ matrix is computed twice for each selected tuning parameter.
K0_list	(optional between this and sigma_list) A list of kernel matrices, which can be calculated using K0_fn.
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.
log_weights	(optional) A vector of length $N$ containing the logged weights of the samples. The default is equal weights. The weights are only used in estimating the cross-validation error. This method is not implemented for the case where est_inds is specified becausing specifying est_inds typically indicates a desire for an unbiased estimator and using self-normalised importance weights introduces bias.
one_in_denom	(optional) Whether or not to include a $1 +$ in the denominator of the control functionals estimator, as in equation 2 on p703 of Oates et al (2017). The $1 +$ in the denominator is an arbitrary choice so we set it to zero by default.
folds	(optional) The number of folds for cross-validation. The default is five.
diagnostics	(optional) A flag for whether to return the necessary outputs for plotting or estimating using the fitted model. The default is false since this requires some additional computation when est_inds is NULL.

Value

A list with the following elements:

• expectation: The estimate(s) of the ($k$) expectation(s).

• mse: A matrix of the cross-validation mean square prediction errors. The number of columns is the number of tuning options given and the number of rows is $k$, the number of integrands of interest.

• optinds: The optimal indices from the list for each expectation.

• f_true: (Only if est_inds is not NULL) The integrands for the evaluation set. This should be the same as integrands[setdiff(1:N,est_inds),].

• f_hat: (Only if est_inds is not NULL) The fitted values for the integrands in the evaluation set. This can be used to help assess the performance of the Gaussian process model.

• a: (Only if diagnostics = TRUE) The value of $a$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + 1*b$ for heldout K0 and estimators using heldout samples are of the form $mean(f - f_hat) + b$.

• b: (Only if diagnostics = TRUE) The value of $b$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + 1*b$ for heldout K0 and estimators using heldout samples are of the form $mean(f - f_hat) + b$.

• ksd: (Only if diagnostics = TRUE) An estimated kernel Stein discrepancy based on the fitted model that can be used for diagnostic purposes. See South et al (2020) for further details.

• bound_const: (Only if diagnostics = TRUE and est_inds=NULL) This is such that the absolute error for the estimator should be less than $ksd \times bound_const$.

Warning

Solving the linear system in CF has $O(N^3)$ complexity and is therefore not suited to large $N$. Using $est_inds$ will instead have an $O(N_0^3)$ cost in solving the linear system and an $O((N-N_0)^2)$ cost in handling the remaining samples, where $N_0$ is the length of $est_inds$. This can be much cheaper for large $N$.

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

Oates, C. J., Girolami, M. & Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3), 695-718.

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

See Also

See ZVCV for examples and related functions. See CF for a function to perform control functionals with fixed kernel specifications.

---

Evidence estimation with ZV-CV

Description

The functions evidence_CTI and evidence_CTI_CF can be used to improve upon the thermodynamic integration (TI) estimate of the normalising constant with ZV-CV and CF, respectively. The functions evidence_SMC and evidence_SMC_CF do the same thing for the sequential Monte Carlo (SMC) normalising constant identity.

Usage

evidence_CTI(
  samples,
  loglike,
  der_loglike,
  der_logprior,
  temperatures,
  temperatures_all,
  most_recent,
  est_inds,
  options,
  folds = 5
)

evidence_CTI_CF(
  samples,
  loglike,
  der_loglike,
  der_logprior,
  temperatures,
  temperatures_all,
  most_recent,
  est_inds,
  steinOrder,
  kernel_function,
  sigma_list,
  folds = 5
)

evidence_SMC(
  samples,
  loglike,
  der_loglike,
  der_logprior,
  temperatures,
  temperatures_all,
  most_recent,
  est_inds,
  options,
  folds = 5
)

evidence_SMC_CF(
  samples,
  loglike,
  der_loglike,
  der_logprior,
  temperatures,
  temperatures_all,
  most_recent,
  est_inds,
  steinOrder,
  kernel_function,
  sigma_list,
  folds = 5
)

Arguments

samples	An $N$ by $d$ by $T$ matrix of samples from the $T$ power posteriors, where $N$ is the number of samples and $d$ is the dimension of the target distribution
loglike	An $N$ by $T$ matrix of log likelihood values corresponding to samples
der_loglike	An $N$ by $d$ by $T$ matrix of the derivatives of the log likelihood with respect to the parameters, with parameter values corresponding to samples
der_logprior	An $N$ by $d$ by $T$ matrix of the derivatives of the log prior with respect to the parameters, with parameter values corresponding to samples
temperatures	A vector of length $T$ of temperatures for the power posterior temperatures
temperatures_all	An adjusted vector of length $tau$ of temperatures. Better performance should be obtained with a more conservative temperature schedule. See Expand_Temperatures for a function to adjust the temperatures.
most_recent	A vector of length $tau$ which gives the indices in the original temperatures that the new temperatures correspond to.
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 for ZV-CV. 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.
folds	The number of folds used in k-fold cross-validation for selecting the optimal control variate. For ZV-CV, this may include selection of the optimal polynomial order, regression type and subset of parameters depending on options. For CF, this includes the selection of the optimal tuning parameters in sigma_list. The default is five.
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma_list	(optional between this and K0_list) A list of tuning parameters for the specified kernel. This involves a list of single length-scale parameter in "gaussian" and "RQ", a list of vectors containing length-scale and smoothness parameters in "matern" and a list of vectors of the two parameters in "product" and "prodsim". See below for further details. When sigma_list is specified and not K0_list, the $K0$ matrix is computed twice for each selected tuning parameter.

Value

The function evidence_CTI returns a list, containing the following components:

• log_evidence_PS1: The 1st order quadrature estimate for the log normalising constant

• log_evidence_PS2: The 2nd order quadrature estimate for the log normalising constant

• regression_LL: The set of $tau$ zvcv type returns for the 1st order quadrature expectations

• regression_vLL: The set of $tau$ zvcv type returns for the 2nd order quadrature expectations

The function evidence_CTI_CF returns a list, containing the following components:

• log_evidence_PS1: The 1st order quadrature estimate for the log normalising constant

• log_evidence_PS2: The 2nd order quadrature estimate for the log normalising constant

• regression_LL: The set of $tau$ CF_crossval type returns for the 1st order quadrature expectations

• regression_vLL: The set of $tau$ CF_crossval type returns for the 2nd order quadrature expectations

• selected_LL_CF: The set of $tau$ selected tuning parameters from sigma_list for the 1st order quadrature expectations.

• selected_vLL_CF: The set of $tau$ selected tuning parameters from sigma_list for the 2nd order quadrature expectations.

The function evidence_SMC returns a list, containing the following components:

• log_evidence: The logged SMC estimate for the normalising constant

• regression_SMC: The set of $tau$ zvcv type returns for the expectations

The function evidence_SMC_CF returns a list, containing the following components:

• log_evidence: The logged SMC estimate for the normalising constant

• regression_SMC: The set of $tau$ CF_crossval type returns for the expectations

• selected_CF: The set of $tau$ selected tuning parameters from sigma_list for the expectations

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

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 an example at VDP and see ZVCV for more package details. See Expand_Temperatures for a function that can be used to find stricter (or less stricter) temperature schedules based on the conditional effective sample size.

---

Adjusting the temperature schedule

Description

This function is used to adjust the temperature schedule so that it is more (or less) strict than the original.

Usage

Expand_Temperatures(
  temperatures,
  loglike,
  rho,
  bisec_tol = .Machine$double.eps^0.25
)

Arguments

temperatures	A vector of length $T$ temperatures for the power posterior temperatures.
loglike	An $N$ by $T$ matrix of log likelihood values corresponding to the samples.
rho	The tolerance for the new temperatures. Temperatures are selected so that the conditional effective sample size (CESS) at each temperature is $\rho*N$ where $N$ is the population size.
bisec_tol	The tolerance for the bisection method used in selecting temperatures. The default is .Machine$double.eps^0.25

Value

A list is returned, containing the following components:

• temperatures_all: The new set of temperatures of length $tau$.

• relevant_samples: A vector of length $tau$ containing indices to show which particle sets the new temperatures are based on.

• logw: An $N$ by $tau$ matrix of log normalised weights of the particles

Author(s)

Leah F. South

References

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 evidence for functions to estimate the evidence, VDP for an example and ZVCV for more package details.

---

ZV-CV design matrix

Description

The function getX is used to get the matrix of covariates for the regression based on a specified polynomial order.

Usage

getX(samples, derivatives, polyorder)

Arguments

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
polyorder	The order of the polynomial.

Value

The design matrix for the regression (except for the column of 1's for the intercept).

See Also

Phi_fn for a very similar function for use in semi-exact control functionals. The function Phi_fn essentially gets the same matrix but with a column of ones added.

---

Kernel matrix calculation

Description

This function calculates the full $K_0$ matrix, which is a first or second order Stein operator applied to a standard kernel. The output of this function can be used as an argument to CF, CF_crossval, SECF, SECF_crossval, aSECF and aSECF_crossval. The kernel matrix is automatically computed in all of the above methods, but it is faster to calculate in advance when using more than one of the above functions and when using any of the crossval functions.

Usage

K0_fn(
  samples,
  derivatives,
  sigma,
  steinOrder,
  kernel_function,
  Z = NULL,
  nystrom_inds = NULL
)

Arguments

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
sigma	The tuning parameters of the specified kernel. This involves a single length-scale parameter in "gaussian" and "RQ", a length-scale and a smoothness parameter in "matern" and two parameters in "product" and "prodsim". See below for further details.
steinOrder	This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
Z	(optional) An $N$ by $N$ (or $N$ by $m$ where $m$ is the length of nystrom_inds). This can be calculated using squareNorm.
nystrom_inds	(optional) The sample indices to be used in the Nystrom approximation (for when using aSECF).

Value

An $N$ by $N$ kernel matrix (or $N$ by $m$ where $m$ is the length of nystrom_inds).

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

Oates, C. J., Girolami, M. & Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3), 695-718.

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

---

Stable log sum of exponential calculations

Description

The function logsumexp is used for stable computation of log(sum(exp(x))), which is useful when summing weights for example.

Usage

logsumexp(x)

Arguments

x	The values for which you want to compute log(sum(exp(x)))

Value

The stable result of log(sum(exp(x)))

See Also

See ZVCV for more package details.

---

Median heuristic

Description

This function calculates the median heuristic for use in e.g. the Gaussian, Matern and rational quadratic kernels.

Usage

medianTune(samples, Z = NULL)

Arguments

samples	An $N$ by $d$ matrix of samples from the target
Z	(optional) An NxN matrix of square norms, which can be calculated using squareNorm, as long as the nystrom_inds are NULL.

Value

The median heuristic, which can then be used as the length-scale parameter in the Gaussian, Matern and rational quadratic kernels

Author(s)

Leah F. South

References

Garreau, D., Jitkrittum, W. and Kanagawa, M. (2017). Large sample analysis of the median heuristic. https://arxiv.org/abs/1707.07269

See Also

See medianTune and K0_fn for functions which use this.

---

Nearest symmetric positive definite matrix

Description

This function finds the nearest symmetric positive definite matrix to the given matrix. It is used throughout the package to handle numerical issues in matrix inverses and cholesky decompositions.

Usage

nearPD(K0)

Arguments

K0	A square matrix

Value

The closest symmetric positive definite matrix to K0.

Author(s)

Adapted from Matlab code by John D'Errico

References

Higham, N. J. (1988). Computing a nearest symmetric positive semidefinite matrix. Linear Algebra and its Applications, 103, 103-118.

D'Errico, J. (2013). nearestSPD Matlab function. https://uk.mathworks.com/matlabcentral/fileexchange/42885-nearestspd.

---

Phi matrix calculation

Description

This function calculates the $\Phi$ matrix, which is a second order Stein operator applied to a polynomial. See South et al (2020) for further details. This function is not required for estimation but may be useful when evaluation samples are not initially available since estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$ where $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices.

Usage

Phi_fn(samples, derivatives, polyorder = NULL, apriori = NULL)

Arguments

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
polyorder	(optional) The order of the polynomial to be used in the parametric component, with a default of $1$. We recommend keeping this value low (e.g. only 1-2).
apriori	(optional) 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.

Value

An $N$ by $Q$ matrix (where Q is determined by the polynomial order and the apriori).

Author(s)

Leah F. South

References

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

---

Semi-exact control functionals (SECF)

Description

This function performs semi-exact control functionals as described in South et al (2020). To choose between different kernels using cross-validation, use SECF_crossval.

Usage

SECF(
  integrands,
  samples,
  derivatives,
  polyorder = NULL,
  steinOrder = NULL,
  kernel_function = NULL,
  sigma = NULL,
  K0 = NULL,
  est_inds = NULL,
  apriori = NULL,
  diagnostics = FALSE
)

Arguments

integrands	An $N$ by $k$ matrix of integrands (evaluations of the function 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
polyorder	(optional) The order of the polynomial to be used in the parametric component, with a default of $1$. We recommend keeping this value low (e.g. only 1-2).
steinOrder	(optional) This is the order of the Stein operator. The default is 1 in the control functionals paper (Oates et al, 2017) and 2 in the semi-exact control functionals paper (South et al, 2020). The following values are currently available: 1 for all kernels and 2 for "gaussian", "matern" and "RQ". See below for further details.
kernel_function	(optional) Choose between "gaussian", "matern", "RQ", "product" or "prodsim". See below for further details.
sigma	(optional) The tuning parameters of the specified kernel. This involves a single length-scale parameter in "gaussian" and "RQ", a length-scale and a smoothness parameter in "matern" and two parameters in "product" and "prodsim". See below for further details.
K0	(optional) The kernel matrix. One can specify either this or all of sigma, steinOrder and kernel_function. The former involves pre-computing the kernel matrix using K0_fn and is more efficient when using multiple estimators out of CF, SECF and aSECF or when using the cross-validation functions.
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.
apriori	(optional) 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.
diagnostics	(optional) A flag for whether to return the necessary outputs for plotting or estimating using the fitted model. The default is false since this requires some additional computation when est_inds is NULL.

Value

A list with the following elements:

• expectation: The estimate(s) of the ($k$) expectation(s).

• f_true: (Only if est_inds is not NULL) The integrands for the evaluation set. This should be the same as integrands[setdiff(1:N,est_inds),].

• f_hat: (Only if est_inds is not NULL) The fitted values for the integrands in the evaluation set. This can be used to help assess the performance of the Gaussian process model.

• a: (Only if diagnostics = TRUE) The value of $a$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• b: (Only if diagnostics = TRUE) The value of $b$ as described in South et al (2020), where predictions are of the form $f_hat = K0*a + Phi*b$ for heldout K0 and Phi matrices and estimators using heldout samples are of the form $mean(f - f_hat) + b[1]$.

• ksd: (Only if diagnostics = TRUE) An estimated kernel Stein discrepancy based on the fitted model that can be used for diagnostic purposes. See South et al (2020) for further details.

• bound_const: (Only if diagnostics = TRUE and est_inds=NULL) This is such that the absolute error for the estimator should be less than $ksd \times bound_const$.

Warning

Solving the linear system in SECF has $O(N^3+Q^3)$ complexity where $N$ is the sample size and $Q$ is the number of terms in the polynomial. Standard SECF is therefore not suited to large $N$. The method aSECF is designed for larger $N$ and details can be found at aSECF and in South et al (2020). An alternative would be to use $est_inds$ which has $O(N_0^3 + Q^3)$ complexity in solving the linear system and $O((N-N_0)^2)$ complexity in handling the remaining samples, where $N_0$ is the length of $est_inds$. This can be much cheaper for small $N_0$ but the estimation of the Gaussian process model is only done using $N_0$ samples and the evaluation of the integral only uses $N-N_0$ samples.

On the choice of $\sigma$, the kernel and the Stein order

The kernel in Stein-based kernel methods is $L_x L_y k(x,y)$ where $L_x$ is a first or second order Stein operator in $x$ and $k(x,y)$ is some generic kernel to be specified.

The Stein operators for distribution $p(x)$ are defined as:

• steinOrder=1: $L_x g(x) = \nabla_x^T g(x) + \nabla_x \log p(x)^T g(x)$ (see e.g. Oates el al (2017))

• steinOrder=2: $L_x g(x) = \Delta_x g(x) + \nabla_x log p(x)^T \nabla_x g(x)$ (see e.g. South el al (2020))

Here $\nabla_x$ is the first order derivative wrt $x$ and $\Delta_x = \nabla_x^T \nabla_x$ is the Laplacian operator.

The generic kernels which are implemented in this package are listed below. Note that the input parameter sigma defines the kernel parameters $\sigma$.

• "gaussian": A Gaussian kernel,

  $k(x,y) = exp(-z(x,y)/\sigma^2)$

• "matern": A Matern kernel with $\sigma = (\lambda,\nu)$,

  $k(x,y) = bc^{\nu}z(x,y)^{\nu/2}K_{\nu}(c z(x,y)^{0.5})$

  where $b=2^{1-\nu}(\Gamma(\nu))^{-1}$, $c=(2\nu)^{0.5}\lambda^{-1}$ and $K_{\nu}(x)$ is the modified Bessel function of the second kind. Note that $\lambda$ is the length-scale parameter and $\nu$ is the smoothness parameter (which defaults to 2.5 for $steinOrder=1$ and 4.5 for $steinOrder=2$).

• "RQ": A rational quadratic kernel,

  $k(x,y) = (1+\sigma^{-2}z(x,y))^{-1}$

• "product": The product kernel that appears in Oates et al (2017) with $\sigma = (a,b)$

  $k(x,y) = (1+a z(x) + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

• "prodsim": A slightly different product kernel with $\sigma = (a,b)$ (see e.g. https://www.imperial.ac.uk/inference-group/projects/monte-carlo-methods/control-functionals/),

  $k(x,y) = (1+a z(x))^{-1}(1 + a z(y))^{-1} exp(-0.5 b^{-2} z(x,y))$

In the above equations, $z(x) = \sum_j x[j]^2$ and $z(x,y) = \sum_j (x[j] - y[j])^2$. For the last two kernels, the code only has implementations for steinOrder=1. Each combination of steinOrder and kernel_function above is currently hard-coded but it may be possible to extend this to other kernels in future versions using autodiff. The calculations for the first three kernels above are detailed in South et al (2020).

Author(s)

Leah F. South

References

Oates, C. J., Girolami, M. & Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3), 695-718.

South, L. F., Karvonen, T., Nemeth, C., Girolami, M. and Oates, C. J. (2020). Semi-Exact Control Functionals From Sard's Method. https://arxiv.org/abs/2002.00033

See Also

See ZVCV for examples and related functions. See SECF_crossval for a function to choose between different kernels for this estimator.

---