\section{Introduction}

In many situations occurring in applied mathematics (for example mathematical models or  numerical simulation), when working with an input-output system with uncertain (random)  inputs, it is crucial to understand the global influence of one of the inputs on the output. This problem is generally called global sensitivity analysis (or in short sensitivity analysis). We refer, for example to
\cite{tarantola2008} and  \cite{saltelli2000} for an overview on practical aspects of sensitivity analysis.  Sensitivity analysis aims to give some quantitative indicator allowing to score the global influence of each input variable of the model. A very popular index well tailored for the case of a real valued output  is the so-called Sobol index early proposed in \cite{sobol1993}.   It is based on second order moments of the distributions related to the input-output system.
Different strategies have been implemented for its statistical estimation (see for instance \cite{gamboa2013_s} or \cite{rugama2018}).  This kind of second order indices may also be considered in the more general  frame where the output is a  real vector or a real function (see for instance \cite{klein2014} or \cite{marrel2011}).  
Nevertheless, the variance-based indices have some drawbacks. Indeed, they only study the impact on the variance of the output and so provide  only restricted summary of the output distribution (see \cite{daveiga2015}). Hence, 
considering higher order indices involving the whole distributions and not only their second order properties may give more accurate information on the system. 
In \cite{pianosi2015},  a sensitivity measure based on the Kolmogorov-Smirnov statistic is proposed and studied while entropy-based sensitivity measures are considered for instance  
in \cite{liu2006r}  and \cite{borgonovo2007}. 
For a general overview we also refer to  \cite{borgonovo2017}.
When the output $Z=F(X)$ ($X=(X_i)$ is the input vector), is valued in a more general space, the authors of  \cite{borgonovo2016} have proposed the following general sensitivity index for the input $X_i$:
\begin{displaymath}
S_i \mathrel{\mathop:}= \mbox{$\mathbb{E}$}_{X_i} \left[ d \left(P_Z, P_{Z \vert X_i} \right) \right].
\end{displaymath} 
Here,  $d(\cdot, \cdot)$ is a given dissimilarity measure between  probability measures, $\mbox{$\mathbb{E}$}_{X_i}(\cdot)$ denotes the expectation with respect to $X_i$ and $P_Z$  (resp. $P_{Z \vert X_i}$) is the unconditional (resp. conditional) probability. Notice that in general, the study of an accurate statistical estimation of such index is not obvious.


Up to our knowledge, the case where the output takes its values on a Riemannian manifold has not yet been studied. 
Beside, in \cite{gamboa2015} a  sensitivity measure in the case of a scalar output and based on the Cram\'er von Mises distance is proposed an studied. This paper extends this last approach to the more general frame of an output valued on a Riemannian manifold.

In the last decades starting from the pioneer 1945' work of Rao \cite{rao1945}, the statistical theory for data valued on a Riemannian manifold have received a lot of interest and contributions. One of the reason of the development of such theory is the spectacular increase of the computation power allowing the treatment of more and more complex objects and structures on computers. References on the subject are numerous. We refer to \cite{mardia1972}, \cite{bhattacharya2012}  and  \cite{ellingson2015} and the references therein for an overview. In this paper, our aim is to bridge this theory to sensitivity analysis. We build a general sensitivity index. On the one hand this index takes into account the whole characteristic of the involved distributions and not only their second order moments.  On the other hand the working frame is a general system with an output valued on a Riemannian manifold. One of the main ingredient of our work is a generalization of the Cram\'er von Mises criterion replacing, before integrating, the half lines by geodesic balls. Consequently, our new index has the nice property to be invariant with respect to the isometric maps of the Riemannian manifold. Moreover, we show that this generalized index can be easily estimated using a Pick Freeze like estimator. This estimator involves a $U$-statistics. Without loss of generality,  we will assume that the random input variables are real and independent. Sensitivity analysis for dependent inputs may be performed but is generally less readable (see for example  \cite{gamboa2015b} and references therein).
  

 The paper is organized as follows. In Section  \ref{conceptosbasicos}, in order to be self contained, we recall some key facts and properties on Riemannian manifolds. In Section \ref{indicedef}, we define our new sensitivity index for  Riemannian manifold. Further, we will build a Pick an Freeze estimator based on a $U$-statistic and study its asymptotic (consistency), and non asymptotic (concentration inequality) properties.    In Section \ref{datosreales}, we study the sensitivity of a real life model coming from mechanic. The output of the system is the rigidness matrix. This matrix is involved to model the linear elasticity of a solid object.  The set of these matrices is as a sub-manifold of the symmetric positive matrices. In the last section, we bold the main advantages of our method and discuss its possible extensions. The codes in  Julia and R languages (see \cite{julia2012} and \cite{team2017} respectively) are available upon request to the authors. All proofs are postponed to the Appendix.

\section{Basic concepts}

\label{conceptosbasicos}
Let us begin with  some useful tools and facts in Riemannian geometry.  A \textit{Riemannian metric} $g$ on a manifold $\mbox{$\mathcal{M}$}$ allows to define for every point $p \in \mbox{$\mathcal{M}$}$  a scalar product $g_p(\cdot,\cdot)$ acting on the tangent space (at $p$), $T_p \mbox{$\mathcal{M}$}$. This scalar product depends smoothly  on $p$.   The \textit{Riemannian manifold} $(\mbox{$\mathcal{M}$},g)$ is the manifold equipped with the Riemannian metric  $g$. 
For any $v \in T_p \mbox{$\mathcal{M}$}$ the Riemannian norm of $v$ is given by $\Vert v \Vert \mathrel{\mathop:}= \sqrt{g_p(v,v)}$.
Let $x,y\in M$ and  $\gamma: I:=[0, 1] \rightarrow \mbox{$\mathcal{M}$}$ be a continuously differentiable curve contained in the Riemannian  manifold $\mbox{$\mathcal{M}$}$ (that is assumed to be connected) with 
$\gamma(0)=x$ and $\gamma(1)=y$. Let denote by  $L(\gamma)$ its length. We may now define the \textit{induced distance}  $d_g(x, y)$ between the points $x$ and $y$ of $\mbox{$\mathcal{M}$}$  by setting
    $$d_g(x,y) = \inf \{ L(\gamma) : \gamma \, 
 \}.$$
