\section{Introduction}
It is difficult to overstate the importance and utility of resampling methods and the bootstrap in particular for determining properties of estimators whenever exact, explicit sampling distributions cannot be readily determined.  Given a sample $X_1, ..., X_n \sim \mathrm{P}$, a parameter of interest $\theta$, and an estimator $\hat{\theta} = s(X_1, ..., X_n)$, it is often of interest to estimate, for example, $\mathrm{Var}(\hat{\theta})$. Let $\bm{x}=(x_1, ..., x_n)$ denote the observed values of the sample so that for a particular realization, $\hat{\theta} = s(\bm{x})$.  To provide bootstrap estimate of the variability of the estimator, we can draw $B$ (re)samples of size $n$ with replacement from $\{x_1, ..., x_n\}$ to form bootstrap samples $\bm{x}^{*}_{1}, ..., \bm{x}^{*}_{B}$ from which we calculate bootstrap estimates $\hat{\theta}_1, ..., \hat{\theta}_B$.  The nonparametric bootstrap variance estimate of $\hat{\theta}$ is then taken as the empirical variance of $\hat{\theta}_1, ..., \hat{\theta}_B$ \cite{Efron1979,Efron2014}.

Within this context, given the necessity of calculating $\hat{\theta}_1, ..., \hat{\theta}_B$, it is natural to instead consider the estimator
\begin{equation}
	\label{eqn:bagest}
	\tilde{\theta}_B = \frac{1}{B} \sum_{b=1}^{B} \hat{\theta}_b
\end{equation}
\noindent as a ``bootstrap smoothed'' alternative to the original $\hat{\theta}$ \cite{EfronBootstrap}.  This sort of bootstrap aggregation (bagging) was also proposed by \cite{bagging} as a means by which predictive variance may be reduced when each bootstrap sample is used to construct an individual model, frequently a classification or regression tree.

The standard bootstrap approach -- referred to recently as the \emph{brute force} approach by \cite{Efron2014} -- to assessing the variability of $\tilde{\theta}_B$, though straightforward, is computationally quite burdensome, requiring several bootstrap replicates of not only the original data, but also from within the bootstrap samples themselves.  This double bootstrap, first proposed in \cite{Beran1988}, is especially costly whenever the original statistic $s$ is computationally costly.

