This section investigates the application of the previous results to Monte Carlo estimates constructed with the help of control variates in order to reduce the variance.
We start by presenting several applications in which uniform bounds for Monte Carlo methods are of interest. Then we give the mathematical background and finally provide sharp uniform error bounds for Monte Carlo estimates that use control variates.

\subsection{Uniformity in Monte Carlo procedures}

Proving uniform bounds on the error of Monte Carlo methods is motivated by the following three applications.

\paragraph{Latent variable models.}

Suppose we are interested in estimating the distribution of the variable $X$ whose density is assumed to lie in the model
\[
p_\theta (y) = \int p_\theta (y|z) \, p(z) \, \diff z
\]
where $\theta\in \Theta$ is the parameter to estimate, $p(z)$ is the known density of the so-called latent variable and $ \{(y,z)\mapsto p_\theta (y|z)\}_{\theta\in \Theta}$ is a model of conditional densities ($Y$ conditionally on the latent variable). This is actually a frequent situation in economics \citep{mcfadden2001economic} and medicine \citep[example 4, 6 and 9]{mcculloch2005generalized}. Given independent and identically distributed random variables $Y_1,\ldots,Y_n$ observed from the previous model with parameter $\theta_0$, the log-likelihood function takes the form
\[
\theta \mapsto \sum_{i=1}^n\log\left(  \int p_\theta(Y_i|z) \, p(z) \, \diff z\right).
\]
In most cases, each term in the previous sum is intractable and the approach proposed in \citet{mcfadden1994estimation} consists in replacing the unknown integrals by Monte Carlo estimates. In such a procedure, there is an additional estimation error compared to the statistical error of the standard maximum likelihood estimator. This additional error can be controlled in terms of the approximation error of $\int p_\theta(y|z) \, p(z) \, \diff z$ by the Monte Carlo estimate uniformly in $(y, \theta)$.

\paragraph{Stochastic programming.}

Consider the stochastic optimization problem
\[
	\min_{\theta\in \Theta} F(\theta) \qquad \text{with} \qquad F(\theta) = \expec [ f(\theta, X) ],
\]
where $X$ is a random variable in some space $\Xc$ with distribution $P$ and where $\Theta$ is a Euclidean set. The response functions are thus the maps $x \mapsto f(\theta, x)$ as $\theta$ ranges over $\Theta$. This problem is different from standard optimization because it takes into account some uncertainty in the output of the function $f$. One might think of the following toy example: $f$ is the output of a laboratory experiment, e.g., the amount of salt in a solution,  $\theta$ is the input of the experiment, e.g., the temperature of the solution and $X$ gathers unobserved random factors that influence the output. Optimizing such kind of functions is of interest in many different fields and we refer the reader to \citet{shapiro2014lectures} for more concrete examples (see also the example below that deals with quantile estimation). 
This problem might be solved using two competitive approaches: the sample average approximation (SAA) and stochastic approximation techniques such as gradient descent---see \citet{nemirovski2009robust} for a comparison between both approaches.
The SAA approach follows from approximating the function $F$ by a functional Monte Carlo estimate defined as
\[
   F_n (\theta ) = \frac{1}{n} \sum_{i=1}^n f(\theta, X_i) ,
\]
where $X_1,\ldots,X_n$ is a random sample from $P$. Then the minimizer of the first stochastic optimization problem is approximated by the minimizer $\theta_n$ of the functional estimate $F_n$. Following the reference textbook \citet{shapiro2014lectures}, the analysis of $\theta_n$ is carried out through error bounds for the Monte Carlo estimate $F_n(\theta)$ uniformly in $\theta \in \Theta$.
When sampling from $P$ is expensive, one may want to reduce the variance of the Monte Carlo estimate $F_n(\theta)$ by the method of control variates or some other method. In that case, a control on the error uniformly in $\theta$ is required.

\paragraph{Quantile estimation in simulation modeling.}