$d_g$ is the Riemannian distance on $\mbox{$\mathcal{M}$}$  with respect to the $g$ metric. A \textit{geodesic} (with speed $s \in \mathbb{R}^{+}_0$) is a smooth map $\alpha:I \rightarrow \mbox{$\mathcal{M}$}$, such that $\Vert \alpha^{'}(t) \Vert= s$  for all $t \in I$ and which is {\it locally length minimizing}. For $p \in \mbox{$\mathcal{M}$}$ and  $v\in T_p \mbox{$\mathcal{M}$}$, there exists a unique geodesic  $\alpha_{(p,v)}(t)$ starting from that point with initial tangent vector $v$. The \textit{exponential map} is the map $\exp_p$ given by $\exp_p (v) := \alpha_{(p,v)}(1)$. 


Notice further that the geodesic that joins two points is not necessarily unique.
The \textit{cut locus}  of $p$ in the tangent space is defined to be the set of all vectors $v \in T_{p} \mbox{$\mathcal{M}$}$ such that  $\exp _{p}(tv)$ is a minimizing geodesic for  $0 \leq t \leq 1$  but fails to be one for $0 \leq t \leq 1+ \epsilon$, $\epsilon >0$. The \textit{cut locus of $p$ in $\mbox{$\mathcal{M}$}$}, denoted  by $C_{\mbox{$\mathcal{M}$}}(p)$, is defined as the image of the cut locus of  $p$ in the tangent space under the exponential map at $p$. The \textit{injectivity radius of $p$} is the maximal radius of centered balls on which the exponential map is a diffeomorphism. The \textit{injectivity radius of the manifold}  $r_{iny}$ is the infimum of the injectivity over the manifold. For example, in the sphere $S_{d}$ of $\mbox{$\mathbb{R}$}^{d+1}$ the cut locus of a point $p$ is its antipodal point  $-p$. There are infinite minimizing geodesics that connect a point with its antipodal (Figure  \ref{bola_diametro} shows two geodesics) but this configuration has zero probability of being drawn when dealing with the uniform probability measure. 

\begin{figure}[!ht]
\centering
\subfigure{\includegraphics[width=50mm]{bola_diametro3.png}}
\subfigure{\includegraphics[width=40mm]{asas2.pdf}}
\caption{\textit{Left Panel}: Representation of a geodesic ball of diameter $\overline{PQ}$ and center $O$ on a hyperbolic paraboloid. \textit{Right Panel}: A graph of two minimizing geodesics that join two antipodal points of a sphere.}
\label{bola_diametro}
\end{figure}


Let $(\mbox{$\mathcal{M}$},g)$ be a connected and orientable Riemannian manifold,(see \cite{docarmo1992}, page $18$). We will assume that  $(\mbox{$\mathcal{M}$},d_g)$ is a complete  separable metric space. Since $(\mbox{$\mathcal{M}$},d_g)$ is complete, the Hopf--Rinow theorem (see \cite{docarmo1992}, p.146) implies that for any pair of points $p,q \in \mbox{$\mathcal{M}$}$  there exist at least one geodesic path in $\mbox{$\mathcal{M}$}$ connecting $p$ and $q$. 
If the manifold $\mbox{$\mathcal{M}$}$ fulfills the assumptions of the Cartan--Hadamard Theorem (see \cite{petersen2006}, p.162)--that is, if $\mbox{$\mathcal{M}$}$ is a simply connected, complete Riemannian manifold with non--positive curvature (a Hadamard manifold), then the geodesic is unique.  
 Let $X$ be a random element taking values in $\mbox{$\mathcal{M}$}$, with distribution $P$. In many useful examples (eg. the sphere),  it occurs that geodesic uniqueness fails.  But in what follows, we will only need that this failure not occurs too often. Roughly speaking, we will assume a condition on $P$ ensuring that there is a unique geodesic between each pair of points $p,q \in \mbox{$\mathcal{M}$}$ with  probability one. More precisely, we will assume that the random element $Y$ has a density $f_Y$ with respect to the volume measure $d \nu(y)$ on $\mbox{$\mathcal{M}$}$, (which exists since $\mbox{$\mathcal{M}$}$ is orientable - see  the last section of \cite{folland2013}) fulfilling the following condition:
Given $q \in \mbox{$\mathcal{M}$}$ let 
\begin{equation}
A_q\mathrel{\mathop:}= \left\{y \in \mbox{$\mathcal{M}$}  \Big/
\begin{array}{ll}
\textrm{ there are more than one different } \\ 
\textrm{minimizing geodesic connecting $y$ and $q$}
\end{array} \right \}.
\end{equation}
Then, there exist a Borel set $B_q \subset \mbox{$\mathcal{M}$}$ with $A_q \subset B_q$ and such that for any  $q \in \mbox{$\mathcal{M}$}$
\begin{equation}
\label{unicidadd}
\int_{B_q} f_Y(y) d \nu(y)=0.
\end{equation}
\begin{remark}
Let $C_{\mbox{$\mathcal{M}$}}(p)$ stand for the  cut locus of $p$ in a complete manifold $\mbox{$\mathcal{M}$}$, (see \cite{docarmo1992}, p. $267$).  If $P \left( Y \in C_{\mbox{$\mathcal{M}$}}(p)\right) = 0$,  for all   $p \in \mbox{$\mathcal{M}$}$, then condition  (\ref{unicidadd}) is fulfilled (see \cite{pennec2006}). For instance, if the Riemannian manifold is Hadamard, then the cut locus of any point $p$ is empty and the condition is obviously fulfilled. In the unit sphere $S_{d} \mathrel{\mathop:}= \{u \in \mathbb{R}^d / \Vert u \Vert= 1 \}$  the cut locus of a point $p$ is its opposite point  $-p$ and any  probability measure having a density supported by the sphere fulfills condition  (\ref{unicidadd}).
\end{remark}
In this paper,  we will assume that the Riemannian manifold $(\mbox{$\mathcal{M}$},g)$ with the induced distance $d_g$ is connected and oriented, and that the metric space $(\mbox{$\mathcal{M}$},d_g)$ is separable and complete. We will also assume  that given two points $p,q \in \mbox{$\mathcal{M}$}$  there is a unique geodesic determined by $p,q$ with probability one with respect to the tensorial probability measure $d\xi(p,q) \mathrel{\mathop:}= f_Y(p)f_Y(q)d \nu(p)d\nu(q)$, (see Figure \ref{geodesica}).

\begin{figure}[!ht]
\centering
\subfigure{\includegraphics[width=38mm]{cara.png}}
\subfigure{\includegraphics[width=50mm]{geodesica_dona.png}}
\caption{\textit{Left Panel}: A graph of a geodesic that determines two points on the representation
of a face. \textit{Right Panel}: A graph of a geodesic that determines two points on the representation
of a Torus.}
\label{geodesica}
\end{figure}


\subsection{Geodesic balls with diameter $\overline{pq}$}
\label{bolasgeodesicas}
For any pair $p,q \in \mbox{$\mathcal{M}$}$ that determine a unique geodesic   $\overline{pq}$,  we define the ball of diameter  $\overline{pq}$  as the closed ball whose center is the middle point of the geodesic joining $p$ and $q$ with radius $d_g(p,q)/2$. It will be denoted by $B_{pq}$. Figure  \ref{bola_diametro}  depicts  a geodesic ball of diameter $p$ and $q$, where $o$ is the midpoint of the geodesic.
In \cite{christensen1970} it is shown that the family  $\mathcal B_p:= \{{B}_{p}, p \in \mbox{$\mathcal{M}$}\}$ of geodesic balls with center $p$ is a determining class; that is, if two probability measures $\eta, \nu$ on $\mbox{$\mathcal{M}$}$  coincide on any event of this class, then $\eta=\nu$.
We assume that $\mbox{$\mathcal{M}$}$ is a compact Riemannian manifold. We will start by proving that the family of balls $$\mathcal{B}_{pq}:= \{B_{pq}, p,q \in \mbox{$\mathcal{M}$}\ \textrm{and there is a unique geodesic between $p$ and $q$} \}$$ is also a determining class if the radius of injectivity  $r_{iny}$ of $\mbox{$\mathcal{M}$}$ is positive.  We will need the following weaker assumption:

\textbf{HB) B-continuity} A probability measure $\nu$ defined on a Riemannian manifold $\mbox{$\mathcal{M}$}$ fulfils \bf HB \rm if $\nu(\partial A )=0$  for all closed Borel sets  $A$ on $\mbox{$\mathcal{M}$}$ (Here $\partial A$ denotes the topological boundary of the set $A$).
 
 
\begin{property}
\label{propiedad1}
Let $(\mbox{$\mathcal{M}$},g)$ be a compact Riemannian manifold and $\nu$ be a probability distribution fulfilling \textbf{HB}, such that $\nu \left( \partial \mbox{$\mathcal{M}$} \right)=0$. If the injectivity radius is positive  then  $\mathcal{B}_{pq}$ is a determining class of $\nu$.
\end{property}
 In  \cite{fraiman2017} we showed what if $K \in \mbox{$\mathcal{M}$}$ is a compact set, then the family of balls 
$\mathcal{B}_{pq}$ is also a Glivenko--Cantelli class in $K$; that is, $$\sup_{p,q \in K}  \vert  P \left(B_{pq} \right) - P_n\left(B_{pq} \right) \vert 
 \rightarrow 0  \quad \textrm{a.s. \ \  as  $n \rightarrow + \infty$}.$$

In the next section, we define a sensitivity index using this family of sets.