A variety of approaches have been suggested to reduce the computational burden of these sorts of problems.  Some of these \citep{White2000,Davidson2000,Davidson2002,Davidson2007} employ what is now referred to as the fast double bootstrap whereby only a single second-level bootstrap sample is collected.  This approach was utilized by \cite{Giacomini2013} in running Monte Carlo experiments and a more careful analysis is given in \cite{Chang2015}.   Sexton and Laake \cite{SextonLaake} propose some alternative nonparametric means by which $\mathrm{Var}(\tilde{\theta}_B)$ may be estimated, suggesting also that the number of second-level bootstrap replicates $B^{'}$ need only be a fraction of the original resample size $B$.  In lieu of full bootstrap samples, subsampling, or $m$-out-of-$n$ bootstrap sampling, was proposed by \cite{Politis1994} and \cite{Bickel1997}.  More recently, \cite{Sengupta2016} proposed a combination of these approaches, first subsampling and then employing a single second-level resample.  Similarly, \cite{Kleiner2014} proposed the bag of little bootstraps which involves splitting the original dataset into a number of subsamples and then taking bootstrap samples on each subset allowing the process to more easily scale by being capable of efficiently running in parallel.

Though the above approaches can substantially reduce the computational complexity in estimating the variance of estimators based on resampling procedures, each nonetheless involves further resampling in order to obtain such an estimate.  Recently, motivated by the problem of taking into account not only the sampling variability but also the variability in model selection, Efron \cite{Efron2014} alleviated this issue by developing an algebraically compact, closed-form estimator for the variance of a bagged estimate.  Instead of additional resampling, Efron's proposal required only additional bookkeeping to recall which samples in the original data appeared how many times in each bootstrap sample.  This development was particularly beneficial in estimating the variance in predictions generated via supervised learning ensembles that are relatively computationally expensive.  A number of recent works, for example, have successfully applied this type of estimator in the context of variance estimation for supervised learning ensembles like random forests \cite{WagerIJ,wager2018estimation,zhou2019asymptotic}.

Though its final form is algebraically simple, the derivation of Efron's variance estimator is fairly involved and may appear somewhat mysterious to many readers.  Its development comes from an application of the original theory for the infinitesimal jackknife involving functional derivatives.  Likely as a result, studies rigorously investigating the statistical properties of this estimator as well as the contexts in which such an estimator would be appropriate are lacking in the current literature.  Efron, for example, notes that the appropriateness of his nonparametric delta method (infinitesimal jackknife) approach follows from the fact that the bagged estimates represent a more smooth function of the data.  Thus, while clearly an extremely significant result in and of itself, these estimates would not appear to apply to more general resampling schemes wherein such smoothness assumptions may not be reasonable.

In this paper, we strive to take a step forward both in better understanding the intuition behind this important estimator as well as in understanding its statistical properties.  In addition to the infinitesimal jackknife derivation utilized by Efron, we consider two alternative approaches that are a bit more straightforward and easily motivated.  The first of these exploits the important fact that conditional on the observed data, the bagged estimate in \cref{eqn:bagest} depends only on the resampling weights.  We consider a linear approximation to this function of bootstrap weights (i.e.\ standard linear regression) and demonstrate that this approach exactly reproduces the infinitesimal jackknife results given in \cite{Efron2014} whenever all bootstrap samples are employed.  As an additional benefit, this setup motivates a more general procedure for estimating the variance of any resampled estimate, not just one based on the bootstrap.

In addition to the linear regression and infinitesimal jackknife approaches, we also consider a classical jackknife motivation and once again demonstrate its equivalence in the full bootstrap context.  Importantly, this alternative representation of the estimator allows us to explore its asymptotic properties and in particular, the bias.  While the variance estimators motivated by the jackknife, infinitesimal jackknife, and linear regression approaches are shown to be identical when all bootstrap samples are used, they differ in practical settings when only a randomly selected subsample are employed, suggesting different bias corrections that might be imposed.

Finally, we derive the form of the infinitesimal jackknife estimate of variance in the U-statistic regime where the resampling is instead done by subsampling without replacement.  We discover that the variance estimators commonly employed in practice for quantifying the predictive uncertainty in supervised learning ensembles are actually something of a ``pseudo" infinitesimal jackknife in that they differ from the ``correct" form when properly derived.
While this difference is minor when the subsample size is small, in practical modern settings where large subsamples (possibly growing with $n$) are typically employed, such deviations become increasingly significant. We then investigate the properties of the pseudo infinitesimal jackknife and provide a consistent estimate of the variance of (generalized )U-statistics.

The remainder this paper proceeds as follows. In Section 2, we provide the background and historical motivation for the infinitesimal jackknife (IJ) method. In Section 3, in addition to Efron's recent derivation in \cite{Efron2014}, we provide three alternative approaches for obtaining variance estimates for any bootstrap smoothed statistic including via both the classical jackknife as well as straightforward linear regression.  We discuss the bias of the estimators here and demonstrate their equivalence when all bootstrap samples are employed. Finally, in Section 4, we derive the IJ estimator for the variance of (generalized) U-statistics and discuss its consistency.  To our knowledge, these results provide the first formal guarantees for the validity of confidence intervals constructed on such estimators in settings where all subsamples are not employed, and therefore apply directly to the kinds of modern supervised learning ensembles like random forests commonly employed in practice.

\section{Background of the Infinitesimal Jackknife (IJ)}
\label{sec:background}

Let $\mathcal{D}_n$ denote a sample of observed values from real-valued random variables $X_1,\dots, X_n$ that are i.i.d. from a distribution $\mathrm{P}$. In practice, we are often interested in estimating statistical functionals -- functions of the underlying distribution $\mathrm{P}$,  often estimated via the empirical distribution $\mathbb{P}_n$.  Denote this statistic as $s(X_1,\dots, X_n) = f(\mathbb{P}_n)$ and assume that $s$ is permutation symmetric in these $n$ arguments.  These ``functions of functions" were first introduced by  \cite{volterra1887} and today are a familiar topic of advanced analysis. Any statistic that treats the samples equivalently can also be viewed as a function of $\mathbb{P}_n$, albeit without always necessarily having an explicit form of $f$.  We can further extend the domain of $f$ to any non-negative functions on $X_1,\dots, X_n$ by defining
\begin{equation}
	f(\mathbb{P}) = f(c\cdot \mathbb{P}), \quad \text{for any } c>0.
\end{equation}

A common task, especially in today's big data era, is to find an appropriate and feasible means of estimating the variance of $f(\mathbb{P}_n)$.  Historically, there have been three primary methods:  the infinitesimal jackknife \cite{Miller1974}, influence curves \cite{hampel1974influence,huber1972}, and the delta method \cite{efron1982jackknife}. Though each method was motivated differently,  \cite{efron1981Nonparametric}  pointed out that the three methods are identical. We thus refer to the common estimator as $\mathrm{IJ}$, which is defined as
\begin{equation}
	\label{eqn:IJdef}
	\mathrm{IJ} =\frac{1}{n^2} \sum_i\mathrm{D}_i^2,
\end{equation}
where
\begin{equation}
	\label{di}
	\mathrm{D}_i = \lim_{\epsilon \to 0} \frac{f((1-\epsilon)\mathbb{P}_n+\epsilon \delta_{X_i}) - f(\mathbb{P}_n)}{\epsilon}
\end{equation}
and $\delta_{x}$ is the Dirac delta function.

\vspace{3mm}

We now briefly review the original derivation of the IJ, following closely to the original constructions given by von Mises in 1947 \cite{mises1947} and Jaeckel in 1972 \cite{Jaeckel1972}.  Let $\mathcal{P}$ be the set of all linear combinations of $\mathrm{P}$ and an arbitrary finite number of the $\delta_x$ measures. Let $\mathcal{P}^+$ be the the set of positive measures in $\mathcal{P}$, not including the zero measures and assume $f$ is defined for the probability measures in $\mathcal{P}^+$.  As above, extend $f$ to all of $\mathcal{P}^+$ by letting $f(c\cdot \mathrm{P}) =f(\mathrm{P})$ for all $c>0$. Note that $\mathcal{P}^+$ is convex and includes $\mathbb{P}_n$.

We now formally define the derivative of $f$. We say $f$ is differentiable at $\mathrm{G}$ in $\mathcal{P}^+$ if there exists a function $f'(\mathrm{G},x)$, defined at all $x$ in $\mathbb{R}$, with the following property:

\begin{definition}[ \cite{Jaeckel1972}]
	Let $\mathrm{H}$ be any member of $\mathcal{P}$ such that $\mathrm{G}+t\mathrm{H}$ is in $\mathcal{P}^+$ for all $t$ in some interval $0\leq t\leq t_{\mathrm{H}}$, $t_{\mathrm{H}}>0$, so that $f(\mathrm{G}+t\mathrm{H})$ is defined for $t$ in this interval. Then for any such $\mathrm{H}$, $f'(\mathrm{G},x)$ satisfies
	\begin{equation}
		\begin{aligned}
			\frac{df(\mathrm{G}+t\mathrm{H})}{dt}|_{t=0} &  :=  \lim_{t\to 0}\frac{f(\mathrm{G}+t\mathrm{H})-f(\mathrm{G})}{t} \\
			& = \int f'(\mathrm{G},x)\,d\mathrm{H}(x).
		\end{aligned}
	\end{equation}
\end{definition}

\vspace{4mm}

If $\mathrm{H}=\mathrm{G}$, we see that   $\int f'(\mathrm{G},x)\,d\mathrm{G}(x) =0$ since $f(c\mathrm{G})=f(\mathrm{G})$. On the other hand, if $\mathrm{H}=\delta_{x}- \mathrm{G}$, we find
\begin{equation}
	\label{ic}
	\lim_{t\to 0}\frac{f((1-t)\mathrm{G}+t\delta_{x})-f(\mathrm{G})}{t} = \int f'(\mathrm{G},x)\, d(\delta_x-\mathrm{G})(x) = f'(\mathrm{G},x).
\end{equation}
Indeed, \cite{hampel1974influence} defined $f'(\mathrm{G}, x)$ by  \cref{ic} and has called it the ``influence curve", since it reflects the influence of $f$ by adding a small mass on $\mathrm{G}$ at $x$.  Additionally, the derivative of $f(\mathrm{G}+t\mathrm{H})$ at arbitrary $t_0$ with $0<t_0<t_\mathrm{H}$ is given by
\begin{equation}
	\label{seg}
	\frac{df(\mathrm{G} + t\mathrm{H})}{dt}|_{t= t_0} = \frac{df(\mathrm{G}+t_0\mathrm{H} + u\mathrm{H})}{du}|_{u=0} =\int f'(\mathrm{G}+t_0\mathrm{H},x)\,d\mathrm{H}(x).
\end{equation}

Now assume that $f$ is differentiable (in the sense defined above) at all $\mathrm{G}$ in some convex neighbor of $\mathrm{P}$ in $\mathcal{P}^+$ such that $\mathbb{P}_n$ lies in the neighborhood with probability approaching one.  We now describe the motivation for addressing when the $\mathrm{IJ}$ estimator might provide a sensible estimate of the variance of $f(\mathbb{P}_n)$.  Parameterizing the segment from $\mathrm{P}$ to $\mathbb{P}_n$ by $\mathrm{P}(t) = \mathrm{P}+ t(\mathbb{P}_n-\mathrm{P})$ for $0\leq t\leq 1$, we have that if $\mathbb{P}_n$ lies in the neighborhood of $\mathrm{P}$, then
\begin{equation}
	\begin{aligned}
		f(\mathbb{P}_n) - f(\mathrm{P}) & =  f(\mathrm{P}(1)) - f(\mathrm{P}(0)) \\
		& = \frac{d f(\mathrm{P}(c))}{dt}|_{t=c} \quad \text{(by the Mean Value Theorem)}\\
		& = \int f'(\mathrm{P}(c))\, d(\mathbb{P}_n - \mathrm{P}) \quad \text{(by \cref{seg})}
	\end{aligned}
\end{equation}
for some $c$ in $[0,1]$.
Now, for large $n$, $\mathbb{P}_n $ is near $\mathrm{P}$ and we would expect that $f'(\mathrm{P}(c))$ is close to $f'(\mathrm{P})$, and thus
\begin{equation}
	\label{app-1}
	\begin{aligned}
		f(\mathbb{P}_n) - f(\mathrm{P}) & = \int f'(\mathrm{P},x)\, d(\mathbb{P}_n-\mathrm{P})(x) + o_p\left(\frac{1}{n}\right) \\
		& = \frac{1}{n}\sum_i f'(\mathrm{P},X_i)  + o_p\left(\frac{1}{n}\right)
	\end{aligned}
\end{equation}
where the last equality is due to the fact that $\int f'(\mathrm{P}, x)\, d\mathrm{P}(x)=0$.
Since the first term on the right-hand side is a sum of i.i.d. random variables, $\sqrt{n}(f(\mathbb{P}_n)-f(\mathrm{P}))$ is  asymptotic normal with mean $0$ and variance $\mathrm{V} = \int [f'(\mathrm{P},x)]^2\,d\mathrm{P}(x)$. Since $\mathrm{P}$ is unknown and $f'(\mathrm{P},x)$ depends on both $f$ and $\mathrm{P}$, we generally will not know $f'(\mathrm{P}, x)$ in advance, and so we could estimate $\mathrm{V}$ by
\begin{equation}
	\label{app-2}
	\begin{aligned}
		\mathrm{V} = \int [f'(\mathrm{P}, x)]^2\, d\mathrm{P}(x) &\approx \int f^{'2}(\mathbb{P}_n,x)\,d\mathbb{P}_n(x) \\
		& = \frac{1}{n}\sum_i  [f'(\mathbb{P}_n, X_i)]^2.
	\end{aligned}
\end{equation}
Then $\mathrm{Var}(f(\mathbb{P}_n))$ can be estimated by $\frac{1}{n^2}\sum  [f'(\mathbb{P}_n, X_i)]^2$, which corresponds to the definition of the IJ estimator given in \cref{eqn:IJdef} with $\mathrm{D}_i = f(\mathbb{P}_n, X_i)$.

In summary, to obtain the final estimate of the variance of $f(\mathbb{P}_n)$, we need to introduce two steps of approximation. First, in \cref{app-1}, we approximate $f(\mathbb{P}_n)$ with a linear statistic at $\mathbb{P}_n$.  Then, in \cref{app-2}, we approximate $\mathrm{P}$ with $\mathbb{P}_n$ and $f'(\mathrm{P})$ with $f'(\mathbb{P}_n)$.  Thus, in evaluating the quality of the IJ estimator, we must determine (i) whether $f$ is close to a linear statistic and (ii) whether $f'(\mathbb{P}_n)$ is close to $f'(\mathrm{P})$.  The following sections provide in-depth investigations into these issues for popular types of statistics formed by resampling.

\section{Infinitesimal Jackknife for Bootstrap ($\mathrm{IJ}_\mathrm{B}$)}
\label{sec: bootstrap}

We now focus on a more particular setting -- bootstrapping -- where the infinitesimal jackknife has seen the most success as a method for estimating variance.  Suppose that $s(X_1,..., X_n)$ is statistic, not necessarily a function of $\mathbb{P}_n$, and we take all possible bootstrap samples $(X_1^*,..., X_n^*)$, plug each into $s$ to  obtain a corresponding bootstrap estimate $s^*$, and then take the average. We call the new statistic the bootstrap smoothed (bagged) alternative of $s$ and denote it as $\mathbb{E}_*[s^*]$, where $\mathbb{E}_*[\cdot]$ denotes the expectation taken over the bootstrap sampling procedure conditional on the data. Note that $\mathbb{E}_*[s^*]$ is now a function of $\mathbb{P}_n$.

The dependence of $f$ on $\mathbb{P}_n$ can be explicitly expressed as
\begin{equation}
	f(\mathbb{P}_n) = \int s\, d\mathbb{P}_n\times \cdots \times \mathbb{P}_n = \int s\,d(\mathbb{P}_n)^n
\end{equation}
and so we see that $f$ depends on $(\mathbb{P}_n)^n$ and the dependence roughly exponentially grows with $n$.  According to the Kolmogorov theorem, $\sqrt{n}\sup_x|\mathbb{P}_n(x)-\mathrm{P}(x)|\to \sup_{t}\mathrm{B}(\mathrm{P}(t))$ in distribution, where $\mathrm{B}(t)$ is the Brownian bridge. Therefore, the distance between $\mathbb{P}_n$ and $\mathrm{P}$ is on the order of $1/\sqrt{n}$.  Note however that the distance between $(\mathbb{P}_n)^n$ and $(\mathrm{P})^n$ is only $\mathcal{O}(1)$ and therefore the distance between $f'(\mathbb{P}_n)$ and $f'(\mathrm{P})$ could be large if $f'(\mathbb{P}_n)$ depends on $\mathbb{P}_n$ exponentially.

\subsection{Three Approaches for Variance Estimation}
We turn now to the question of how to derive an estimator for $\mathrm{Var}(\mathbb{E}_*[s^*])$.  In the following subsections, we lay out three different approaches, including the Efron's infinitesimal jackknife formulation and ultimately demonstrate that all three are equivalent when all bootstrap samples are utilized. \\

\noindent \textbf{Method 1:  The Infinitesimal Jackknife Approach:}
Since $\mathbb{E}_*[s^*]$ can be viewed as a function of $\mathbb{P}_n$, estimating $\mathrm{Var}(\mathbb{E}_*[s^*])= \mathrm{Var}(f(\mathbb{P}_n))$ is a standard problem for the infinitesimal jackknife method and we denote the estimator $\mathrm{IJ_B}$. From the definition of $\mathbb{E}_*[s^*]$, we have
\begin{equation*}
	\begin{aligned}
		f((1-\epsilon)\mathbb{P}_n+\epsilon\delta_{X_i}) & = n^{-n} \sum\frac{s(X_1^*,\dots, X_n^*)n!}{(w_{1}^*!)(w_{2}^*!)\dots (w_{n}^*!)}\left[(1-\epsilon)^{\sum_{k\neq i} w_{k}^*} (1+(n-1)\epsilon)^{w_{i}^*}\right] \\
		&  =  n^{-n} \sum \frac{s(X_1^*,\dots, X_n^*)n!}{(w_{1}^*!)(w_{2}^*!)\dots (w_{n}^*!)}\left[1+n\epsilon(w_{i}^*-1)\right] + o(\epsilon^2) \\
		& = f(\mathbb{P}_n) + \epsilon n\mathrm{Cov}_*(s^*,w_i^*)  +  o(\epsilon^2),
	\end{aligned}
\end{equation*}
where $w_i^* = \#\{j:X_j^*=X_i\}$. By \cref{di}, $\mathrm{D}_i = n \mathrm{Cov}_*(s^*, w_i^*))= n\mathrm{Cov}_i$ and thus
\begin{equation}
	\label{IJ-bootstrap}
	\mathrm{IJ_B} =\sum_i \mathrm{Cov}_i^2 = \sum_i \mathrm{Cov}^2_*(s^*, w_i^*).
\end{equation}

The estimator in \cref{IJ-bootstrap} was first derived in 2014 by Efron \cite{Efron2014} as a straightforward application of the classical infinitesimal jackknife method.  This work, however, did not discuss how well the estimator might be expected to perform in various settings.  In fact, even going back to the original classical infinitesimal jackknife results, there has been little to no general theory for addressing such questions.  One possible explanation for the lack of answers in this area is that the process itself is quite abstract; it involves taking functional derivatives, which obfuscates the connection to classical probability theory.  Thus, we now provide an alternative derivation of the same estimator that avoids the need for functional derivatives. In particular, we formulate an explicit expression of $ f'(\mathrm{P},x)$ and a different interpretation of $\mathrm{Cov}_*(s^*,w_i^*)$. \\

\noindent \textbf{Method 1a:  An Alternative Infinitesimal Jackknife Derivation:}  First, note that $ t= \mathbb{E}_*[s^*]$ is a symmetric function of $X_1,\dots, X_n$ and therefore there exists an H-decomposition \cite{hoeffding1948central} of $t$.
To set up this H-decomposition, we first need to introduce following notation for kernels $t^{1}, \dots, t^{n}$ of degrees $1,\dots, n$. These kernels are defined recursively as follows
\begin{equation}
	t^{1}(x_1) = t_1(x_1)
\end{equation}
and
\begin{equation}
	t^c(x_1,\dots, x_c) = t_c(x_1,x_2,...,x_c) - \sum_{j=1}^c \sum_{i_1,\dots, i_j \in \{1,\dots, c\}} t^j(x_{i_1},\dots, x_{i_j})
\end{equation}
where $t_c(x_1,\dots, x_c) = \mathbb{E}[t(x_1,\dots, x_c, X_{c+1},\dots, X_n)] -\mathbb{E}[t]$. Then we can write $t= \mathbb{E}[t] + \sum_j \sum_{i_1,\dots, i_j}t^j (X_{i_1}, \dots X_{i_j})$.
Now consider using the linear term $l_b = \sum_i t^1(X_i)$ as an approximation of $t - \mathbb{E}[t]$.  Then we would expect
\begin{equation}
	\label{ap-1}
	\begin{aligned}
		t  - \mathbb{E}[t] =   l_b + o_p\left(\frac{1}{n}\right)
	\end{aligned}
\end{equation}
and naturally, we could use $\mathrm{Var}(l_b)$ as an estimate of $\mathrm{Var}(\mathbb{E}_*[s^*])$. Note however that $\mathrm{Var}(l_b) = n \int (t^1)^2\,d\mathrm{P}$ and both $\mathrm{P}$ and $t_1$ are unknown.  First, $\mathbb{P}_n$ is a good candidate for estimating $\mathrm{P}$, and we have
\begin{equation}
	\int (t^1)^2\,d\mathrm{P} \approx \int (t^1)^2\,d\mathbb{P}_n= \frac{1}{n}\sum_{i=1}^n (t^1(X_i))^2.
\end{equation}
Second, as for $t^1(X_1) = \mathbb{E}[t|X_1] - \mathbb{E}[t]$, $\mathbb{E}[\cdot]$ is again unknown, but we could substitute it with  $\mathbb{E}_*[\cdot]$ to obtain
\begin{equation}
	\label{ap-2}
	\begin{aligned}
		t^1(X_1) &  = \mathbb{E}[t(X_1,\dots, X_n)| X_1] - \mathbb{E}[t] \\
		& \approx  \mathbb{E}_*[t(X_1,X_2^*,\dots, X_n^*)] - \mathbb{E}_*[t] \\
		& =   \mathbb{E}_*[s^*(X_1,X_2^*, \dots, X_n^*)] - \mathbb{E}_*[s^*] \\
		& = e_1 - s_0
	\end{aligned}
\end{equation}
where $e_1=\mathbb{E}_*[s^*(X_1,X_2^*,\dots, X_n^*)]$ and $s_0 = \sum e_i/n$. The second equality in \cref{ap-2} holds because $t = \mathbb{E}_*[s^*]$ is already a bagged estimator.
Putting all of the above approximations together, we have
\begin{equation}
	\widehat{\mathrm{Var}(l_b)} = \sum_i (e_i-s_0)^2.
\end{equation}

It's worth pausing here to emphasize how the classical ideas behind the infinitesimal jackknife discussed in the previous sections (\cref{app-1} and \cref{app-2}) coincide with the idea of a linear approximation via the H-Decomposition.  In \cref{app-1}, we approximate $t = f(\mathbb{P}) $ by $f(\mathrm{P}) + \int  f'(\mathrm{P},x)\,d\mathbb{P}_n$ whereas in \cref{ap-1}, we approximate $t$ by $\mathbb{E}[t] + \sum_i t^1(X_i)$.  In \cref{app-2}, we approximate $f'(\mathrm{P},X_i)$ with $f'(\mathbb{P}_n, X_i)$; in \cref{ap-2} we approximate $t^1(X_i) = \mathbb{E}[t|X_1 = X_i] -\mathbb{E}[t]$ with $\mathbb{E}_*[t|X_1^* = X_i]- \mathbb{E}_*[t] = e_i-s_0$.  It follows that $f'(\mathrm{P}, x) = nt^1(x)$ and $f'(\mathbb{P}_n, X_i) = n(e_i-s_0)$.  Finally, this alternative derivation also suggests that $\mathrm{IJ_B}$ may only be appropriate for estimating $\mathrm{Var}(\mathbb{E}_*[s^*])$ in limited cases when the approximations in \cref{ap-1} and \cref{ap-2} are reasonable. \\

\noindent \textbf{Method 2:  The Jackknife Approach:}  We now propose an estimator for $\mathrm{Var}(\mathbb{E}_*[s^*])$ based on the classical jackknife procedure. Recall that delete-1 jackknife samples are selected by taking the original data vector and deleting one observation from the set. Thus, there are $n$ unique jackknife samples, the $i^{th}$ of which is given by $\mathcal{D}_n[i] = (X_1,\dots, X_{i-1},X_{i+1},\dots, X_n)$. The $i^{th}$ jackknife replicate is then defined as the value of the estimator evaluated on the $i^{th}$ jackknife sample. Letting $t$ denote $\mathbb{E}_*[s^*]$, the jackknife estimate of variance is then defined as
\begin{equation}
	\begin{aligned}
		\mathrm{JK} & = \frac{n-1}{n} \sum (t(\mathcal{D}_n[i])-\bar{t})^2\\
		& \approx n\mathrm{Var}_*\left(\sum_i t(\mathcal{D}_n[i])\cdot \mathbf{1}_{\{X_i\text{ is deleted}\}}\right).
	\end{aligned}
\end{equation}
Here, the key idea is that we expect $t(\mathcal{D}_n[i])$ be closed to $t(\mathcal{D}_n)$ and thus use the sample variance $t(\mathcal{D}_n[i])$ to estimate the variance of $t(\mathcal{D}_n)$.  Since those $t(\mathcal{D}_n[1]),\dots, t(\mathcal{D}_n[n])$ are strongly correlated, we scale the sample variance by $n$. Note that $t=\int s(x_1,\dots, x_n)\,d\mathbb{P}_n(x_1)\times \cdots \times \mathbb{P}_n(x_n)$.  Now suppose that we consider fixing the $i^{th}$ position rather than deleting the $i^{th}$ sample.  We would thus replace $t(\mathcal{D}_n[i])$ by
\begin{equation*}
	\begin{aligned}
		t_{(i,j)} & =\int s(x_1, x_2,\dots,x_{i-1}, X_j, x_{i+1},\dots, x_n)\,d\mathbb{P}_n(x_1)\cdots\mathbb{P}_n(x_{i-1})\times\mathbb{P}_n(x_{i+1})\cdots \mathbb{P}_n(x_n) \\
	\end{aligned}
\end{equation*}
for $ 1\leq i,j\leq n$ and propose
\begin{equation}
	\label{var}
	\begin{aligned}
		\mathrm{JK_B}
		& = n \mathrm{Var}_*\left(\sum_j t_{(i,j)}\mathbf{1}_{\{X_i^*=X_j\}}\right)  = \sum_j (e_j-s_0)^2\\
	\end{aligned}
\end{equation}
as an estimate of the variance of $\mathbb{E}_*[s^*]$. The third equality in \cref{var} is simply due to the fact that  $t_{(i,j)} = e_j$.\\

\noindent \textbf{Method 3:  The OLS Linear Regression Approach:} Recall that in a standard bootstrap procedure, $B$ resamples of size $n$ are sampled uniformly at random (with replacement) from the rows of $\mathcal{D}_n$.  Each observation in the original dataset receives equal weight and so we can write that weight vector as $\mathbf{w}^* \sim \text{Multinomial}(\frac{1}{n}, \dots, \frac{1}{n})$.  Now realize that conditional on the original data $\mathcal{D}_n$, each bootstrap estimate $s^*$ is a function of only those (empirical) weights $(w_1^*,\dots, w_n^*)$ corresponding to the observations actually selected in the bootstrap sample.  Consider the linear space spanned by $\mathbf{w}^*=(w_1^*,\dots, w_n^*)$ and denote the $l^*$ as the projection of $s(\mathbf{w}^*)$ onto that linear space.  We can use $\mathrm{Var}_*(l^*)$ as an estimate of $\mathrm{Var}(\mathbb{E}_*[s^*])$ which corresponds exactly to the setup where the relationship between the bootstrap estimates and observation weights is estimated via ordinary least squares linear regression.

\begin{remark}
	In the sections that follow, we rely more heavily on Methods 1 and 2 for assessing the properties of the infinitesimal jackknife variance estimator.  Nonetheless, the straightforward OLS approach in Method 3 reveals helpful and important insight into how the infinitesimal jackknife can be viewed.  A much more explicit walkthrough of how the familiar linear regression setup can be used to derive the infinitesimal jackknife estimator is provided in \cref{app:OLSconnection}.
\end{remark}

\subsection{Comparing the Three Approaches}
We now examine how the three estimators derived above compare to each other.  The first result below gives that the three estimators are identical whenever all bootstrap samples are used.
\begin{theorem}
	\label{3rivers}
	Suppose that we have data $\mathcal{D}_n$ and a statistic $s$. Let $(X_1^*, \dots, X_n^*)$ be a general bootstrap sample of $\mathcal{D}_n$, $s^*=s(X_1^*,\dots, X_n^*)$ and $w_j^*=\#\{i:X_i^*=X_j\}$, then
	\begin{enumerate}
		\item $\mathbb{E}_*[s^*w_j^*] = e_j$;
		\item $l^* =    \sum_j w_j^*\beta_j$, where $\beta_j = (e_j-s_0)$;
		\item 	$\mathrm{Var}_*(l^*) = \mathrm{JK_B}  = \mathrm{IJ_B}$.
	\end{enumerate}
	where $e_j=\mathbb{E}_*[s^*|X_1^*=X_j]$ and $s_0=\mathbb{E}_*[s^*]$.
\end{theorem}

The proof of \cref{3rivers} can be found in \cref{app:bootstrap}.

Now consider the more practical setting in which we draw only $B$ bootstrap samples $(X^*_{b1},\dots, X^*_{bn})$ and calculate $s_b^*= s(X_{b1}^*, \dots, X_{bn}^*)$ for each $b=1,\dots, B$.  Consider the bagged estimate $ \overline{s^*} = \frac{1}{B}\sum_{b=1}^B s_b^*$. By law of total variance, we have
\begin{equation}
	\label{dec}
	\begin{aligned}
		\mathrm{Var}(\overline{s^*}) & = \mathrm{Var}(\mathbb{E}[\overline{s^*}|\mathcal{D}_n])  +  \mathbb{E}[\mathrm{Var}(\overline{s^*}|\mathcal{D}_n)]  \\
		&= \mathrm{Var}(\mathbb{E}_*[s^*]) + \frac{1}{B}\mathbb{E}[\mathrm{Var}_*(s^*)].
	\end{aligned}
\end{equation}

\noindent The first term in \cref{dec} is dominant and so to provide a good estimate of $\mathrm{Var}(\overline{s^*})$, we must provide a good estimate of $\mathrm{Var}(\mathbb{E}_*[s^*])$. Since we do not employ all bootstrap samples, we cannot use $\mathrm{IJ}_\mathrm{B}$ directly but we can use the $B$ bootstrap samples to provide an estimate. First, a natural estimate of $\mathrm{Cov}_j$ is simply $\widehat{\mathrm{Cov}}_j$, the sample covariance of $(s_1^*, \dots, s_B^*)$ and $(w_{1j}^*, \dots, w_{Bj}^*)$.  Next, for $e_j-s_0$, we can first estimate $s_0$ with $\sum_{b=1}^B s_b^*/B$. Then, since $e_j$ the expected value of $s^*$ given $X_i^*=X_j$ for $i = 1,\dots, n$, a natural estimate is the weighted average of the mean of $s_b^*$ where $X_i^*=X_j$. Here, the weights are the proportion of times when $X_i^* = X_j$ across the $B$ bootstrap samples. A straightforward calculation gives that
\begin{equation}
	\widehat{e_j} = \sum_{b=1}^B\frac{w_{bj}^*}{\sum w_{bj}^*} s_b^*  \quad\text{and}\quad  \widehat{e_j}-\widehat{s_0} =  \sum_{b=1}^B
	\left(\frac{w^*_{bj} }{\sum w^*_{bj}} - \frac{1}{B}\right) s_b^*.
\end{equation}

Finally, the natural estimate of $\mathrm{Var}_*(l^*)$ is given by $\widehat{\mathrm{Var}}(\hat{l}) $, where $\hat{l}= (\hat{l}_1, \dots, \hat{l}_B)$ is the projection of $(s_1^*,\dots, s_B^*)$ onto the linear space spanned by $(w_{1j}^*, \dots, w_{Bj}^*)$ for $j=1,\dots, n$ and $\widehat{\mathrm{Var}}(\hat{l}) = \frac{1}{B}\sum_{b=1}^B ({\hat{l}_b - \bar{\hat{l}}})^2$ is the sample variance of $\hat{l}$.

Putting all of this together, we have the following three finite sample (i.e.\ not all bootstrap samples employed) variance estimators corresponding to the infinitesimal jackknife, jackknife, and OLS linear regression methods:

\begin{equation}
	\label{mc}
	\begin{aligned}
		\hat{\sigma}^{2}_{\text{IJ}} &:= \widehat{\mathrm{IJ}}_{\mathrm{B}} = \sum_{j=1}^n \widehat{\mathrm{Cov}}_j^2 = \sum_{j=1}^n\left[ \frac{1}{B-1} \sum_{b=1}^B (s_b^*-\bar{s^*}) (w^*_{bj}-\bar{w^*_{j}})\right]^2 \\
		\hat{\sigma}^{2}_{\text{JK}} &:= \widehat{\mathrm{JK}}_\mathrm{B} =\sum_{j=1}^n (\widehat{e_j}-\widehat{s_0})^2 = \sum_{j=1}^n \left[\sum_{b=1}^B
		\left(\frac{w^*_{bj} }{\sum w^*_{bj}} - \frac{1}{B}\right) s_b^*\right]^2\\
		\hat{\sigma}^{2}_{\text{OLS}} &:= \widehat{\mathrm{Var}}(\hat{l}) = \frac{1}{B}\sum_{b=1}^B ({\hat{l}_b - \bar{\hat{l}}})^2. \\
	\end{aligned}
\end{equation}

While these alternative derivations and estimators are helpful for providing insight into the infinitesimal jackknife in general, $\hat{\sigma}^{2}_{\text{IJ}} = \widehat{\mathrm{IJ}}_{\mathrm{B}}$ seems to be the estimator that has garnered the most practical interest.  In the following subsections, we focus on assessing its properties.

\subsection{The Bias of $\widehat{\mathrm{IJ}}_\mathrm{B}$}
\label{sec:biasIJ}

We now begin to analyze the properties of $\widehat{\mathrm{IJ}}_{\mathrm{B}}$ as an estimator.  As discussed above in reference to \cref{dec}, in order to understand how well $\widehat{\mathrm{IJ}}_{\mathrm{B}}$ estimates $\mathrm{Var}(\overline{s^*})$, we need only understand how well it estimates $\mathrm{Var}(\mathbb{E}_*[s^*])$.  Here we consider both the Monte Carlo bias and sampling bias of $\widehat{\mathrm{IJ}}_{\mathrm{B}}$ where the sampling bias is considered with respect to variation in the data, whereas the Monte Carlo bias is considered with respect to bootstrap process conditional on the data. We combine these two sources of bias with
\begin{equation}
	\widehat{\mathrm{IJ}}_\mathrm{B} - \mathrm{Var}(\mathbb{E}_*[s^*])   \propto  \underbrace{\mathbb{E}_*[\widehat{\mathrm{IJ}}_\mathrm{B}] -\mathrm{IJ}_\mathrm{B}}_{\textbf{Monte Carlo Bias}} + \underbrace{ \mathbb{E}[\mathrm{IJ}_\mathrm{B}] - \mathrm{Var}(\mathbb{E}_*[s^*])}_{\textbf{Sampling Bias}}.
\end{equation}

\noindent and consider their impact separately in the following subsections. \\

\subsubsection{Monte Carlo Bias}
We first consider the Monte Carlo bias of  $\widehat{\mathrm{IJ}}_\mathrm{B}$.  Note that $\mathbb{E}_*[\widehat{\mathrm{IJ}}_\mathrm{B}] -\mathrm{IJ}_\mathrm{B} = \sum_j \mathbb{E}_*[\widehat{\mathrm{Cov}}_j^2] - \mathrm{Cov}_j^2$ and $\mathbb{E}_*[\widehat{\mathrm{Cov}}_j^2 ] - \mathrm{Cov}_j^2
= \mathbb{E}_*[\widehat{\mathrm{Cov}}_j^2 ] - \mathbb{E}_*^2[\widehat{\mathrm{Cov}}_j]   = \mathrm{Var}_*(\widehat{\mathrm{Cov}}_j)$. Next, we have
\begin{equation*}
	\begin{aligned}
		\mathrm{Var}_*(\widehat{\mathrm{Cov}}_j)
		&  =  -\frac{B-2}{B(B-1)}\mathrm{Cov}^2_*(s^*,w_j^*)  + \frac{\mathrm{Var}_*(s^*)\mathrm{Var}_*(w_j^*)}{B(B-1)} +\\ & ~~~~\frac{\mathbb{E}_*[(s^*-\mathbb{E}_*[s^*])^2(w_j^*-\mathbb{E}_*[w_j^*])]^2}{B} \\
		& := \mathrm{I} + \mathrm{II}
	\end{aligned}
\end{equation*}
where
\begin{equation}
	\begin{aligned}
		\mathrm{I} &= \frac{1}{B}\mathrm{Var}_*\left((s^*-\mathbb{E}_*[s^*])(w_j^*-\mathbb{E}_*[w_j^*])\right) \\
		\mathrm{II} & = \frac{1}{B(B-1)}\left[\mathrm{Var}_*(s^*)\mathrm{Var}_*(w_j^*) + \mathrm{Cov}_*^2(s^*,w_j^*)\right].
	\end{aligned}
\end{equation}
The first term $\mathrm{I}$ is the dominant term and is $\mathcal{O}(1/B)$. Essentially, we see here that using $\widehat{\mathrm{Cov}}_j^2$ to estimate $\mathrm{Cov}_j^2$ is analogous to using $\bar{X}^2$ to estimate $\mathbb{E}^2[X]$, which is biased since $\mathbb{E}[\bar{X}^2] - \mathbb{E}^2[X] = \mathrm{Var}(X)/B$.  Indeed, $\mathrm{Var}(X)/B$ may not be negligible, especially when $B$ is small and  the coefficient of variation of $X$ is large. The variance estimator defined below offers a bias correction for $\widehat{\mathrm{IJ}}_\mathrm{B}$.

\begin{definition}
	A Monte Carlo bias corrected version of $\widehat{\mathrm{IJ}}_\mathrm{B}$ is given by
	\begin{equation}
		\label{mc-ij}
		\widehat{\mathrm{IJ}}_\mathrm{B}^{mc} = \widehat{\mathrm{IJ}}_\mathrm{B} -  \frac{1}{B} \sum_j \widehat{\mathrm{Var}}\left((s^*-\bar{s^*})(w_j^*-\bar{w_j^*})\right),
	\end{equation}
	where $\widehat{\mathrm{Var}}$ denotes sample variance.
\end{definition}

Note that the bias correction term above is a sum over $n$ terms. Then if $B$ is small, the bias correction term will be significant.
\begin{remark}
	In recent work, \cite{WagerIJ} proposed the Monte Carlo bias corrected estimator
	\begin{equation}
		\label{mc-efron}
		\widehat{\mathrm{IJ}}_\mathrm{B}^{whe}= \widehat{\mathrm{IJ}}_\mathrm{B}-  \frac{n}{B}\widehat{\mathrm{Var}}(s^*).

	\end{equation}
	We see that if $\mathrm{Var}_*((s^*-\mathbb{E}_*[s^*])(w_j^*-\mathbb{E}_*[w_j^*])$ is close to $\mathrm{Var}_*(s^*-\mathbb{E}_*[s^*])\mathrm{Var}_*(w_j^*-\mathbb{E}^*[w_j^*]) = (1-\frac{1}{n})\mathrm{Var}_*(s^*-\mathbb{E}_*[s^*])$, then \cref{mc-ij} is close to \cref{mc-efron}.  Simulations comparing \cref{mc-ij} and  \cref{mc-efron} are provided in the next section. \\
\end{remark}

\subsubsection{Sampling Bias}
Before looking more generally at the sampling bias of the infinitesimal jackknife, let's first examine how $\mathrm{IJ}_\mathrm{B}$ behaves on some simple examples.

\vspace{3mm}

\noindent \textbf{Example 1: Sample Mean}
Consider $s=s(x_1,\dots, x_n) = \frac{1}{n}\sum_{i=1}^n x_i$. We have
\begin{equation}
	s^* = \frac{1}{n}\sum_{i=1}^n X_i^*, \quad l^* = \sum_{i=1}^n \sum_{j=1}^n 1_{X_i^*=X_j} (e_j-s_0).
\end{equation}
Then, $\mathbb{E}_*[s^*] = \frac{1}{n}\sum_{i=1}^n X_i$ and  $\mathrm{Var}_*(l^*) =\frac{1}{n^2}\sum_{i=1}^n (X_i-\bar{X})^2$. Therefore,
\begin{equation}
	\mathrm{Var}(\mathbb{E}_*[s^*]) = \sigma^2/n, \quad \mathbb{E}[\mathrm{Var}_*(l^*)] = (n-1)\sigma^2/ n^2
\end{equation}
and thus, we have  $ \frac{\mathbb{E}[\mathrm{Var}_*(l^*)]}{\mathrm{Var}(\mathbb{E}_*[s^*])} = \frac{n-1}{n}\to 1$ as $n\to \infty$. In  Figures 1 and 2, $X_1,\dots, X_n$ follow $\mathcal{N}(0, \sigma^2)$, $n=100$ and $\sigma^2 =1$. Since we know that $\mathrm{Var}(\mathbb{E}_*[s^*])=\sigma^2/n$ , an oracle estimate would be $\widehat{\sigma^2}/n$, where $\widehat{\sigma^2}$ is the sample variance. The gray dashed line denotes the true value of $\mathrm{Var}(\mathbb{E}_*[s^*])$.  We find that $\widehat{\mathrm{IJ}}_\mathrm{B}^{mc}$ and $\widehat{\mathrm{IJ}}_\mathrm{B}^{whe}$ are quite close as expected and both perform well. The original $\widehat{\mathrm{IJ}}_\mathrm{B}$ seems to overestimate substantially when $B=100$.

\begin{figure}[h]
	\centering
	\begin{minipage}{0.45\linewidth}
		\label{bpmean}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{mean-100.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample mean (B=100).}
	\end{minipage}
	\quad
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{mean-1000.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample mean (B=1000).}
	\end{minipage}
\end{figure}

\vspace{3mm}
\noindent \textbf{Example 2: Sample Variance}
Consider  $s= {n \choose 2}^{-1}\sum \limits_{i<j} (x_i-x_j)^2$. We have
\begin{equation*}
	\mathbb{E}_*[s^*] = \frac{1}{n}\sum_{i=1}^n(X_i-\bar{X})^2\quad \text{and}\quad \mathrm{Var}_*(l^*) = \frac{1}{n^2} \sum_{i}\left[(X_i-\bar{X})^2- \frac{1}{n}\sum_{i}(X_i-\bar{X})^2\right]^2.
\end{equation*}
Then we have
\begin{equation}
	\begin{aligned}
		\mathrm{Var}(\mathbb{E}_*[s^*]) & =\left(\frac{n-1}{n}\right)^2\left[\frac{\mu_4}{n} - \frac{\mu_2^2}{n}\frac{n-3}{n-1}\right] \\
		& = a_n\mu_4 - b_n\mu_2^2,
	\end{aligned}
\end{equation}
where $\mu_i$ is the $i$th central moment of $X$ for $i=2,4$. Let $\mathbf{X}= (X_1,\dots, X_n)^T$, then
$\mathbb{E}[\mathrm{Var}_*(l^*)] $ can be written as  $ \frac{1}{n} \mathbb{E}[\mathbf{X}^T\mathrm{A}\mathbf{X}]^2$,
where  $\mathrm{A} = \Sigma_1 - \frac{1}{n}\sum_i \Sigma_i$, $\Sigma_i = (e_i-\frac{1}{n}\mathbf{1}_n)(e_i^T -\frac{1}{n}\mathbf{1}_n^T)$ and $e_i=(0,\dots, 0,1,0,\dots, 0)$. After some calculation, we obtain
\begin{equation*}
	\begin{aligned}
		& ~~~~ \mathbb{E}[\mathrm{Var}_*(l^*)]  \\
		&= \frac{(n-1)}{n^2}\left[\mathbb{E}[(X_1-\bar{X})^4] -\mathbb{E}[(X_1-\bar{X})^2(X_2-\bar{X})^2]\right] \\
		& = \left(\frac{n-1}{n}\right)^2 \left[ \left(\frac{n^3-(n-1)^2}{n^2(n-1)^2} + \frac{n}{(n-1)^5}\right)\mu_4 - \left(\frac{n^2-2n+3}{(n-1)n^2} - \frac{3n^2(2n-3)}{(n-1)^5}\right)\mu_2^2\right]\\
		& = a_n'\mu_4 - b_n' \mu_2^2.
	\end{aligned}
\end{equation*}
Thus, we have
\begin{align}
	\nonumber
	\frac{a_n'}{a_n} &=   1 + \frac{n^2+n-1}{n(n-1)^2} + \frac{n^2}{(n-1)^5}  = 1 + \frac{1}{n} + o(\frac{1}{n}) \\
	\nonumber
	\frac{b_n'}{b_n} & = 1 + \frac{n+3}{n(n-3)}  -  \frac{3n^3(2n-3)}{(n-1)^4(n-3)} = 1- \frac{5}{n} + o(\frac{1}{n}).
\end{align}
Since $a_n'/a_n\to 1 $ and $b_n'/b_n\to 1$ as $n\to \infty$, we have $\frac{\mathbb{E}[\mathrm{Var}_*(l^*)]}{\mathrm{Var}(\mathbb{E}_*[s^*])} \to 1$.  $\mathrm{IJ}_\mathrm{B}$ is therefore asymptotically unbiased for estimating  the variance of the sample variance. Since the sample variance is close to a linear statistic, the result is not surprising. In Figures 3 and 4,  $X_1,\dots, X_n$ follow $\mathcal{N}(0, \sigma^2)$, $n=100$,  and $\sigma^2 =1$. Since we know $\mathrm{Var}(\mathbb{E}_*[s^*]) = 2\sigma^4/n$, an oracle estimate would be $2(\widehat{\sigma^2})^2/n$, where $\widehat{\sigma^2}$ is  the sample variance. The gray dashed line denotes the true value of $\mathrm{Var}(\mathbb{E}_*[s^*])$. As in the first example, $\widehat{\mathrm{IJ}}_\mathrm{B}^{mc}$ and $\widehat{\mathrm{IJ}}_\mathrm{B}^{whe}$ are both quite close and perform well. The original $\widehat{\mathrm{IJ}}_\mathrm{B}$ again seems to suffer from overestimation when $B=100$.

\begin{figure}[h]
	\centering
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{var-100.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample variance (B=100).}
	\end{minipage}
	\quad
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{var-1000.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample variance (B=1000).}
	\end{minipage}
\end{figure}

\vspace{3mm}
\noindent \textbf{Example 3: Sample Maximum}
Consider $s=\max \{X_1,\dots, X_n\}$, where $X_1,\dots, X_n$ are uniformly distributed in $[0,1]$. Then
\begin{equation}
	\mathrm{Var}(\mathbb{E}_*[s^*]) = \frac{(n-\sum_{j=1}^n (\frac{j-1}{n})^n) (1+ \sum_{j=1}^n (\frac{j-1}{n})^n)}{(n+1)^2(n+2) }
\end{equation}
and
\begin{equation}
	\mathbb{E}[\mathrm{Var}_*(l^*)] = \frac{\sum_i [ \sum_{j=1}^{n} (\frac{j-1}{n})^n - \sum_{j=i+1}^{n}(\frac{j-1}{n})^{n-1}]^2}{(n+1)(n+2)}.
\end{equation}
The details of the calculation can be found in \cref{app:bootstrap}. We have
\begin{equation}
	\begin{aligned}
		\frac{\mathbb{E}[\mathrm{Var}_*(l^*)]}{\mathrm{Var}(\mathbb{E}_*[s^*])}
		& =  \frac{(n+1)\sum_i [ \sum_{j=1}^{n} (\frac{j-1}{n})^n - \sum_{j=i+1}^{n}(\frac{j-1}{n})^{n-1}]^2}{(n-\sum_j  (\frac{j-1}{n})^n) (1+ \sum_j (\frac{j-1}{n})^n)}\\
		& \to c  \in [0.24,0.25]\quad  (n \to \infty)
	\end{aligned}
\end{equation}
and thus we see that $\mathrm{IJ}_\mathrm{B}$ is underestimating of $\mathrm{Var}(\mathbb{E}_*[s^*])$ by a considerable margin. In this case, $\mathbb{E}_*[s^*]$ is not close to a linear statistic, so $\mathrm{IJ}_\mathrm{B}$ should not be expected to perform well. In Figures 5 and 6, $X_1,\dots, X_n$ follow Unif(0, 1) and $n=100$ with the dashed line corresponding to the true value of $\mathrm{Var}(\mathbb{E}_*[s^*])$. In this case, there is no obvious oracle estimator for $\mathbb{E}_*[s^*]$.  Unlike the previous two examples, although $\widehat{\mathrm{IJ}}_\mathrm{B}^{mc}$ and $\widehat{\mathrm{IJ}}_\mathrm{B}^{whe}$ remain quite similar, all three estimators suffer from considerable underestimation even when $B=1000$.

\begin{figure}[h]
	\centering
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{max-100.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample maximum (B=100).}
	\end{minipage}
	\quad
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{max-1000.pdf}
		\caption{Performance of the infinitesimal jackknife and its bias-corrected alternatives on estimating the variance of the  bagged sample maximum (B=1000).}
	\end{minipage}
\end{figure}

The following theorem provides an equivalent condition for $\mathrm{IJ}_\mathrm{B}$ to be asymptotically unbiased and suggests that $\mathbb{E}_*[s^*]$ needs to be asymptotically linear.

\begin{theorem}
	\label{BOOT-2}
	Let $\mathbb{E}_*[s^*]$ be the bootstrap smoothed alternative of $s$, then
	\begin{equation}
		\lim_{n\to \infty} \frac{\mathbb{E}[\mathrm{Var}_*(l^*)]}{\mathrm{Var}(\mathbb{E}_*[s^*])} = 1 \iff  \lim_{n\to \infty} n(1-\rho) = 1,
	\end{equation}
	where $\rho$ is the correlation coefficient between $e_1$ and $e_2$ and $e_i = \mathbb{E}_*[s^*|X_1^*=X_i]$.
\end{theorem}

To help provide some intuition, consider the case where $s=\bar{X}$, which is linear. Here, $\rho= \frac{n^2-1}{n^2+n-1} = 1 -1/n + o(1/n)$ and $\mathbb{E}[\mathrm{Var}_*(l^*)] / \mathrm{Var}(\mathbb{E}_*[s^*]) = \frac{n-1}{n}\to 1$. We suspect that to make $\rho = 1-1/n + o(1/n)$, $\mathbb{E}_*[s^*]$ is required to be equal to $l_b + o_p(1/n)$. From the original derivative of IJ and  the above examples, we can see that it is fairly hard for $\mathrm{IJ}_\mathrm{B}$ to estimate $\mathrm{Var}_*(\mathbb{E}_*[s^*])$ well. In general, it is difficult to understand the bias and other statistical properties of $\mathrm{IJ}_\mathrm{B}$ comprehensively, largely due to the replicates of bootstrap samples.

\section{U-statistics and Infinitesimal Jackknife}
\subsection{Infinitesimal Jackknife for U-statistics ($\mathrm{IJ}_\mathrm{U}$)}

While the infinitesimal jackknife is naturally appealing for estimating the variance of bagged statistics, the current resurgence in interest is due in large part to its potential for quantifying the uncertainty of complex supervised machine learning ensembles often formed by subsampling.  Immediately following Efron's work in \cite{Efron2014}, \cite{WagerIJ} showed that the same infinitesimal jackknife approach could be used to generate confidence intervals for random forests.  In recent years, a number of works have expanded on this idea \cite{Mentch2016,Wager2018,peng2019asymptotic}.  Among the key breakthroughs was the realization that predictions generated by averaging across ensembles of regression estimates, each formed with subsamples of the original data, could be seen as akin to classical U-statistics.

We thus now explore how the infinitesimal jackknife can be extended to the case of subsampling without replacement. In this case,  $\mathbb{E}_*[s^*]$ is a U-statistic, which is often more convenient for theoretical analysis and also more likely to be close to linear.  Here $s$  is a (permutation-symmetric) function of $k$ i.i.d. random variables and the U-statistic can be written as
\begin{equation}
	\begin{aligned}
		\mathrm{U} & = {n\choose k}^{-1}\sum_{(n,k)} s(X_{i_1}, \dots, X_{i_k}) \\
	\end{aligned}
\end{equation}
where the sum is taken over all ${n\choose k}$ subsamples of size $k$.
How does U depend on $\mathbb{P}_n$, such that $\mathrm{U} = f(\mathbb{P}_n)$ for some $f$?  The dependence is abstract so that the subsampling proceeds according to the probabilities determined by $\mathbb{P}_n$.  Following directly from the original definition of the IJ, we arrive at the following theorem.
\begin{theorem}
	\label{U-ij-1}
	The IJ estimator of the variance of a U-statistic is given by
	\begin{equation}
		\label{def:iju}
		\mathrm{IJ}_\mathrm{U}= \frac{k^2}{n^2}\sum_{j=1}^n[\alpha e_j-\beta s_0]^2,
	\end{equation}
	where $e_j = \mathbb{E}_*[s^*|X_1^*=X_j]$, $s_0=\mathbb{E}_*[s^*]$ and
	$$\alpha = 1 + \frac{1}{n}\left\{\frac{k-1}{2} - \frac{1}{k}\sum_{j=0}^{k-1}\frac{j^2}{(n-j)}\right\}, \quad
	\beta = 1 + \frac{1}{k}\sum_{j=0}^{k-1}\frac{j}{n-j}.$$
\end{theorem}
To understand the bias of $\mathrm{IJ}_\mathrm{U}$, the H-decomposition will prove quite useful. To set up this H-decomposition, we first need to introduce following notation for kernels $s^{1}, \dots, s^{k}$ of degrees $1,\dots, k$. These kernels
are defined recursively as follows
\begin{equation}
	s^1(x_1) = s_1(x_1)
\end{equation}
and
\begin{equation}
	s^c(x_1,\dots, x_c) = s_c(x_1,x_2,...,x_c) - \sum_{j=1}^c \sum_{i_1,\dots, i_j \in \{1,\dots, c\}} s^j(x_{i_1},\dots, x_{i_j})
\end{equation}
where $s_c(x_1,\dots, x_c) = \mathbb{E}[s(x_1,\dots, x_c, X_{c+1},\dots, X_k)] -\mathbb{E}[s]$. Let $V_j = \mathrm{Var}(s^{j})$ for $j=1,\dots, k$.  Then $\mathbb{E}[\mathrm{IJ}_\mathrm{U}]$ can both be written as a linear combination of those $V_j$. In particular, we have the following theorem.
\begin{theorem}
	\label{U-ij-2}
	Let $\theta=\mathbb{E}[s]$ and $\mathrm{IJ}_\mathrm{U}$ be as defined in \cref{def:iju}.  Then
	\begin{equation}
		\mathbb{E}[\mathrm{IJ}_\mathrm{U}] = \sum_{j=1}^k r_j {k\choose j}^2{n\choose j}^{-1} V_j + \frac{k^2}{n}(\alpha-\beta)^2\theta^2,
	\end{equation}
	where
	\begin{equation}
		r_j=   \frac{(n-k)^2}{n^2}\left[\frac{j}{1-j/n} \alpha^2\right]+\frac{k^2}{n} (\alpha-\beta)^2 , \quad \text{for } j=1,\dots,k.
	\end{equation}

\end{theorem}

\begin{remark}
	Note that $\mathrm{Var}(\mathrm{U}) = \sum_{j=1}^k {k\choose j}^2{n\choose j}^{-1}V_j$. If $k$ is fixed, then $\alpha, \beta \to 1$ and thus  $r_j\to j$ for $j=1,\dots, k$. Since $\mathrm{Var}(\mathrm{U})$ is dominated by the $V_1$ term, $\mathrm{IJ}_\mathrm{U}$ is asymptotically unbiased.
\end{remark}

\subsection{The Pseudo Infinitesimal Jackknife for U-statistics ($\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$)}

In recent work, \cite{Wager2018} investigated the consistency and asymptotic normality of random forests where the individual trees were constructed with subsamples of the original data.  As part of this work, the authors proposed another estimate of variance wherein the format of the IJ for bootstrap was simply copied over to this subsampling regime to arrive at
\begin{equation}
	\label{pseudo-IJ-U}
	\sum_j \mathrm{Cov}^2_*(s^*, w_j^*),
\end{equation}
where $*$ refers to the subsampling procedure.  We refer to this estimator as $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$, since it is not derived from the definition of the infinitesimal jackknife.   It is, however, possible to provide a more rigorous motivation for $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$.
Recall from Section 2 (see \cref{app-1}), we assume that $f(\mathbb{P}_n)-f(\mathrm{P})$ can be written as $\frac{1}{n}\sum_i f'(\mathrm{P},X_i) + o_p(1/n)$, where the dominant term is a sum of i.i.d. random variables and we estimate the variance of $f'(\mathrm{P},X_i)$ by $\frac{1}{n}\sum_i f^{'2}(\mathbb{P}_n,X_i)$. Now suppose that we write $f(\mathbb{P}_n) - f(\mathrm{P})$ as  $\sum_i g(X_i) + o_p(1/n)$, where $g(X_i)$ is not necessarily $\frac{1}{n}f'(\mathrm{P}, X_i)$. From classical U-statistic theory, we know that there is a natural candidate for $g(X_i)$:  the H\'ajek projection - $\mathbb{E}[f(\mathbb{P}_n)-f(\mathrm{P})|X_i] = \frac{k}{n}\mathbb{E}[s-\mathbb{E}[s]|X_i]$. Since $\mathrm{Var}(\sum g(X_i)) = \frac{k^2}{n}V_1$, where $V_1 = \mathrm{Var}(\mathbb{E}[s|X_1])$, we need only to propose a reasonable estimate for $V_1$ and we can then use
\begin{equation}
	\label{l-u}
	\frac{k^2}{n}\hat{V}_1
\end{equation}
as an estimate of the variance of the U-statistic. Since $V_1 = \mathbb{E}[\mathbb{E}[s|X_1]-\mathbb{E}[s]]^2$, a natural candidate of $\hat{V}_1$ would be
\begin{equation}
	\frac{1}{n}\sum_j (\mathbb{E}_*[s^*|X_1^*=X_j] - \mathbb{E}_*[s^*])^2 = \frac{1}{n}\sum_j (e_j-s_0)^2.
\end{equation}
As it turns out that \cref{l-u} is the same as \cref{pseudo-IJ-U}.
\begin{proposition}
	\label{U-sij-pre}
	Let $\mathcal{D}_n^*=(X_1^*,\dots, X_k^*)$ denote a subsample of size $k$ from the original data $\mathcal{D}_n$ and define $w_j^* =  \mathbf{1}_{X_j\in\mathcal{D}_n^* }$.  Then
	\[\mathrm{Cov}_*(s^*, w_j^*) =  \frac{k}{n}(e_j-s_0),\]
	where $*$ refers to the subsampling procedure, $e_j = \mathbb{E}_*[s^*|X_1^*=X_j]$ and $s_0 = \mathbb{E}_*[s^*]$.
	Thus,
	\begin{equation}
		\mathrm{ps}\text{-}\mathrm{IJ} _\mathrm{U} =  \sum_j \mathrm{Cov}^2_*(s^*, w_j^*)= \left(\frac{{k\choose 1}}{{n\choose 1}}\right)^2\sum_j  (e_j-s_0)^2.
	\end{equation}
\end{proposition}

\noindent In an analogous fashion to $\mathrm{IJ}_\mathrm{U}$, we can establish the following theorem.
\begin{theorem}
	\label{U-sij}
	The pseudo-IJ estimator of the variance of a U-statistic is defined as
	\begin{equation}
		\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}= \frac{k^2}{n^2}\sum_{j=1}^n[ e_j-s_0]^2
	\end{equation}
	where $e_j = \mathbb{E}_*[s^*|X_1^*=X_j]$ and $s_0=\mathbb{E}_*[s^*]$. Then
	\begin{equation}
		\mathbb{E}[\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}] = \sum_{j=1}^k r_j{k\choose j}^2{n\choose j}^{-1} V_j,
	\end{equation}
	where
	\begin{equation}
		r_j = \left(\frac{n-k}{n}\right)^2 \frac{j}{1-j/n}, \quad \textit{for } i = 1,\dots, k.
	\end{equation}
\end{theorem}

Note that although our goal here is to use $\frac{k^2}{n}\hat{V}_1$ to estimate $\frac{k^2}{n}V_1$, $\mathbb{E}[\frac{k^2}{n}\hat{V}_1] = \mathbb{E}[\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}]$ involves higher order terms of $V_2,\dots, V_k$.  Though not ideal, it is unavoidable since we do not have new data generated from the underlying distribution. If we simply multiply $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$ by $(\frac{n}{n-k})^2\frac{n-1}{n}$ as proposed in \cite{Wager2018}, then only the first term is unbiased, but it doubles the quadratic term, triples the cubic term etc. (a similar phenomenon was discovered by \cite{efron1981Nonparametric} for the jackknife variance estimator). This explains why this estimation procedure is inflated in practice.  In many applications, $k$ is not small and so the higher order terms of $\mathrm{Var}(\mathrm{U})$ are not negligible and the effect of $r_j$ cannot be ignored. Indeed, as alluded to at the beginning of this section, in many modern machine learning applications like random forests, $k$ corresponds to the number subsamples utilized in the construction of each base learner (trees, in the case of random forests),and so $k$ is generally best chosen to be as large as possible.

\subsection{The Consistency of $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$ and Its Derivatives}

Comparing Theorem's \ref{U-ij-1} and \ref{U-sij}, it is straightforward to see that if $k/n\to 0$, then $\alpha \to 1$ and $\alpha-\beta\to 0$, and thus $\mathrm{IJ}_\mathrm{U} \to \mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$.  More generally, $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$ has a slighter simpler expression and in this subsection, we investigate its statistical properties, focusing in particular on its consistency.  We begin by introducing the idea of generalized U-statistics, recently defined in \cite{peng2019asymptotic}.
\begin{definition}[Generalized U-statistic \cite{peng2019asymptotic}]
	\label{gu}
	Suppose $X_1, \dots, X_n$ are i.i.d.\ samples from $\mathrm{P}$ and let $s$ denote a (possibly randomized) real-valued function that is permutation symmetric in its $k \leq n$ arguments. A generalized U-statistic with kernel $s$ of order (rank) $k$ refers to any estimator of the form
	\begin{equation}
		\label{eqn:gu}
		\mathrm{U}_{n,k,N,\omega} = \frac{1}{\hat{N}} \sum_{(n,k)} \rho s(X_{i_1}, \dots, X_{i_k}; \omega)
	\end{equation}
	where $\omega$ denotes i.i.d.\ randomness, independent of the original data. The $\rho$ denotes i.i.d.\ Bernoulli random variables determining which subsamples are selected where $\mathbb{P}(\rho=1) = N/{n\choose k}$ and $\hat{N}$ corresponds to the sum of the $\rho$. When $N={n\choose k}$, the estimator in  \cref{eqn:gu} is a generalized complete U-statistic and is denoted as $\mathrm{U}_{n,k,\omega}$.  When $N < {n\choose k}$, these estimators are generalized incomplete U-statistics.
\end{definition}
Generalized U-statistics are essentially incomplete U-statistics with potentially extra randomness and where the order of the kernel may grow with $n$.  Random forests, for example, are generalized U-statistics in which the additional randomness $\omega$ determines which features are eligible for splitting at each node in each tree.
In recent work, \cite{peng2019asymptotic} proved that if $\frac{k}{n}({\zeta_k}/{k\zeta_{1,\omega}}-1)\to 0$, then the complete generalized U-statistic $\mathrm{U}_{n,k,\omega}$ is asymptoticly normal with variance $\frac{k^2}{n}\zeta_{1,\omega}$, where $\zeta_k= \mathrm{Var}(s)$ and $\zeta_{1,\omega} =V_1=  \mathrm{Var}(\mathbb{E}[s|X_1])$. Let $e^{\omega}_i = {n-1\choose k-1}^{-1} \sum s(X_i, \dots;\omega)$ and let $s_0^{\omega} = {n \choose k}^{-1}\sum s(\dots; \omega)$. Note that each collection of subsamples is paired with an i.i.d. $\omega$. Then the corresponding $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$ for generalized U-statistics is defined as
\begin{equation}
	\label{eq:psIJomega}
	\mathrm{ps}\text{-}\mathrm{IJ}^{\omega}_\mathrm{U} = \frac{k^2}{n^2}\sum [e_i^\omega - s_0^\omega]^2.
\end{equation}
The following theorem gives that the same conditions used by \cite{peng2019asymptotic} to establish asymptotic normality are sufficient to establish the consistency of $\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}$.  In other words, if the (generalized) U-statistic is nearly linear, then  $\mathrm{ps}\text{-}\mathrm{IJ} ^\omega_\mathrm{U}/ \mathrm{Var}(\mathrm{U}_{n,k,\omega})$ converges to 1 in probability.

\begin{theorem}
	\label{U-1}
	Let $X_1,\dots, X_n$ be i.i.d.\ from $\mathrm{P}$ and  $\mathrm{U}_{n,k,\omega}$ be a generalized complete U-statistic with kernel $s(X_1,\dots,X_{k};\omega)$. Let $\theta=\mathbb{E}[s]$, $\zeta_{1,\omega}=\mathrm{Var}(\mathbb{E}[s|X_1])$ and $\zeta_k=\mathrm{Var}(s)$.
	If $\frac{k}{n}(\frac{\zeta_k}{k\zeta_{1,\omega}}-1)\to 0 $, then
	\begin{equation}
		\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U} /\mathrm{Var}(\mathrm{U}_{n,k, \omega})\xrightarrow[]{p}1.
	\end{equation}
\end{theorem}

\begin{corollary}
	If the conditions in \cref{U-1} are met, then
	\begin{equation}
		\begin{aligned}
			\label{eq:corpsij}
			\mathrm{U}_{n,k,\omega} \pm z_{\alpha/2}\frac{n}{n-k}\sqrt{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}}
			& =  \mathrm{U}_{n,k,\omega} \pm z_{\alpha/2}\frac{n}{n-k}\sqrt{\sum \mathrm{Cov}^{\omega}_*(s^*,w_i^*)^2}\\
			& =   \mathrm{U}_{n,k,\omega} \pm z_{\alpha/2}\frac{k}{n-k}\sqrt{\sum(e^\omega_i-s^\omega_0)^2}
		\end{aligned}
	\end{equation}
	provides an asymptotically valid confidence interval for $\theta$ with confidence level $1-\alpha$. Note that $\mathrm{Cov}^\omega_*(s^*, w_i^*)$ can be defined in the same fashion as \cref{eq:psIJomega} to include the extra randomness and equals $\frac{k}{n}(e_i^\omega-s_0^\omega)$.
\end{corollary}

\noindent The $\frac{n}{n-k}$ appearing in \cref{eq:corpsij} is there to correcting for finite sample bias.

Because of their inherent computational burden in calculating ${n\choose k}$ base estimates, the complete forms of these estimators are almost never utilized in practice.  Fortunately, asymptotic normality for incomplete generalized U-statistics was also established in \cite{peng2019asymptotic}. Thus, to establish asymptotically valid confidence intervals for these incomplete counterparts, we need only establish a consistent means of estimating the asymptotic variance of $\mathrm{U}_{n,k,N,\omega}$.

Let
\begin{equation}
	\begin{aligned}
		\label{sij-dagger}
		\widehat{\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}^\omega} &= \frac{k^2}{n^2}\sum_i [\hat{e}^\omega_i - \hat{s}^\omega_0]^2,
	\end{aligned}
\end{equation}
where
\begin{equation}
	\hat{s}^\omega_0 = \frac{1}{N}\sum s(...;\omega), \quad \hat{e}^\omega_i = \frac{n}{Nk}\sum s(X_i,...;\omega).
\end{equation}
Here, $\sum s(...;\omega)$ denotes the sum of all kernels that make up the incomplete U-statistic, whereas $\sum s(X_i,...;\omega)$ denotes the sum of all kernels that make up the incomplete U-statistic and include $X_i$ in their respective subsamples. The following theorem shows that we can obtain a consistent estimate of $\zeta_{1,\omega}$ so long as $N$ is large enough to ensure that $\frac{n}{Nk\zeta_{1,\omega}}\to 0$.
Importantly, this means that $N$ need not be on the order of ${n\choose k}$.

\begin{theorem}
	\label{U-2}
	Let $X_1,\dots,X_n$ be i.i.d.\ from $\mathrm{P}$ and $\widehat{\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}^\omega}$ be as  \cref{sij-dagger}. Let  $\zeta_{1,\omega}=\mathrm{Var}(\mathbb{E}[s|X_1])$ and  $\zeta_k=\mathrm{Var}(s)$. Assume that $\mathbb{E}[s]$ and $\zeta_k$ are bounded.
	Then if $\frac{k}{n} (\frac{\zeta_k}{k\zeta_{1,\omega}}-1)  \to  0 $ and
	$\frac{n}{Nk\zeta_{1,\omega}}\to 0 $,
	we have
	\begin{equation}
		\widehat{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}} / \frac{k^2}{n}\zeta_{1,\omega}\xrightarrow[]{p}1.
	\end{equation}