Quantiles are of prime importance when it comes to measure the uncertainty of random models \citep{law2000}. When the stochastic experiments are costly, variance reduction techniques such as the use of control variates are helpful \citep{hesterberg1998control,cannamela2008controlled}.
Let $F(y) = \pr \{ g(X) \leq y \}$, for $y\in \reals$, be the cumulative distribution function of a transformation $g : \Xc \to \reals$ of a random element $X$ in some space $\Xc$.
Suppose the interest is in the quantile $F^{-}(u) = \inf\{y \in \reals\,:\, F(y)\geq u\}$ for $u\in (0,1)$.
The functions $f_y : \Xc \to \reals$ to be integrated with respect to the distribution $P$ of $X$ are thus the indicators $x \mapsto f_y(x) = I\{g(x) \le y\}$, indexed by $y \in \reals$.
It is necessary to control the accuracy of an estimate of the probability $F(y) = \expec[f_y(X)]$ uniformly in $y \in \reals$ in order to have a control on the accuracy of an estimate of the quantile $F^{-}(u)$, even for a single $u \in (0, 1)$; see for instance Lemma~12 in \cite{portier2018weak}.
If drawing samples from $P$ or evaluating $g$ is expensive, it may be of interest to limit the number of Monte Carlo draws $X_i$ and function evaluations $g(X_i)$.

Finally, note that due to the formulation of a quantile as the minimiser of the expectation of the check function, see e.g., \citet{hjort2011asymptotics}, this example is an instance of the stochastic programming framework described before.

\paragraph{Related results.}
The previous examples underline the need of error bounds for Monte Carlo methods that are uniform over a family of response functions. The uniform consistency of standard Monte Carlo estimates over certain collections of functions can easily be shown by relying on Glivenko-Cantelli classes \citep{vdvaart+w:1996}; see for instance \citet{shapiro2014lectures} for applications to stochastic programming problems. Similarly, uniform convergence rates for standard Monte Carlo estimates can be derived from classical empirical process theory. We refer to \cite{gine+g:02} and the references therein for suprema over VC-type classes and to \cite{kloeckner2020empirical} and the references therein for suprema over H\"older-type classes. 
	
For variance reduction methods based on adaptive importance sampling, uniform consistency has been proven recently in \citet{delyon+p:2018} and \citet{feng2018uniform}. For control variates, however, we are not aware of any uniform error bounds and we believe the next results to be the first of their kind.

\subsection{Mathematical background for control variates}
\label{sec:CVback}