\section{Sensitivity Index in Geodesic Balls}
\label{indicedef}
In order to build a sensitivity index, we use an idea previously developed in \cite{gamboa2015}. In our frame we  compare the distribution $P_Y$ and the conditional distribution $P_{Y \vert X_i}$ on the family of geodesic balls on the manifold (instead of half--lines in \cite{gamboa2015}. So that, we build an intrinsic index related to the manifold. This index depends only on the dimension and structure of the manifold and not on the dimension of the space where it is immersed, unlike in \cite{klein2014} and \cite{gamboa2015}. 


\subsection{Constructing the index}

Let $\mathbb{X}=(X_1, \ldots, X_d) \in \mbox{$\mathbb{R}$}^d$ be an random vector and assume that $\mathbb{P} \mathrel{\mathop:}=  \mathbb{P}_1 \times \ldots \times \mathbb{P}_d$ is the probability law of $\mathbb{X}$. Further, let
\begin{displaymath}
Z=f(X_1, \ldots, X_d),
\end{displaymath}
\noindent where $f: \mathbb{R}^d \rightarrow  \mbox{$\mathcal{M}$}$ is a continuous function. Here, $\mbox{$\mathcal{M}$}$ a Riemannian manifold of dimension $k$. We wish to understand how sensitive is the  output $Z$ to perturbations in some of the input variables $(X_1,\ldots, X_d)$. As any manifold can be immersed in $\mathbb{R}^p$ for a sufficiently large $p$ (see \cite{nash1954}), this problem could be considered as a particular case of a multivariate output. Nevertheless,  in such a way the geometry and \textit{``minor dimensionality''} of the output are not taken into account. We begin by studying a very particular case and then extend our results to any Riemannian manifold, that satisfies   the previous conditions. Let $h$  a measurable function such that $h(X_1,X_2)$ is integrable, we set $$\mbox{$\mathbb{E}$}_{X_2}\left(h(X_1,X_2) \right) \mathrel{\mathop:}= \mbox{$\mathbb{E}$} \left( h(X_1,X_2)  / X_1 \right).$$




\subsection{A very particular case: The real line.}

Let us first focus on the very particular case  $\mbox{$\mathcal{M}$}=\mbox{$\mathbb{R}$}$. Let $F$ be the distribution function of $Z$, $F(t) \mathrel{\mathop:}= P(Z \leq t) \, (t\in \mathbb{R})$, and $F^{\nu}$ be  the distribution of $Z$ conditioned by a subset of the variable $(X_1, \ldots, X_d)$. That is, let $\nu= \{i_1,\ldots, i_k\}  \subset \{1, \dots, d\}$, and


\begin{align*}
& \hspace{-0.5cm} X_{\nu}\mathrel{\mathop:}= (X_{i_1}, \ldots, X_{i_k}),\\
&\hspace{-0.5cm} F^{\nu}(t)\mathrel{\mathop:}=   \mbox{$\mathbb{E}$} \left( \mathbbm{1}_{(\infty,t]} (Z) \big \vert X_{\nu} \right).
\end{align*}

\noindent In \cite{gamboa2015} the following Cram\'er von Mises sensitivity index is considered. Assume that $F$ is absolutely continuous. This normalized index (denoted  by  $C^{\nu}_{2}$) is defined as 

\begin{equation}
C^{\nu}_{2} \mathrel{\mathop:}= \frac{N^{\nu}_{2,}}{\int_{\mathbb R} F(x)\left(1-F(x) \right) dF(x)} = 6 N^{\nu}_{2},
\end{equation}

where
\begin{align*}
N^{\nu}_{2} =\int  \mbox{$\mathbb{E}$} \left( \left[  F(t) -F^{\nu}(t) \right ]^2 \right) dF(t) \\
&\hspace{-4.7cm}= \int \left ( \mbox{$\mathbb{E}$}  \left[   \int \mathbbm{1}_{(\infty,t]}(z)  d(F -F^{\nu})(z) \right ]^2\right )  dF(t)\\
&\hspace{-4.7cm}= \mbox{$\mathbb{E}$} \left\{ \mbox{$\mathbb{E}$}  \left(  \left[ \mbox{$\mathbb{E}$} \left(\mathbbm{1}_{(\infty,t]}(Z) \right) - \mbox{$\mathbb{E}$} \left(\mathbbm{1}_{(\infty,t]}(Z) \big \vert X_{\nu} \right)  \right]^2 \right) \right\}.
\end{align*}

\noindent To build  of our new index, we replace the  function $\mathbbm{1}_{(\infty,t]}$  by the indicator function of the interval $[\min(s,t), \max(s,t)]$ (denoted by $h_{s,t}$),

\begin{displaymath}
h_{s,t}(x) \mathrel{\mathop:}= \mathbbm{1}_{ \{s \leq x \leq t\} }   +   \mathbbm{1}_{\{t \leq x \leq s\}}=   \mathbbm{1}_{ \{\min\{s,t\} \leq x \leq \max\{s,t\}\}}.  
\end{displaymath}





\noindent The function $h_{s,t}$, can be thought of as the indicator function of the ball of diameter $\overline{st}$ in $\mathbb{R}$. Further, let

\begin{equation}
\label{bla}
H(s,t)\mathrel{\mathop:}= \mbox{$\mathbb{E}$}\left[ h_{s,t} (Z)\right] \quad \textrm{and}  \quad  H^{\nu}(s,t)= \mbox{$\mathbb{E}$}\left[ h_{s,t} (Z) \Big \vert X_{\nu} \right]. 
\end{equation}
\noindent  Obviously, $\mbox{$\mathbb{E}$}_{X_{\nu}}\left[ H^{\nu}(s,t) \right]= H(s,t)$. 

\begin{definition}
\label{sensibilidad}
The normalized ball sensitivity index (denoted by $B^{\nu}_{2}$) is then defined as,
  \begin{equation}
  B^{\nu}_{2} \mathrel{\mathop:}= \frac{S^{\nu}_{2}}{\int_{\mathbb R^2} H(x,y)\left(1-H(x,y) \right) dF(y)dF(x)} = 6 S^{\nu}_{2},
  \end{equation}
  where 
  \begin{equation}
  S^{\nu}_{2} \mathrel{\mathop:}= \mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ \mbox{$\mathbb{E}$}_{X_{\nu}} \left \{  \left[ H(Z_1,Z_2) - H^{\nu}(Z_1,Z_2)\right]^2 \right \} \right],
  \end{equation}
\noindent  and $Z_1$ and $Z_2$ are two independent copies of $Z$.
(The constant 6 follows from the well known fact that $\int_{t>s} (t-s)(1-(t-s))dtds = \int_{t\leq s} (t-s)(1-(t-s))dtds=1/12$.)

  
  \end{definition}

We may can rewrite $S^{\nu}_{2}$ as

\begin{displaymath}
S^{\nu}_{2}=\int_{ \mbox{$\mathcal{M}$} \times \mbox{$\mathcal{M}$}}  \mbox{$\mathbb{E}$}_{X_{\nu}} \left \{  \left( H(z_1, z_{2})-  H^{\nu}(z_1,z_{2})\right)^2 \right \} dF(z_1) dF(z_{2}).
\end{displaymath}




\subsection{Generalization for a  Riemannian manifold}

As discussed before, $\mbox{$\mathcal{M}$}$ is a Riemannian manifold of dimension $k$. So that, under some regularity conditions, given two points  $\{z_1, z_2\} \subset \mbox{$\mathcal{M}$}$, the following  function $h_{z_1, z_2 }: \mbox{$\mathcal{M}$} \rightarrow \{0,1\}$ is generically well defined,

\begin{displaymath}
h_{ z_1, z_2  }(t) \mathrel{\mathop:}= \mathbbm{1}_{ B_{z_1z_2}} (t).
\end{displaymath} 


\noindent Let $Z_1,Z_2$ be independent copies of $Z$. Our normalized  sensitivity index is built by means of balls as in Definition \ref{sensibilidad}, 

 \begin{equation}
 \label{indic}
  B^{\nu}_{2} \mathrel{\mathop:}= \frac{S^{\nu}_{2}}{D^{\nu}_{2}},
  \end{equation}
where 
\begin{displaymath}
S^{\nu}_{2} \mathrel{\mathop:}=  \mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ \mbox{Var}_{X_{\nu}} \left \{ \mbox{$\mathbb{E}$}_Z \left[ h_{Z_1,Z_2}(Z) \Big \vert X_{\nu} \right] \right\} \right],
\end{displaymath}


and