\end{theorem}
\begin{remark}
	Consider the case that $ \zeta_k/k\zeta_{1,\omega}\leq c_1$ for and $k\zeta_k\geq c_2$ for some constants $c_1$ and $c_2$. \cref{U-2} states that we need only have $n\gg k$ and $N\gg nk$ in order to ensure that $\widehat{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}}/ \frac{k^2}{n}\zeta_{1,\omega}  \xrightarrow[]{p}1$.
\end{remark}

In the recent literature on random forests discussed above, $\sum_i \widehat{\mathrm{Cov}}^2(s^*,w_i^*)$ is often the quantity used to estimate the variance predictions in practice.  The estimator $\widehat{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}}$ has a strong connection to this quantity. Indeed,
\begin{equation}
	\begin{aligned}
		\frac{k^2}{n^2}\sum_i [\hat{e}_i^\omega - \hat{s}_0^\omega]^2
		& = \frac{k^2}{n^2}\cdot \frac{n^2}{k^2}\sum [\frac{1}{N}s(...;\omega) w_i^* - \frac{1}{N}\sum s(...;\omega) \frac{k}{n}]^2 \\

		& \approx  \frac{\hat{N}^2}{N^2}  \sum_i \left[ \frac{1}{\hat{N}}\sum s(\dots;\omega)w_i^* - \frac{1}{\hat{N}} \sum s(\dots;\omega)\frac{\hat{N}_i}{\hat{N}}\right]^2\\
		& =  \frac{\hat{N}^2}{N^2} \sum_i \widehat{\mathrm{Cov}}^2(s^*, w_i^*)	\approx \sum_i \widehat{\mathrm{Cov}}^2(s^*,w_i^*).
	\end{aligned}