Let $\Fc \subset L^2(P)$ be a collection of square-integrable, real-valued functions $f$ on a probability space $(\mathcal{X}, \mathcal{A}, P)$ of which we would like to calculate the integral $P(f) = \int_{\mathcal{X}} f(x) \, \diff P(x)$. Let $X_1,\ldots,X_n $ be  independent random variables  taking values in $\mathcal{X}$ and with common distribution $P$.
The standard Monte Carlo estimate of $P(f)$ simply takes the form $P_n(f)=\frac{1}{n}\sum_{i=1}^n f(X_i)$. However, this estimator may converge slowly to $P(f)$ due to a high variance. To tackle this issue, it is common practice to use control variates, which are functions in $L^2(P)$ with known integrals. Without loss of generality, we can center the control variates $g_{n,1}\ldots,g_{n,d_n}$ and assume they have zero expectation, that is, $P(g_{n,k})=0$ for all $k\in \{1,\ldots,d_n\}$.
Let $g_n = (g_{n,1},\ldots,g_{n,d_n}) $ denote the $\mathbb R^{d_n}$-valued function with the $d_n$ control variates as elements and put $h_n = (1,g_n^\T)^\T$. Similarly as before, we assume that the Gram matrix $P(g_n g_n^\T)$ is invertible.  The control variate Monte Carlo estimate of $P(f)$ is given by $\hat \alpha_{n,f}$ defined as \citep[see for instance][Section~1]{portier+s:2019}, 
\begin{equation}
\label{eq:alphanf}
(\hat \alpha_{n,f}, \hat \beta_{n,f}) \in \argmin _ {\alpha \in \mathbb R, \, \beta \in \mathbb R^{d_n}}  P_n ( f   - \alpha -  g_n  ^\T  \beta)^2.
\end{equation}
The vector $\hat \beta_{n,f}$ contains the regression coefficients for the prediction of $f$ based on the covariates $h_n$. Remark that the control variate integral estimate $\hat \alpha_{n,f}$ coincides with the integral of the least square estimate of $f$, i.e., $\hat \alpha_{n,f} = P(\hat \alpha_{n,f}+ g_n  ^\T \hat \beta_{n,f})$. In addition, since $\hat \alpha_{n,f}$ can be expressed as a weighted estimate $\sum_{i=1}^n w_{i} f(X_i)$ where the weights $(w_i)_{i=1,\ldots, n}$ do not depend on the integrand $f$, there is a computational benefit to integrating multiple functions \citep[Remark 4]{leluc2019control}. It is useful to define
\begin{equation*}
(\alpha_{n,f},  \beta_{n,f}) \in \argmin _ {\alpha \in \mathbb R, \, \beta \in \mathbb R^d}  P  ( f   - \alpha -  g_n   ^\T  \beta)^2,
\end{equation*}
as well as the residual function 
\[
	\eps_{n,f}  = f - \alpha_{n,f} - g_n^\T  \beta_{n,f}.
\]
Note that $ \alpha_{n,f}  = P(f)$.
If $ \beta_{n,f} $ would be known, the resulting oracle estimator would be
\begin{equation}
\label{eq:alphanforacle}
\hat \alpha_{n,f}^{\mathrm{or}} = P_n[f-   g_n  ^\T \beta_{n,f} ].
\end{equation}
The question raised in the next section is whether the control variate estimate $\hat \alpha_{n,f}$ can achieve a similar accuracy uniformly in $f \in \Fc$ as the oracle estimator $\hat \alpha_{n,f}^{\mathrm{or}}$.


\subsection{Uniform error bounds}

Motivated by the examples above, we provide an error bound for the control variate Monte Carlo estimate 
$ \hat \alpha_{n,f}$ in~\eqref{eq:alphanf} uniformly in $f \in \Fc$. Before doing so, we give a uniform error bound for the oracle estimate $\hat \alpha_{n,f}^{\mathrm{or}}$ in~\eqref{eq:alphanforacle}. This serves two purposes: first, it will be useful in the analysis of $ \hat \alpha_{n,f}$ and second, it will provide sufficient conditions for the two estimates to achieve the same level of performance.
Recall $M_n$, $\gamma_n^2$ and $L_n^2$ in \eqref{eq:Mn} and \eqref{eq:Ln} and put
\[
	\sigma_n^2 = \supf P(\eps_{n,f}^2).
\]
A new assumption, $M_n^2 = \Oh \left(\sigma_n ^2  n / \log (n ) \right)$ as $n \to \infty$, in the same vein but weaker\footnote{Since $L_n^2 = M_n^2 \ninf{q_n}$ and $\gamma_n^2 \le \sigma_n^2 \ninf{q_n}$.} than $ L_n^2 = \Oh \left(\gamma_n ^2 \sqrt{ n / \log (n ) }\right)$, turns out to be useful to obtain the result.

\begin{theorem}[Uniform bound on error of oracle estimator]
	\label{prop:oracle_bound}
	Assume the framework of Section~\ref{sec:CVback} and suppose that Conditions \ref{cond:qn}, \ref{cond:F} and \ref{cond:VC} hold. If $\liminf _ {n\to \infty}  n^{\alpha} M_n> 0 $ for some $\alpha >0$ and if $ M_n^2 = \Oh \left(\sigma_n ^2  n / \log (n ) \right) $, then
\begin{align*}
\sup_{f\in \mathcal F}  \left|\hat \alpha_{n,f}^{\mathrm{or}} - P(f) \right| 
 &= \Op\left(\sigma_n \sqrt{  n^{-1} \log( n )} \right), \qquad n \to \infty.
\end{align*}
\end{theorem}