\begin{displaymath}
D^{\nu}_{2} \mathrel{\mathop:}= \mbox{$\mathbb{E}$} \left[ H(Z_1,Z_2) \left(1-H(Z_1,Z_2) \right) \right].
\end{displaymath}

\noindent  Notice that the index defined in the previous subsection is a particular case with  $\mbox{$\mathcal{M}$}=\mbox{$\mathbb{R}$}$.

\begin{remark}
 If   $B^{\nu}_{2}=0$ we have that  $\mbox{$\mathbb{E}$}_Z \left[ h_{z_1,z_2}(Z) \Big \vert X_{\nu} \right]= \mbox{$\mathbb{E}$}_Z \left[ h_{z_1,z_2}(Z) \right] \, a.s.$ for all $(z_1,z_2) \in \Omega \subset \mathcal{M}\times\mathcal{M}$ with $P(\Omega)=1$. Therefore, under the assumptions of Property \ref{propiedad1},  the probability measures  $P_{Z \vert X_{\nu}}$ and $P_{Z}$ are the same. 
\end{remark}

\subsection{Estimation}

In the next subsection, an estimator of the index $B^{\nu}_{2}$ is proposed.   The variance of the conditional mean is estimated using the pick and freeze method (see \cite{sobol1993} and \cite{sobol2001}). Further, as the expected value  $\mbox{$\mathbb{E}$}_{Z_1,Z_2}(\cdot)$  is a symmetric function of $Z_1$ and $Z_2$ we will use an $U$--statistic. At the end  of the next subsection, the steps for the construction of the estimator will be detailed.

\subsubsection{Estimation by \textit{``Pick and Freeze''} method}
The method consists in writing the variance of the conditional mean as a covariance. For $\nu \in \{1, \ldots, d \}$, let $X^{\nu}$ be the random vector that coincides with $X$ on its $\nu$ components and is independently regenerated on the others components. That is, $X^{\nu}_{\nu}=  X_{\nu}$  and $X^{\nu}_{\bar{\nu}}=  X^{'}_{\bar{\nu}}$, where $X^{'}$ is an independent copy of $X$ and $\bar{\nu}$  is the complementary set of $\nu$ ($\bar{\nu}= \{1, \ldots, d \} \setminus \nu $). We set
\begin{displaymath}
Z^{\nu} \mathrel{\mathop:}= f(X^{\nu}).
\end{displaymath}


\noindent For sake of simplicity assume first that  $\mbox{$\mathbb{E}$}(Z)=0$. Then,

\begin{equation}
\mbox{Var}\left(  \mbox{$\mathbb{E}$}(Z \big \vert X_{\nu}) \right)= \mbox{$\mathbb{E}$} \left(  \mbox{$\mathbb{E}$}^2(Z \big \vert X_{\nu}) \right).
\label{pick}
\end{equation}

\noindent Since $Z$ and $Z^{\nu}$ are conditionally independent of $X_{\nu}$ (see \cite{klein2013}), so that

 \begin{align*}
 \mbox{cov} (Z, Z^{\nu})=  \mbox{$\mathbb{E}$}(Z Z^{\nu}) - \mbox{$\mathbb{E}$}(Z)\mbox{$\mathbb{E}$}( Z^{\nu})=  \mbox{$\mathbb{E}$}(Z Z^{\nu}) - \mbox{$\mathbb{E}$}^2(Z)\\
 &\hspace{-6.95cm} = \mbox{$\mathbb{E}$} \left[ \mbox{$\mathbb{E}$}(Z Z^{\nu} \big \vert X_{\nu}) \right] = \mbox{$\mathbb{E}$} \left[ \mbox{$\mathbb{E}$}( Z \big \vert X_{\nu})  \mbox{$\mathbb{E}$}( Z^{\nu} \big \vert X_{\nu}) \right] \\
  &\hspace{-6.95cm} =  \mbox{$\mathbb{E}$} \left[  \mbox{$\mathbb{E}$}^2  (Z \vert X_{\nu}) \right].
 \end{align*}

\noindent  Following  \cite{klein2013}, estimating the covariance  by a Monte Carlo Method we obtain the following  Pick Freeze estimator of  $\mbox{Var} \left[  \mbox{$\mathbb{E}$}  (Z \vert X_{\nu}) \right]$,

\begin{equation}
T^{\nu} \mathrel{\mathop:}= \frac{1}{N}\sum_{j=1}^N Z_j Z_j^{\nu} - \left( \frac{1}{2N} \sum_{j=1}^N (Z_j + Z_j^{\nu}) \right)^2.
\end{equation}

\noindent Here, $Z_j$ and $Z^{\nu}_j$,  are $N$ independent copies of $Z$  and  $Z^{\nu}$ respectively.

\vspace{0.2 cm}
\underline{\emph{First step}.} Using the previous idea, we build  a consistent estimator of the numerator $S^{\nu}_{2}$. 
Let $\mathcal{P}_{N,p}$ be the set of all possible subsets of $\{1, \ldots,N\}$ with $p$ elements, that is,

\begin{displaymath}
\mathcal{P}_{N,p}= \{(i_1, \ldots, i_p) \in \{1, \ldots, N \}^p / i_1< \ldots < i_p  \}.
\end{displaymath}


\noindent Let $\tau=(i_1, \ldots, i_p) \in \mathcal{P}_{N,p}$. We denote by $\textbf{W}_{\tau}= \left( W_{i_1},\ldots ,W_{i_{p}} \right)$ a sample of $p$ independent copy of $Z$.
The estimator is obtained as follows,

\noindent  Further,  $D^{\nu}_{2}$  is  estimated analogously setting


\begin{align*}
\hat{D}^{\nu}_{2} \mathrel{\mathop:}= \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2}} \Big \{ \frac{1}{2N} \sum_{j=1}^N   \left( h_{\textbf{W}_{\tau}} \left( Z_j\right)+ h_{\textbf{W}_{\tau}} \left( Z_j^{\nu}\right) \right)  -  \\
& \hspace{-7.5cm}  - \left( \frac{1}{2N} \sum_{i=1}^N \left[ h_{\textbf{W}_{\tau}} \left( Z_i\right)+ h_{\textbf{W}_{\tau}} \left( Z_i^{\nu}\right) \right] \right)^2   \Big \}.
\end{align*}

\noindent The computation  of the  estimator is simple. Indeed, it only involves computation of indicator functions. Nevertheless, the computational time to perform  this estimator could be high. Indeed,  the number of sums is of  order $N^3$. Each term involves the determination of a geodesic. In the following section, we give some asymptotic properties of this estimator. We provide an exponential inequality that leads, from the Borel--Cantelli Lemma, to strong consistency.


\subsection{Asymptotic properties of  $\hat{B}^{\nu}_{2}$.}
We analyzed separately the strong convergence of $\hat{S}^{\nu}_{2}$ and $\hat{D}^{\nu}_{2}$. The strong convergence of  $\hat{S}^{\nu}_{2}$  follows from the next Theorem. Notice that another proof may be built using the fourth moment and the Rosenthal inequality for $U$--statistic of order $2$ developed in \cite{ibragimov1999}.


\textbf{Theorem (Exponential inequality )}
\label{exponencial}
Let $s>0$, there exist $N_0$ such that if $N>N_0$,
\begin{equation}
P \left( \left \vert \hat{S}^{\nu}_{2} -  S^{\nu}_{2} \right \vert > s   \right) \leq 16 \exp \left\{  -\frac{N \left(\frac{s}{9}\right)^2}{8} \right\}.
 \end{equation}

From the previous Theorem and the Borel--Cantelli Lemma the strong consistency of $\hat{S}^{\nu}_{2}$ follows.

\textbf{Corollary (Consistency of the estimator)}
\label{consistencia}
$\hat{S}^{\nu}_{2}$ is a consistent estimator of $S^{\nu}_{2}$.


Notice that we may analogously show the convergence of the denominator involved in (\ref{indic}).
\section{Simulations}
\label{simulaciones}
This simulation   section is made up on three examples. The first one shows the accuracy of the estimator  when the output is real valued. In the second example the output lies on a circle immersed in $\mathbb{R}^2$. The estimates are compared with those obtained in \cite{gamboa2015}. The last example shows that the index proposed in \cite{gamboa2015} may fail as sensitivity indicator.
 