\end{equation}
where $\hat{N}_i = \sum w_i$. And since $\hat{N}$ concentrates around $N$ and $\hat{N}_i/\hat{N}$ concentrates around $\mathbb{E}_*[w^*_i] =\frac{k}{n}$, $\widehat{\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}^\omega} $ should be close to $\sum_i \widehat{\mathrm{Cov}}^2(s^*,w_i^*)$.

To our knowledge, \cref{U-2} is the first result that provides a consistency of estimate the variance of predictions generated by subsampled random forests, where the number of trees in the random forest need not to be ${n\choose k}$ and the trees themselves may be randomized.

Recent work has established that $\mathrm{U}_{n,k,N,\omega}$ has different asymptotic distributions depending on the number of subsamples $N$ that are employed \cite{peng2019asymptotic}. When $N\ll n/k$, $\mathrm{U}_{n,k,N,\omega} \sim \mathcal{N}(0, \zeta_k/N)$, when $N = O(n/k)$, $\mathrm{U}_{n,k,N,\omega} \sim \mathcal{N}(0, \frac{k^2}{n}\zeta_{1,\omega} + \frac{\zeta_k}{N})$, and when $N \gg n/k$ and $\mathcal{N}(0, \frac{k^2}{n}\zeta_{1,\omega})$.  Thankfully, regardless of the setting, there are only two variance parameters that may need to be estimated:  $\zeta_k$ and $\zeta_{1,\omega}$.  We can estimate  $\zeta_k$ simply by calculating sample variance of base learners built on non-overlapping subsamples. Based on the above arguments, $\zeta_{1,\omega}$ can be estimated by $\widehat{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}}$. Therefore, it is guaranteed that
\begin{equation}
	\widehat{\mathrm{ps}\text{-}\mathrm{IJ}^\omega_\mathrm{U}} \xrightarrow[]{p} \frac{k^2}{n}\zeta_{1,\omega} \quad \text{and}\quad \widehat{\zeta_k}\xrightarrow[]{p} \zeta_k.