The proof of Theorem \ref{prop:oracle_bound} is provided in Appendix \ref{app:proofMC1}. The derivation of the stated rate relies on the property that the residual class $\Ec_n = \{ \eps_{n,f} : f \in \Fc\}$ is a VC-class of functions (as detailed in the proof of Proposition~\ref{prop:cover}). Indeed, noticing that 
$\hat \alpha_{n,f}^{\mathrm{or}} - P(f) =  P_n  (\eps_{n,f})$ allows to rely on the next proposition, dedicated to suprema of empirical process.

\begin{proposition}[Bound of supremum of empirical process]
	\label{prop:general term_Pn f}
	On the probability space $(\Xc, \Ac, P)$, let $\Sc$ be a VC-class of parameters $(w, B)$ with respect to the constant envelope $U \ge \sup_{\scfunc \in \Sc}\ninf{\scfunc}$. Suppose the following two conditions hold:
	\begin{compactenum}[(i)]
		\item
		$\tau^2 \ge \sup_{\scfunc\in\Sc} \Var_P(\scfunc)$ and $\tau \le 2U$;
		\item
		$w \ge 1$ and $B \ge 1 $. 
	\end{compactenum}
	Then, for $P_n$ the empirical distribution of an independent random sample $X_1, \ldots, X_n$ from $P$, we have with probability $1-\delta$: 
	\begin{align*}
		\sup_{\scfunc \in\Sc} \left| P_n(\scfunc) - P(\scfunc) \right|
		&\leq L \left(\tau \sqrt{w n^{-1} \log(L \theta/\delta)} + U w n^{-1} \log(L \theta/\delta ) \right),
	\end{align*}	
	with $\theta = B U / \tau $ and $L>0$ a universal constant.
\end{proposition}

The proof of Proposition~\ref{prop:general term_Pn f} is given in Appendix~\ref{supp:prop:general term_Pn f}. In the proof, we bound the expectation of the supremum by combining a well-known symmetrization inequality \citep[Lemma~2.3.1]{vdvaart+w:1996} with Proposition~2.1 in \citet{gine2001consistency}, and we will find a rate on the deviation of the supremum around its expectation by Theorem~1.4 in \citet{talagrand1996new}. Compared to existing results such as Proposition 2.2 stated in \cite{gine2001consistency}, our version is more precise due to the explicit role played by the VC constants in the bound.

The next result follows from an application of Corollary \ref{cor:simple-rate} and Theorem \ref{prop:oracle_bound} combined with some other bounds that are standard when analyzing control variates estimates.

\begin{theorem}[Uniform error bound on control variate Monte Carlo estimator]
	\label{th:unif_cv}
Assume the framework of Section~\ref{sec:CVback} and suppose that Conditions \ref{cond:qn}, \ref{cond:F} and \ref{cond:VC} hold. If $\liminf _ {n\to \infty}  n^{\alpha} M_n> 0 $ for some $\alpha >0$ and if $ L_n^2 = \Oh \left( \gamma_n ^2 \sqrt{ n / \log (n ) } \right) $, then
\begin{equation*}
\supf  \left|\hat \alpha _ f - P(f) \right| 
= \Op\left(   \sigma_n \sqrt{  n^{-1} \log( n )} \left( 1  + \sqrt {  d_n n^{-1} \ninf {q_n}       }   \right)  \right),
\qquad n \to \infty.
\end{equation*}
\end{theorem}

The proof of Theorem \ref{th:unif_cv} is provided in Appendix \ref{app:proofMC2}.  Compared to the error bound given in Theorem~\ref{prop:oracle_bound} for the oracle estimator, the error bound in Theorem~\ref{th:unif_cv} for the control variate estimator has an additional term. This term, which is due to the additional learning step that is needed to estimate the optimal control variate, vanishes as soon as $d_n \ninf{q_n} = \oh(n)$ as $n \to \infty$. This condition, which was used in \cite{newey1997convergence} as well as in \cite{portier+s:2019}, is meaningful as it relates the model complexity to the sample size, i.e., the computing time of the experiment.
{"config": "arxiv", "file": "2006.09223/MC_unif.tex"}