\subsection{Example 1: Output on the real line}
\label{ejemplo1}
We study here an example where the Sobol index does not provide information on sensitivity, while $D^{\nu}_{2}$  and $B^{\nu}_{2}$ do. Let $Z=\alpha X_1 + X_2$  where $\alpha>0$, and $X_1,X_2$ are independent. Assume further that $X_1 \sim \textrm{Bernoulli}(p)$, $X_2 \sim F$, with $m=\mbox{$\mathbb{E}$}(X_2)$, $\sigma^2=\textrm{Var}(X_2)=\alpha^2 p(1-p)$.

We calculate $B^{1}_{2}$ in the Appendix while  $C^{1}_{2}$ has been calculated in \cite{gamboa2015}. In the case when $X_2 \sim U(0,b)$ with $b= \sqrt{12\alpha^2 p(1-p)}$. We find
\begin{displaymath}
B^{1}_{2}= 12p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}
\begin{displaymath}
C^{1}_{2}= 6p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}

In the Figure \ref{index}, the estimates of both indices are compared for different values of $p \in [0,1]$. Samples sizes $100$, $500$ and $1000$ are considered.  For each $p$, through resampling Bootstrap for $U$--statistics (see \cite{gine1992b}), $95\%$ confidence intervals are also given. In all cases, the mean squared deviation (MSD) of the estimators are computed.
The MSD of $\hat{B}^{1}_{2}$ is slightly smaller than the one for $\hat{C}^{1}_{2}$ (see Table \ref{MSDrecta}). Similar results are obtained when comparing  $\hat{B}^{2}_{2}$  and  $\hat{C}^{2}_{2}$. 


\begin{table}[ht]
\centering
\begin{tabular}[t]{ccc}
\toprule
 Size &   $MSD_{\hat{B}^{1}_{2}}$ &  $MSD_{\hat{C}^{1}_{2}}$ \\
\midrule
$N=100$    & 0.051  & 0.067 \\
$N=500$    & 0.022 & 0.028 \\
$N=1000$   & 0.013  & 0.018  \\
\bottomrule
\end{tabular}
\caption{The MSD of the  $\hat{B}^{1}_{2}$ and $\hat{C}^{1}_{2}$ estimators for sample sizes $N=100,500$ and $1000$.}
\label{MSDrecta}
\end{table}%

We may conclude that the new proposed index has a similar behavior as the Cram\'er von Mises one for  the real line.




 
\begin{figure}[ht]\centering 
\includegraphics[width=.3\textwidth]{toyN100.pdf}
\includegraphics[width=.3\textwidth]{toyN500.pdf}
\includegraphics[width=.3\textwidth]{toyN1000.pdf}
\caption{Cram\'er von Mises indexes and our new indexes for Example \ref{ejemplo1}. In all cases,   ($\textendash \, \, \textendash$) are the true index  and ($\edge$) their estimate. Violet color depicts CVM sensitivity index while orange color depicts Ball sensitivity index. The $95\%$ bootstrap confidence intervals are represented with shading.  Left--hand Panel: $N=100$. Center--hand Panel: $N=500$. Right--hand Panel: $N=1000$.}
\label{index}
\end{figure}



\subsection{Example 2: Output on a simple manifold immersed in $\mathbb{R}^2$} 
\label{ejemplo2}
The case where $Z$ ranges on the unit circle $S_1$ of  $\mathbb{R}^2$ is considered. We assume that the input vector has distribution  $$X=\binom{X_1}{X_2}\sim N \left[ \binom{\mu_1}{\mu_2},  \left(
\begin{array}{cc}
\sigma^2_1 & 0 \\
0 & \sigma^2_2
\end{array} \right)  \right].$$

Here, we study the case when the output $Z$ is the normalized version of $X$

\begin{displaymath}
\label{circu}
Z \mathrel{\mathop:}= \frac{X}{\Vert X \Vert}.
\end{displaymath}

The distribution of $Z$ has been widely studied, we refer for example to \cite{mardia1972}, \cite{kendall1974} and \cite{watson1983}. In Figure \ref{circulo} the  index proposed  in \cite{gamboa2015} and the one studied here are depicted for the last system. These sensitivity indices are computed for $\mu_1 \in [-5,0]$. The other values are setted to  $\mu_2=0$ and $\sigma^2_1=\sigma^2_2=1$. We observe that the variability of $\hat{B}^{\nu}_{2}$ is  smaller that  the one of $\hat{D}^{\nu}_{2}$  for a sample size $N=300$.

\begin{figure}[ht]\centering 
\includegraphics[width=.49\textwidth]{circuloB1.pdf}
\includegraphics[width=.49\textwidth]{circuloB2.pdf}
\caption{Calculation of the indices $\hat{B}^{\nu}_{2}$ and $D^{\nu}_{2}$ in the equation \ref{circu} with $\mu_1 \in [-5,0]$ and $N=300$. For all cases, it is represented using  ($\edge$) the index estimation functions for $N=300$. In violet color depicts CVM sensitivity index and orange color depicts Ball sensitivity index. The $95\%$ bootstrap confidence intervals are represented with shading. Left--hand Panel: $\nu=1$. Right--hand Panel: $\nu=2$.}
\label{circulo}
\end{figure}


\subsection{Example 3: Output on a manifold immersed in $\mathbb{R}^3$ }

Let us now consider a system with output in the following manifold $\mbox{$\mathcal{M}$}$ immersed in $\mathbb{R}^3$


\begin{displaymath}
\mbox{$\mathcal{M}$}= \left \{(x,y,z) \in \mathbb{R}^3 /  xyz=1, \, \, x,y,z >0  \right \}.
\end{displaymath}

More precisely, let $Z=f(X,Y)\mathrel{\mathop:}= (X+Y, \frac{1}{X}, \frac{X}{X+Y})$  where $X$ and $Y$ are iid random variables with Gamma distribution with parameters $(\mu_1,1)$, $(\mu_1>0)$. Since the function    $\mathbbm{1}_{\{z \leq w\}} = 0$ for all $z, w \in M, z\neq w$, the index estimation $\hat{C}^{\nu}_{2}$ based on the left lower quadrants introduced in \cite{gamboa2015}, does not provide any information about the sensitivity. Also the second moment of $1/X$ does not exist, therefore it is not possible to compute the Sobol' index defined in \cite{klein2014}. By varying the parameter $\mu_1$, we calculate the estimation  of the ball sensitivity index $B^{1}_{2}$ and $B^{2}_{2}$ . For this purpose,  $N=1000$ pick-freeze samples were generated for each value of $Z$, $(Z_j, Z^{\nu}_j)$, $j=1,\ldots,N$; and other $1000$ samples of $Z$,    independent of $(Z_j, Z^{\nu}_j)$ for the corresponding $W_k$.
In Figure \ref{ej3} the values of the indices are observed by varying the parameter $\mu_1$ between  $0$ and $5$. It is clear from Figure \ref{ej3} that the new index is able to detect the effect of each variable on the output for different values of the parameter $\mu$. 


\begin{figure}[ht]\centering 
\includegraphics[width=.80\textwidth]{ejemplo3.pdf}
\caption{Calculation of the index $\hat{B}^{\nu}_{2}$ }
\label{ej3}
\end{figure}


\section{Sensitivity for isotropic matrix}
\label{datosreales}
We consider now the stiffness matrix $Z$ of an  isotropic materials as a function of the constants of Lam\'e  $\lambda$ and $\mu$ (see \cite{landau1965}, p. 13). Hence, if the temperature is constant we have 

\begin{equation*}
Z = \left(
\begin{array}{cccccc}
K+ 4\mu/3   & K-2\mu/3  & K-2\mu/3 & 0 & 0 & 0
 \\
K-2\mu/3  & K+ 4\mu/3  & K-2\mu/3 & 0 & 0 & 0 \\
K-2\mu/3  & K-2\mu/3  & K+ 4\mu/3 & 0 & 0 & 0 \\
0&0&0&\mu&0&0 \\
0&0&0&0&\mu&0 \\
0&0&0&0&0& \mu\\
\end{array} \right)
\end{equation*}

\noindent where $K= \lambda +2\mu/3$ is the volumetric modulus. The parameter $\mu$ is called the stiffness modulus, $K$ and $\mu$ take nonnegative values and are are modeled with a distribution supported by $\mathbb{R}^{+}$. The set of stiffness matrices is considered as a sub--manifold of the Riemannian manifold of the symmetric positive-definite matrices with the metric $g$. This manifold is usually denoted by $(\mathbb{P}_d,g)$, we refer to \cite{moakher2005} for more on the subject. Given two matrices $A$ and $B$  there exists an unique geodesic joining $A$ and $B$ given by,


\begin{equation}
\gamma(t) \mathrel{\mathop:}=  A ^{1/2} \left(A ^{-1/2} B A ^{-1/2}  \right) ^ t A ^{1/2}. 
\end{equation}  

\noindent So that, we can calculate the midpoint between $A$ and $B$. We will denote by $A\#B$ this midpoint. We have

\begin{align}
A\#B = A ^{1/2} \left(A ^{-1/2} B A ^{-1/2}  \right) ^{1/2} A ^{1/2} \\
& \hspace{-6.5 cm} d(A,B)= \Vert \textrm{ln}  \left(A ^{-1/2} B A ^{-1/2}  \right) \Vert,
\end{align}
where $\Vert \cdot \Vert$  is the Hilbert--Schmidt norm. Note that the midpoint is simply the geometric mean of the matrices. 

\vspace{2 mm}

We will focus on the case when  $K$ and $\mu$ are independent random variables and on two scenarios:  

\begin{tabular}{lll}
\textbf{Case 1} & $\quad K \sim \gamma( 1/\lambda_K, \lambda_K )$ & $ \mu \sim \gamma(1/\lambda_{\mu}, \lambda_{\mu} )$ \\ 
\textbf{Case 2} & $\quad K \sim \textrm{U}(1-\lambda_K , 1+\lambda_K )$ & $\mu \sim \textrm{U}(1-\lambda_{\mu},1+\lambda_{\mu})$ \\ 
\end{tabular} 
 
Tables \ref{indice_ug1} and \ref{indice_ug2}  show the values of  $\hat{B}^{1}_{2}$ and $\hat{B}^{2}_{2}$  for  various values of $(\lambda_K,\lambda_{\mu})$. The sample size was $N=500$. In both tables we can see that the new index provides significant information regarding the sensibility of $Z$ with respect to the input variables $K$ and $\mu$.  In particular, when the variance of the input variables increases (larger $\lambda_{\mu}$ or $\lambda_{K}$) the value of the index increases, i.e. a a greater effect is observed on the output.





\begin{table}[ht]
\begin{center}
\caption{Index on isotropic matrix with Gamma distributions parameter}
\scalebox{0.8}{
\begin{tabular}{c | cccc|| cccc }
  \toprule[0.6 mm]
  Distribution  &  \multicolumn{4}{c}{Case 1: $B^1$ } & \multicolumn{4}{c}{Case 1: $B^2$ }             \\
\toprule[0.3 mm]
 $ \lambda_{\mu} \Big \backslash \lambda_K$ & \textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} &\textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} \\
  \toprule[0.3 mm]
  \textbf{0.001} & 0.625   & 0.212 & 0.016 & 0.001           & 0.083 & 0.435 & 0.855 & 0.980  \\
  \textbf{0.01} &0.925 & 0.593 & 0.215 & 0.033           & 0.004 & 0.072 & 0.458 & 0.865  \\
 \textbf{0.1} & 0.987 & 0.912 & 0.587 & 0.184          & 0.001 &0.006& 0.137 & 0.518         \\