\end{equation}

Interestingly, if a random forest is built with $N = O(n)$ decision trees, we cannot estimate the variance of the random forest consistently using only the trees contained in the random forest; this requires $\gg \frac{n}{k\zeta_{1,\omega}}$ decision trees. Note that $k\zeta_{1,\omega}\leq \zeta_k$ and $\zeta_k$ typically tends toward $0$ as $k$ grows. Therefore, the number of trees required for consistent variance estimation is $\gg n$.
These results shed light on the intuition established throughout the machine learning literature that it is always significantly more computationally intensive to estimate the variance of ensembles than to obtain the ensemble itself.

\section{Discussion}
The work above provides an in-depth examination of the infinitesimal jackknife estimate of variance for resampled statistics.  We provided alternative perspectives and derivations of the estimator, perhaps most notably demonstrating its equivalence to an OLS linear regression of the bootstrap estimates on their respective sampling weights.  Ultimately we derived three alternative estimation methods under the bootstrap regime and demonstrated their equivalence when all bootstrap samples are employed.
We also gave an in-depth examination of both the Monte Carlo and sampling bias of the IJ estimator in the bootstrap setting and proposed a novel bias-corrected estimator, ultimately providing conditions related to the asymptotic linearity of the statistic under which the IJ estimator is asymptotically unbiased.
In the latter portion of the work, we examine how these preliminary results translate outside the bootstrap setup by looking at setups where the estimator is formed via subsampling without replacement.  Here the statistics admit a familiar U-statistic structure and we derived corresponding results for generalized U-statistics, now a popular tool for analyzing modern supervised learning ensembles like random forests.  Here we also provided a formal motivation for the pseudo IJ often employed in such settings and establish its consistency under similar linearity conditions.  Importantly, for the first time, we further established consistency for the finite sample (incomplete) versions of these estimators.  In particular, this means that one needn't utilize all possible subsamples in order for the empirical $\mathrm{ps}\text{-}\mathrm{IJ}_\mathrm{U}$ to form a consistent estimator.

Finally, recall from the previous section that in dealing with the U-statistic estimators, we assume that the statistics are approximately linear so that when we write $\mathrm{Var}(\mathrm{U}) =\sum_{j=1}^k {k\choose j}^2{n\choose j}^{-1}V_j$, the variance is dominated by the first order term $k^2/n V_1$ and we propose an estimate of $V_1$ accordingly.  One may wonder whether, if the remaining terms are not negligible, an improved estimate can be provided by also including estimates for  $V_j$ for $j=2,\dots, k$.  While such estimates can be provided, establishing the superiority of the resulting estimator is a far more in depth undertaking that we reserve for future work.  Further detailed discussion is provided in Appendix \ref{app:higherorderIJ}.

\bibliographystyle{Chicago}