\textbf{1} & 0.999 & 0.990 & 0.930 & 0.600       & 0.000 & 0.007 & 0.210 & 0.311   \\

  \toprule[0.6 mm]
\end{tabular}}
\label{indice_ug1}
\end{center}
\end{table}


\begin{table}[ht]
\begin{center}
\caption{ Index on isotropic matrix with Uniform distributions parameter}
\scalebox{0.8}{
\begin{tabular}{c | cccc || cccc }
  \toprule[0.6 mm]
  Distribution  &  \multicolumn{4}{c}{Case 1: $B^1$ } & \multicolumn{4}{c}{Case 1: $B^2$ }             \\
\toprule[0.3 mm]
  $ \lambda_{\mu} \Big \backslash \lambda_K$ & \textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} &\textbf{0.001} & \textbf{0.01}  & \textbf{0.1} & \textbf{1} \\
  \toprule[0.3 mm]
  \textbf{0.001} & 0.620   & 0.008  & 0.001 & 0.000          & 0.109 & 0.848 & 0.997 & 1.000  \\
  \textbf{0.01}  &0.989     & 0.623 & 0.003 & 0.001           & 0.018 & 0.092 & 0.849 & 0.997  \\
 \textbf{0.1} & 1.000      & 0.989 & 0.623 & 0.624          & 0.016 &0.016& 0.102 & 0.846        \\
\textbf{1} & 1.000         & 1.000 & 0.990 & 0.987            & 0.015 & 0.016 & 0.100 & 0.211     \\

  \toprule[0.6 mm]
\end{tabular}}
\label{indice_ug2}
\end{center}
\end{table}


\section{Conclusions and Future Work}
\label{conclusiones}



We introduce a new framework to measure the  global sensitivity  for system outputs valued on a Riemannian manifolds. This sensitivity index has nice properties:
\begin{enumerate}
\item It is a generalization of the one proposed by \cite{gamboa2015}.
\item  It is built on the whole distributions and not only on moments.

\item The construction of the index lies on the geometry of the support of the output through the  geodesic distance.

\item A pick-freeze like estimator of the index based on an $U$--statistic is  proposed. It is easy to compute. The consistency of the estimator is shown. 

\item Considering  various simulation scenarios, the desirable properties of the estimator are illustrated. Using a Bootstrap resampling method for $U$--statistics we construct confidence intervals.

\item  The sensitivity index for an isotropic matrix model is studied.

\item In future research we will challenge to study new applications, for instance, the impact of inputs generated by a dynamic system, when the output is supported by a Riemannian Lie  Group. 


\end{enumerate}





\section*{Appendix}
\begin{enumerate}[label=(\Alph*)]
\item \textbf{Proof of Property \ref{propiedad1}.}

In \cite{christensen1970} it is shown that both the sets




\begin{equation}\label{determ}
 \mathcal B \mathrel{\mathop:}= \{{B}_{p,r} \big/ p \in \mbox{$\mathcal{M}$}  \, \textrm{and} \, r>0 \}, \ \mbox{and} \ \mathcal B_{\delta} \mathrel{\mathop:}= \{{B}_{p,r} \big/ p \in \mbox{$\mathcal{M}$}  \, \textrm{and} \, 0<r<\delta \},
 \end{equation}
are determining classes. We aim to  show that the family of balls $$\mathcal{B}_{pq}:= \{B_{pq}, p,q \in \mbox{$\mathcal{M}$}\ \textrm{and there is a unique geodesic between $p$ and $q$} \},$$ is also a determining class if the radius of injectivity  $r_{iny}$ of $\mbox{$\mathcal{M}$}$ is positive.
Recall that in a metric space the family of closed sets $\mathcal{C}$ is a determining class (\cite{billingsley1999}, page $7$). Therefore, it suffices to show that if two probabilities $\nu$ and $\eta$ coincides on $\mathcal{B}_{pq}$, then $\nu(A)=\eta(A)$ for all $A \in \mathcal{C}$. Let  $\partial A$ the topological boundary of $A$ and $\partial A^{\epsilon} \mathrel{\mathop:}= \bigcup_{x \in \partial A} B(x, \epsilon)$ for $\epsilon >0$. So,

$$\nu (A)=  \nu \left( A \setminus (\partial A^{\epsilon} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon} ) \right) + \nu \left( A \cap (\partial A^{\epsilon} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon} ) \right).$$

 
Let  $\epsilon^{*} \mathrel{\mathop:}= \min (\epsilon/2, r_{iny})$.  Since by (\ref{determ}) $\mathcal B_{\epsilon^{*}}$ is a determining class, using the exponential map  we derive that the set of balls in $\mathcal B_{\epsilon^{*}}$ 
that determine the $\nu$-probability of $A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}})$  are in
  $\mathcal{B}_{pq}$. Therefore, $$\nu \left( A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}} ) \right)= \eta \left( A \setminus (\partial A^{\epsilon^{*}} \cup \partial \mbox{$\mathcal{M}$}^{\epsilon^{*}}) \right).$$ 
Finally, from the dominated convergence theorem  and the assumptions, we conclude that
$$\nu (A)=  \eta(A)- \eta(\partial A  \cup \partial \mbox{$\mathcal{M}$}) + \nu(\partial A  \cup \partial \mbox{$\mathcal{M}$}) =\eta(A).$$





\QEDB

\item \textbf{Proof of Property \ref{exponencial}.}


We will show that given $s>0$, there is an $N_0$ such that for any $N>N_0$,
\begin{equation}
\label{exponencial2}
P \left( \left \vert \hat{S}^{\nu}_{2} -  S^{\nu}_{2} \right \vert > 9s   \right) \leq 16 \exp \left\{  -\frac{N^2}{8} \right\}.
 \end{equation}
 
We will prove that $$P \left(  \hat{S}^{\nu}_{2} -  S^{\nu}_{2}  > 9s   \right)  \leq 8 \exp \left\{  -\frac{Ns^2}{8} \right\},$$  and it is analogous for the other tail.

For $1 \leq j,k \leq N$ and $\tau \in \mathcal{P}_{N,2}$ let, 
\begin{itemize}
\item $\textbf{W}_{\tau}= \left( W_{k_1},W_{k_2} \right)$  
\item $\textbf{Z}_{j}= \left( Z_j, Z_j^{\nu}  \right)$
\item $G(\textbf{Z}_{j},\textbf{W}_{\tau})= h_{\textbf{W}_{\tau}} (Z_j) h_{\textbf{W}_{\tau}} (Z_j^{\nu})$
\item $J(\textbf{Z}_{j},\textbf{W}_{\tau})= \frac{1}{2} \left[  h_{\textbf{W}_{\tau}} (Z_j)+ h_{\textbf{W}_{\tau}} (Z_j^{\nu}) \right]$
\item  $H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau})=J(\textbf{Z}_{i},\textbf{W}_{\tau})J(\textbf{Z}_{j},\textbf{W}_{\tau})$
\end{itemize}

The proof is built on three steps:

\begin{itemize}
\item \textbf{Step 1} [Rewrite the difference $\hat{S}^{\nu}_{2}-S^{\nu}_{2}$]
Similarly to  \cite{gamboa2015}  but for an $U$--statistics, we have that

\begin{align*}
\hat{S}^{\nu}_{2} = \frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau})  \\
&\hspace{-11.7cm} =\frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right\}-\\
 &\hspace{-11.7cm}  -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right \}  +  \\
&\hspace{-11.7cm} + \mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right] - \left( 1-\frac{1}{N}\right) \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{2},\textbf{W}_{\tau_1}) \right]  -   \frac{1}{N}\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right],
\end{align*}
and
\begin{align*}
S^{\nu}_{2}= \mbox{$\mathbb{E}$}_{\textbf{W}_1} \left\{ \textrm{Var}_{X_{\nu}}\left( H^{\nu}(\textbf{W}_1)\right) \right \}&\\
&\hspace{-4.5cm}  =  \mbox{$\mathbb{E}$}_{\textbf{W}_1} \left\{  \textrm{Var}_{X_{\nu}}\left( \mbox{$\mathbb{E}$}_Z (h_{\textbf{W}_1} (Z)/ X_{\nu}) \right) \right \}=  \mbox{$\mathbb{E}$}_{\textbf{W}_{1}} \left\{  \textrm{cov} \left( h_{\textbf{W}_{1}} (Z_1), h_{\textbf{W}_{1}} (Z_1^{\nu}) \right) \right \} \\
&\hspace{-4.5cm}= \mbox{$\mathbb{E}$}_{\textbf{W}_{1}} \left\{  \mbox{$\mathbb{E}$}_{\textbf{Z}_1} \left( h_{\textbf{W}_{1}} (Z_1) h_{\textbf{W}_{1}} (Z_1^{\nu}) \right) \right \}-  \mbox{$\mathbb{E}$}_{\textbf{W}_{1}}\left\{ \left[  \mbox{$\mathbb{E}$}_{\textbf{Z}_1}  \left( h_{\textbf{W}_{1}} (Z_1) \right)\right]^2 \right \}.
\end{align*}


Now, we may  decompose the error   into three terms
\begin{align*}
\hat{S}^{\nu}_{2} -  S^{\nu}_{2}= \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\substack{j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right\}}_{(A)}-\\
 &\hspace{-9cm} \underbrace{ -\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }} \left \{ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right] \right \}  }_{(B)} + \\
&\hspace{-9cm} \underbrace{+\frac{1}{N}  \left \{ \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{2},\textbf{W}_{\tau_1}) \right]  - \mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{1},\textbf{Z}_{1},\textbf{W}_{\tau_1}) \right] \right \}.}_{(C)}
\end{align*}

In the following steps, we bounded the terms $(A)$ and $(B)$ for a sufficiently large $N$. 

In  steps  $2$ and $3$ we use Hoeffding inequality (see for example \cite{hoeffding1963}) for independent and real variables. That is, if  $W_1, \ldots, W_N$ are  independent centred random variables, supported by $[a_i,b_i] \quad i=1, \ldots,N$,
\begin{equation}
P\left( \sum_{i=1}^{N} W_i > s \right) \leq  \exp \left\{ - \frac{2s^2}{\sum_{i=1}^{N}(b_i-a_i)^2} \right\} \quad (s>0).
\end{equation}

We also use an extension of the previous inequality for $U$--statistics of order $2$, see \cite{serfling1980}, Theorem A, (5.6).  If  $s=s(X_1, X_2)$ is the $U$--statistics kernel of $U_n$ such that  $\mbox{$\mathbb{E}$}(s(X_1,X_2))=\theta$ and  $a \leq s(x_1, x_2) \leq b$. So, for  $s>0$ and  $N>2$,
\begin{equation}
P\left( U_n -\theta > s \right) \leq  \exp \left\{ - \frac{Ns^2}{(b-a)^2} \right\}
\end{equation}
For us, the kernels are centred bounded indicator functions.

\item \textbf{Step 2} [Bounds for $(A)$.]

Let $G_c(\textbf{Z}_{j},\textbf{W}_{\tau})= G(\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$} \left[ G(\textbf{Z}_{j},\textbf{W}_{\tau}) \right]$ and $\tilde{G}(Z_i, \textbf{W}_{\tau})= G_c(Z_i, \textbf{W}_{\tau})  - \mbox{$\mathbb{E}$}_{Z_i} (G_c (Z_i, \textbf{W}_{\tau}))$,
\begin{displaymath}
A= \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\substack{ i \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{G}(Z_i,\textbf{W}_{\tau})}_{A_1} + \underbrace{\frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ G_c(Z,\textbf{W}_{\tau}) \right] }_{A_2}
\end{displaymath}

\begin{align*}
P(A_1>s)= P \left(\frac{1}{N\binom{N}{2}} \sum_{\substack{ i \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{G}(Z_i,\textbf{W}_{\tau}) >s \right) \\
&\hspace{-7 cm} = P \left( \sum_{i=1}^N  \frac{1}{\binom{N}{2}}\sum_{\tau \in \mathcal{P}_{N,2} } \tilde{G}(Z_i,\textbf{W}_{\tau}) >Ns\right)\\
&\hspace{-7 cm} \leq \exp \left\{  -\frac{s^2N}{2} \right\} \leq  \exp \left\{ -\frac{Ns^2}{8} \right\},
\end{align*}
while
\begin{align*}
P(A_2>s)= P \left( \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ G_c(Z,\textbf{W}_{\tau}) \right] >s \right) \\
&\hspace{-7 cm} \leq   \exp \left\{\frac{-2Ns^2}{8} \right\}\leq \exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}


But, $\left\{ A_1 + A_2> 2s \right\} \subset \{A_1>s\} \cup \{A_2>s \} $, so
\begin{equation}
P(A>2s) \leq 2\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{equation}

\item \textbf{Step 3} [Bounds for $(B)$.]

Let $H_c(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) = H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) -\mbox{$\mathbb{E}$}\left[ H(\textbf{Z}_{i},\textbf{Z}_{j},\textbf{W}_{\tau}) \right]$ and  $\tilde{H}(Z_i, Z_j , \textbf{W}_{\tau})= H_c (Z_i,Z_j, \textbf{W}_{\tau})  - \mbox{$\mathbb{E}$}_{Z_j} (H_c (Z_i, Z_j,\textbf{W}_{\tau}) ) $, 

\begin{align*}
B= \underbrace{\frac{1}{N^2\binom{N}{2}} \sum_{\substack{\{i,j\} \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} }}  \tilde{H}(Z_i,Z_j, \textbf{W}_{\tau})}_{B_1}+ \\
&\hspace{-6.5cm}  + \underbrace{\frac{1}{N\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z,\textbf{W}_{\tau}) \right] }_{B_2} + \\
&\hspace{-6.5 cm}+ \underbrace{\frac{N-1}{N^2\binom{N}{2}} \sum_{\substack{ j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} } } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] }_{B_3}.
\end{align*}

We have $\left\{ \sum_{i=1}^3B_i> 6s \right\} \subset \{B_1>3s\} \cup \{B_2>s \}\cup \{B_3>2s \} $.






So that, 
\begin{align*}
P(B_1>3s) \leq P \left(\sum_{j} \frac{1}{\binom{N}{2}}  \sum_{ \tau \in \mathcal{P}_{N,2} }  \tilde{H}(Z_j,Z_j,\textbf{W}_{\tau}) >N^2s \right)+ \\
&\hspace{-8.5 cm}+ P \left( \sum_{j>i} \frac{1}{\binom{N}{2}}  \sum_{ \tau \in \mathcal{P}_{N,2} }   \tilde{H}(Z_i,Z_j,\textbf{W}_{\tau}) >N^2s  \right)+ \\
&\hspace{-8.5 cm}+ P \left( \sum_{j<i} \frac{1}{\binom{N}{2}} \sum_{ \tau \in \mathcal{P}_{N,2} }   \tilde{H}(Z_i,Z_j,\textbf{W}_{\tau}) >N^2s \right)  \leq \\
&\hspace{-8.5cm} \leq \exp \left\{  -\frac{s^2N^3}{2} \right\} + 2 \exp \left\{  \frac{-Ns^2N^2}{4} \right\}\leq 3\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}
Further,
\begin{align*}
P(B_2>s)  \leq P \left( \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z,\textbf{W}_{\tau}) \right]  > Ns \right) \\
&\hspace{-8 cm} \leq  \exp \left\{  -\frac{Ns^2N^2}{4} \right\} \leq \exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}

Finally,
\begin{align*}
P(B_3>2s) \leq P \left( \frac{N-1}{N^2\binom{N}{2}} \sum_{\substack{ j \in \{1,\ldots, N\} \\ \tau \in \mathcal{P}_{N,2} } } \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] >2s \right)  \\
&\hspace{-10.1 cm} \leq P \left( \sum_{ j \in \{1,\ldots, N\}  } \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} }\mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] >2s \frac{N^2}{N-1}\right) \\
&\hspace{-10.1 cm}  \leq P \Bigg( \sum_{ j \in \{1,\ldots, N\}  } \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} } \Big \{ \mbox{$\mathbb{E}$}_Z \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right]- \\
&\hspace{-10.1 cm} -\mbox{$\mathbb{E}$}_{Z,Z_j} \left[ H_c(Z,Z_j,\textbf{W}_{\tau}) \right] \Big \} >\frac{s N^2}{N-1}\Bigg) + \\
&\hspace{-10.1 cm} + P \left(  \frac{1}{\binom{N}{2}} \sum_{\tau \in \mathcal{P}_{N,2} }\mbox{$\mathbb{E}$}_{Z_1,Z_2} \left[ H_c(Z_1,Z_2,\textbf{W}_{\tau}) \right]  >s \frac{N}{N-1}\right) \\
&\hspace{-10.1 cm} \leq  \exp \left\{  -\frac{s^2N^4}{2N(N-1)^2} \right\} + \exp \left\{  -\frac{2Ns^2 N^2}{8(N-1)^2} \right\} \leq 2\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{align*}


Collecting the previous inequality we obtain,
\begin{equation}
P(B>6s) \leq 6\exp \left\{  -\frac{Ns^2}{8} \right\}.
\end{equation}
\end{itemize}

\noindent  Now, using the bounds obtained in steps $2$ and $3$, we may conclude that,  for large enough  large $N$,  inequality (\ref{exponencial2}) holds.


\QEDB


\item \textbf{Proof to Example \ref{ejemplo1}.}


Let  $P(z_1, z_2)= F(z_2)-F(z_1)$, then we have 
 \begin{align*}
H(s,t)=\mathbbm{1}_{t <s} +F_Z(t) -F_Z(s))\\
 &\hspace{-3.4cm} =  \mathbbm{1}_{t <s} +(1-p)P(s,t)+pP(s-\alpha,t-\alpha),
 \end{align*}
 and,
\begin{displaymath}
H^1(s,t)=\begin{cases}
\mathbbm{1}_{t <s} + P(s,t) & \text{if } X_1=0,\\
\mathbbm{1}_{t <s} +P(s-\alpha,t-\alpha) & \text{if } X_1=1.
\end{cases}
\end{displaymath}
Therefore
\begin{small}
\begin{displaymath}
H(s,t)-H^1(s,t)=\begin{cases}
p \left[ P(s,t)-P(s-\alpha,t-\alpha)\right] & \text{if } X_1=0,\\
(1-p)\left[P(s-\alpha,t-\alpha)- P(s,t)\right]& \text{if } X_1=1,
\end{cases}
\end{displaymath}
\begin{displaymath}
\mbox{$\mathbb{E}$}_1 \left[\left(H(s,t)-H^1(s,t)\right)^2 \right] =p(1-p)\left[ P(s,t)-P(s-\alpha,t-\alpha)\right]^2,
\end{displaymath}
\end{small}
and,
\begin{align*}
S^{1}_{2}=p(1-p) \mbox{$\mathbb{E}$}\left[ \left(P(Z_1-\alpha,Z_1)-P(Z_2-\alpha,Z_2) \right)^2 \right]\\
 &\hspace{-6.9cm} = 2p(1-p) \textrm{Var} \left[ F(Z)-F(Z-\alpha) \right].
\end{align*}

Let  $X_2 \sim U(0,b)$, with $b= \sqrt{12\alpha^2 p(1-p)}$. We  obtain 



\begin{displaymath}
\mbox{$\mathbb{E}$}\left(  F(Z)- F(Z-\alpha)  \right)=\begin{cases}\lim_{N\rightarrow \infty} \frac{1}{N^2}
  \frac{\alpha}{b}   \left(1- \frac{1}{2}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/2 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

\begin{displaymath}
\mbox{$\mathbb{E}$}\left[ \left( F(Z)- F(Z-\alpha) \right)^2 \right]=\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}



\begin{displaymath}
\textrm{Var} \left( F(Z)- F(Z-\alpha) \right)=\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}

So,


\begin{displaymath}
D^{1}_{2,CVM}= p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^2  \left(1- \frac{2}{3}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/3 & \text{if } \alpha > b,
\end{cases}
\end{displaymath}

and

\begin{displaymath}
S^{1}_{2}= 2p(1-p)\begin{cases}
 \left( \frac{\alpha}{b} \right)^3  \left(\frac{1}{3}- \frac{1}{4}\frac{\alpha}{b} \right) & \text{if } \alpha \leq b,\\
 1/12 & \text{if } \alpha > b.
\end{cases}
\end{displaymath}


\end{enumerate}

\QEDB








