# Inspecting differences between multivariate distributions: graphical tool-
kit and related tests

Pratim Guha Niyogi
Johns Hopkins University
Department of Biostatistics
Baltimore, MD 21205, USA
Email ID<EMAIL_ADDRESS>Subhra Sankar Dhar
IIT Kanpur
Department of Mathematics and Statistics
Kanpur 208106, India
Email ID<EMAIL_ADDRESS>
###### Abstract

This article inspects whether a multivariate distribution is different from a
specified distribution or not, and it also tests the equality of two
multivariate distributions. In the course of this study, a graphical tool-kit
based on well-known half-spaced depth is proposed, which is a two-dimensional
plot, regardless of the dimension of the data, and it is even useful in
comparing high-dimensional distributions. The simple interpretability of the
proposed graphical tool-kit motivates us to formulate test statistics to carry
out the corresponding testing of hypothesis problems. It is established that
the proposed tests are consistent, and moreover, the asymptotic distributions
of the test statistics under contiguous alternatives are derived, which enable
us to compute the asymptotic power of these tests. Furthermore, it is observed
that the computations associated with the proposed tests are unburdensome.
Besides, these tests perform better than many other tests available in the
literature when data are generated from various distributions such as heavy
tailed distributions, which indicates that the proposed methodology is robust
as well. Finally, the usefulness of the proposed graphical tool-kit and tests
is shown on two benchmark real data sets.

Key words and phrases: goodness-of-fit test; two-sample test; half-space
depth; data-depth plot; hypothesis testing

##  1 Introduction

In many scientific investigations, data are collected for multiple variables,
and therefore, some techniques are needed to visualize and compare
multivariate observations. In this spirit, this article addresses the problem
of goodness-of-fit to validate the assumption of the underlying multivariate
distribution, which is commonly encountered in many statistical data analyses.
Moreover, we propose statistical tests to compare two multivariate probability
distributions along with appropriate graphical tool-kits. It is needless to
mention that such testing of hypothesis problems have a plethora of
applications. For example, to understand the dynamics of the human physical
activity patterns, the distribution of diurnal activity and the rest period is
assumed to be a double Pareto distribution, and in order to validate it, a
goodness-of-fit test is needed (Paraschiv-Ionescu et al., 2013). In
psychology, the study of maternal depression, examination of maternal and
infant sleep throughout the first year of life is an example of a two-sample
testing problem (Newland et al., 2016).

Let us now formulate the problem with notation. Consider the data
$\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ associated with
unknown distribution $F$, and $F_{0}$ is the specified distribution function
on $\mathbb{R}^{d}~{}(d\geq 1)$. We now want to test that $H_{0}:F=F_{0}$
against $H_{1}:F\neq F_{0}$, and a few more technical assumptions on $F$ and
$F_{0}$ will be stated at appropriate places. Moreover, as said before, the
two-sample problem for comparing multivariate distributions can also be
formulated in the following way. Suppose that
$\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ and
$\mathcal{Y}=\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\\}$ are two independent
data sets associated with two unknown multivariate distributions $F$ and $G$,
respectively. Here, we now want to test $H_{0}:F=G$ against $H_{1}:F\neq G$.
For $d=1$, i.e., in the case of univariate data, the aforesaid problems have
already been investigated in many articles. For example, the readers may refer
to Kolmogorov-Smirnov (KS) (Daniel, 1990), Cramér-von Mises (CvM) (Anderson,
1962), Anderson-Darling (AD) (Anderson and Darling, 1952) and few more tests.
In this context, for comparing univariate distributions, we would like to
mention that a few graphical tool-kits have also been used in literature for a
long time such as probability-probability plot (Michael, 1983) quantile-
quantile plot (Gnanadesikan, 2011; Wilk and Gnanadesikan, 1968; Chambers et
al., 2018), Lorenz curve (Lorenz, 1905).

For $d>1$, such statistical tests have been studied like AD test (Paulson et
al., 1987), CvM test (Koziol, 1982), Doornik-Hansen (DH) test (Doornik and
Hansen, 2008), Henze-Zirkler (HZ) test (Henze and Zirkler, 1990), Royston (R)
test (Royston, 1992) and a series of $\chi^{2}$-type tests (for example
McCulloch (M), Nikulin-Rao-Robson (NRR), Dzhaparidze-Nikulin (DN) (Voinov et
al., 2016)). These tests validate whether the data are obtained from
multivariate normality assumption or not. For two-sample multivariate
problems, similar studies have been investigated by Chen and Friedman (2017)
and see a few references therein. Besides, like univariate data, there have
been a few attempts in visualizing two multivariate distributions in two-
dimensional plots (see for example, Friedman and Rafsky (1981); Marden (1998,
2004); Dhar et al. (2014)). However, none of the above tests and graphical
tool-kits are easily implementable to large-dimensional data.

To overcome aforesaid difficulties, we introduce two test statistics based on
the difference of some functions that are equipped to measure the closeness of
an arbitrary point in a space located to an implicitly defined center of the
data (Tukey, 1975). This function is known as data-depth (for details see Liu
et al. (1999)), and the proposed test procedure is called data-depth
discrepancy (DDD). Here, DDD is large for the majority of the sample points if
and only if, the assumption on the parent distribution is likely wrong (in
goodness-of-fit test) or the samples are likely from different distributions
(in two-sample test). In this work, DDD is defined based on $L_{2}$ and
$L_{\infty}$ distance between two relevant data-depth functions which provide
us corresponding KS and CvM test statistics. However, to measure DDD, one may
use other suitable distance functions (such as $L_{p},p\geq 1$ distance) in
principle. The well-known concept of data-depth for multivariate data is
briefly reviewed next.

To study various features of multivariate distribution, there are some usual
extensions of the moment-based approach in univariate setup, and exemplary
three such features are location, scale, and shape. It is well-known that
these features can be captured by different moments and their certain
functions, if the moments exist. Even for univariate data, one can employ the
concept of the rank of the observations instead of moments to study those
features associated with descriptive statistics. However, the concept of the
rank is not uniquely defined for multivariate data unlike the univariate data.
Data-depth is a general extension of the standard univariate rank of the
distribution that is easy to plot on a plane, and therefore easy to visualize
(Liu et al., 1999). Although the rank in the univariate case assigns a ‘linear
ranking’ from the smallest to the largest data point, depth is a measure of
‘center-outward’ ranking, which starts from middle sample points and increases
outward in all directions. The highest depth is associated with the deepest or
central point, and the smallest depth corresponds to the most outlying point
with respect to the data cloud. There are various types of data-depths that
include Mahalanobis depth (Mahalanobis, 1936) and spatial depth (Gower, 1974;
Brown and Hettmansperger, 1989) based on the two moments, half-space depth
(Hodges, 1955; Tukey, 1975) based on the geometry of the half-spaces among
many others. For an overview on various depth functions, we refer to Zuo and
Serfling (2000) and the references therein.

In this article, we propose DDD based on Tukey’s half-space depth since it
uniquely determines a uniformly absolutely continuous distribution with
compact support under minimal assumptions (Koshevoy, 2002; Hassairi and
Regaieg, 2008). In addition, this can be easily computed since there exist
some packages that are freely available in various statistical softwares such
as R. In order to carry out the test, we replace the distribution function
with data-depth and propose KS and CvM type test statistics for the goodness-
fit and the two-sample tests. Along with the construction of the aforesaid
tests, following are the main contributions of this article. First, we define
a discrepancy measure based on half-space depth, which can be used to create
an alternative graphical tool-kit to visualize the disparity of two
distribution functions. Second, we propose two test statistics based on the
proposed discrepancy measure for both goodness-of-fit and two-sample testing
procedures. Third, most of the existing test statistics are available for the
normality test only; our testing procedure provides a computationally feasible
unified solution for any underlying distribution for any dimension. Fourth, we
have shown that the proposed tests are consistent. Moreover, the asymptotic
distributions of the test statistics under the contiguous alternative are
derived, which enables us to compute the asymptotic power of the tests under
contiguous alternatives. Fifth, a new graphical tool-kit based on discrepancy
measure is introduced here. The diagram related to our proposed graphical
tool-kit lies on the two-dimensional plane irrespective of the dimension of
the data. This fact enables us to visualize features of fairly large-
dimensional data using the proposed graphical tool-kit.

The rest of the paper is organized as follows. In Section 2, we briefly review
the different well-known notions of data-depth, usefulness of Tukey’s half-
space depth and a graphical display, called depth-depth plot. Section 3 is
dedicated to the proposed methods in statistical testing for goodness-of-fit
and two-sample situations. First, we provide the heuristic establishment of
the graphical tool-kit, and in later part of Section 3, we proceed with the
formal definition of DDD and conclude with two multivariate tests based on the
data-depth. In Section 4, we investigate the asymptotic properties of the
proposed testing procedures. The finite sample performance is presented in
Section 5. Finally, in Section 6, we implement the proposed test on two
interesting data sets. Section 7 concludes the article with a discussion. At
the end, Appendix contains all technical details and mathematical proofs.

The codes of all numerical studies are made available at
https://github.com/pratimguhaniyogi.

##  2 Some preliminaries

In the pioneering work by John W. Tukey in 1975, he introduced a novel way to
describe the data cloud. For a given multivariate data cloud
$\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$, a point $\mathbf{x}$
in the same Euclidean space becomes the representative of $\mathcal{X}$
through the function $D_{\mathcal{X}}(\mathbf{x})$, which measures how ‘close’
$\mathbf{x}$ is to the center of $\mathcal{X}$. Mathematically, a function,
viz., depth function, will be bounded non-negative function of the form
$D:\mathbb{R}^{d}\times\mathcal{F}\rightarrow\mathbb{R}$ where $\mathcal{F}$
is the class of all distributions on $\mathbb{R}^{d}$. In most of the examples
of the data-depth, the functions are affine invariant in the scenes
$D_{F_{\mathbf{A}\mathbf{x}+\mathbf{b}}}(\mathbf{A}\mathbf{x}+\mathbf{b})=D_{F_{\mathbf{x}}}(\mathbf{x})$
for any $d\times d$ matrix $\mathbf{A}$ and a vector
$\mathbf{b}\in\mathbb{R}^{d}$ and a quasi-concave function that maximizes
somewhere in the center of the data cloud and vanishes to infinity. Many data-
depth functions are proposed such as Mahalanobis (Mahalanobis, 1936), half-
space depth (Tukey, 1975), simplicial volume (Oja, 1983), simplicial (Liu,
1990), projection (Zuo and Serfling, 2000) etc. Mosler (2013) provides a
detailed survey of such depth measures with their properties.

Half-space depth is one of the depth functions which does not impose any
moment conditions of the data, and can be computed nonparametrically and can
characterize certain family of distributions (Koshevoy, 2003; Hassairi and
Regaieg, 2008; Cuesta-Albertos and Nieto-Reyes, 2008; Kong and Zuo, 2010). The
technical definition of the half-space depth is as follows. Let $F$ be a
probability distribution on $\mathbb{R}^{d}$ for $d\geq 1$ and $\mathcal{H}$
be the class of closed half-spaces $H$ in $\mathbb{R}^{d}$. The Tukey’s depth
(or half-space depth) of a point $\mathbf{x}\in\mathbb{R}^{d}$ with respect to
$F$ is defined by

$D_{F}(\mathbf{x})=\inf\left\\{F(H):H\in\mathcal{H},\mathbf{x}\in H\right\\}.$
(1)

One may also look at $D_{F}(\cdot)$ in the following way. Let
$\mathcal{S}^{d-1}=\\{\mathbf{u}\in\mathbb{R}^{d}:\|\mathbf{u}\|=1\\}$ be the
unit sphere of $\mathbb{R}^{d}$, then for $\mathbf{u}\in\mathcal{S}^{d-1}$ and
$\mathbf{x}\in\mathbb{R}^{d}$, we consider the closed half-space
$H[\mathbf{x},\mathbf{u}]=\\{\mathbf{y}\in\mathbb{R}^{d}:\mathbf{u}^{\text{T}}\mathbf{x}\geq\mathbf{u}^{\text{T}}\mathbf{y}\\}$.
Now, the Tukey’s half-space depth with respect to the distribution function
$F$ can be defined as

$D_{F}(\mathbf{x})=\inf_{\mathbf{u}\in\mathcal{S}^{d-1}}F(H[\mathbf{x},\mathbf{u}]).$
(2)

The sample version of $D_{F}(\mathbf{x})$, which is denoted as
$D_{\mathcal{X}}(\mathbf{x})$, based on the empirical distribution
$\widehat{F}_{n}$ of the sample
$\mathcal{X}=\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}$ and is
defined as

$D_{\mathcal{X}}(\mathbf{x})=\min_{\mathbf{u}\in\mathcal{S}^{d-1}}\frac{1}{n}\sum_{i=1}^{n}\textbf{1}\\{\mathbf{X}_{i}^{\text{T}}\mathbf{u}\leq\mathbf{x}^{\text{T}}\mathbf{u}\\}.$
(3)

Observe that, for $d=1$, $D_{F}(x)=\min\\{F(x),S(x)\\}$, where $S(x)=1-F(x)$.
Since $\widehat{F}_{n}(x)$ is the right-continuous empirical cumulative
distribution function, we define
$\widehat{S}_{n}(x)=1-\widehat{F}_{n}(x^{-})$, and consequently, the sample
version of half-space depth for univariate data becomes
$D_{\mathcal{X}}(x)=\min\\{\widehat{F}_{n}(x),\widehat{S}_{n}(x)\\}$. Based on
this definition, it can easily be observed that the depth function achieves
the maximum at the median of the distribution and monotonically decays to zero
when $x$ deviates from the median. Massé (2004) provided some results on
asymptotics for the half-space depth process, and Rousseeuw and Ruts (1996);
Ruts and Rousseeuw (1996); Rousseeuw and Struyf (1998) pointed out the
computational aspects of the half-space depth.

As a graphical tool-kit, the depth-depth (DD) plot (Liu et al., 1999) is an
important tool to compare two multivariate distributions graphically. For the
sake of completeness, we here briefly describe the methodology of the DD-plot.
We plot the depth values of the combined sample under the two corresponding
empirical distributions. Let $F$ and $G$ be two distributions on
$\mathbb{R}^{d}$, and let $D(\cdot)$ be an affine-invariant depth (for
example, half-space depth). Define the population version of DD-plot,

$DD(F,G)=\left\\{(D_{F}(\mathbf{x}),D_{G}(\mathbf{x}))\;\text{for
all}\;\mathbf{x}\in\mathbb{R}^{d}\right\\}.$ (4)

If the two distributions are identical, then these graphs are segments of the
diagonal line from (0,0) to (1,1) on the $\mathbb{R}^{2}$ plane. Plots that
deviate from the straight line indicate the difference of two underlying
distributions. For goodness-of-fit problem, if the underlying distribution $F$
is unknown, for a given sample
$\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}$, we may determine
whether $F$ is some specified distribution, say $G$, by examining the sample
version of DD-plot

$DD(\mathcal{X},G)=\left\\{(D_{F_{n}}(\mathbf{x}),D_{G}(\mathbf{x}))\;\text{for
all}\;\mathbf{x}\in\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}\right\\}.$
(5)

For two-sample problem, if $F$ and $G$ are the population distributions for
the sample $\mathcal{X}=\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}$
and $\mathcal{Y}=\left\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\right\\}$,
respectively, one can use the DD-plot to determine whether or not the two
distributions are identical.

$DD(\widehat{F}_{n},\widehat{G}_{m})=\left\\{(D_{\widehat{F}_{n}}(\mathbf{x}),D_{\widehat{G}_{m}}(\mathbf{x}))\;\text{for
all}\;\mathbf{x}\in\left\\{\mathcal{X}\cup\mathcal{Y}\right\\}\right\\}.$ (6)

It is easy to check that the DD-plot is affine-invariant if the associated
depth measures are so be. On the one hand, if $d=1$, the Lebesgue measure of
$DD$ vanishes when $F\neq G$, and on the other hand, for $d>1$ and $F$ and $G$
are absolutely continuous, $DD$ is a region with non-zero area. Moreover, Liu
et al. (1999) showed that, for different distributional differences, such as
location, scale, skewness, and kurtosis, different patterns are observed in
the DD-plot.

##  3 Data-depth discrepancy and associated tests

In this section, we propose a data-depth discrepancy (DDD) describing a
graphical tool-kit and associated tests. The DDD is defined in terms of the
difference of the depth functions, and we begin this discussion in Section 3.1
for any probability distribution. We conclude this section with the motivation
for defining a graphical tool-kit to visualize the dissimilarities of the
distributions. We introduce new multivariate goodness-of-fit and two-sample
tests based on half-space depth in Section 3.2.

###  3.1 Data-depth discrepancy (DDD)

Our goal is to construct statistical tests that can answer the following
questions.

######  Problem 1.

(One-sample/Goodness-of-fit multivariate problem) For $d\geq 1$, let
$\mathbf{X}$ be a $d$-dimensional random variable with distribution function
$F$ on $\mathbb{R}^{d}$, where $\mathbf{X}=(X_{1},\cdots,X_{d})^{\text{T}}$.
Suppose that $F_{0}$ is a pre-specified distribution function. Given
observations $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$
independent and identically distributed (i.i.d.) from $F$, can we decide
whether $F\neq F_{0}$?

######  Problem 2.

(Two-sample multivariate problem) For $d\geq 1$, let $\mathbf{X}$ and
$\mathbf{Y}$ be two random variables with distribution functions $F$ and $G$,
respectively, where $\mathbf{X}=(X_{1},\cdots,X_{d})^{\text{T}}$ and
$\mathbf{Y}=(Y_{1},\cdots,Y_{d})^{\text{T}}$. Given the two independent data
sets $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ and
$\mathcal{Y}=\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\\}$ independently and
identically distributed from $F$ and $G$, respectively, can we decide whether
$F\neq G$?

To answer the above-mentioned questions, we propose the data-depth discrepancy
(DDD) between two distributions $F$ and $G$ at $\mathbf{x}\in\mathbb{R}^{d}$
based on the data-depth as follows.

$\text{DDD}(\mathbf{x};F,G)=D_{F}(\mathbf{x})-D_{G}(\mathbf{x})$ (7)

The sample version of these measures can be obtained by replacing $D$ with its
empirical data-depth values based on the specific sample sizes. Moreover, the
discrepancy between the empirical and theoretical distribution is defined as
$\text{DDD}(\mathbf{x};\mathcal{X},F_{0})=D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})$.
Note that, $\text{DDD}(\cdot;\cdot,\cdot)$ is valid criterion to resolve
Problems 1 and 2 as long as the depth functions (here, half-space depth)
characterizes the distribution. For details, see the following propositions.

######  Proposition 1 (Koshevoy (2002); Cuesta-Albertos and Nieto-Reyes
(2008)).

If $F$ and $G$ are two distribution functions defined on $\mathbb{R}^{d}$,
both with finite support and their half-space depth coincide, then $F=G$.
Moreover, if $F$ is a discrete distribution with finite or countable support
defined on $\mathbb{R}^{d}$ and $G$ is a Borel distribution on
$\mathbb{R}^{d}$ such that the half-space depth coincides, then $F=G$.

The characterization property of the half-space depth for discrete
distributions is straightforward. For a continuous distribution function,
suppose that a hyperplane $\Delta$ divides the space $\mathbb{R}^{d}$ in two
closed half-spaces, viz, $H_{1}(\Delta)$ and $H_{2}(\Delta)$ and
$m(\Delta)=\min\\{F(H_{1}(\Delta)),F(H_{2}(\Delta))\\}$. Furthermore, assume
$B$ is the unit ball in $\mathbb{R}^{d}$. For $j=1,\cdots,d$, let
$\mathcal{S}_{j}^{d-1}$ is the relative open subspace of $\mathcal{S}^{d-1}$,
where
$\mathcal{S}_{j}^{d-1}=\left\\{\sum\limits_{i=1}^{d}u_{i}\mathbf{e}_{i}\in\mathcal{S}^{d-1}:u_{j}\neq
0\right\\}$. The map $\mathbf{u}\mapsto\boldsymbol{\beta}=-\mathbf{u}/u_{j}$
is a diffeomorphism from one of the two connected components of
$\mathcal{S}_{j}^{d-1}$ into an open subset $B_{j}$ of $\mathbb{R}^{d-1}$. Let
$\boldsymbol{\beta}\in B_{j}$, define the hyperplane
$\Delta=\Delta(\boldsymbol{\beta})$ containing a point
$\mathbf{a}=(a_{1},\cdots,a_{d})^{\text{T}}$ such that
$\Delta(\boldsymbol{\beta})=\\{\mathbf{x}\in\mathbb{R}^{d}:x_{j}=\boldsymbol{\beta}^{\text{T}}(\mathbf{x}-\mathbf{a})+a_{j}\\}$.
Thus, define two half-spaces that partitioned $\Delta(\boldsymbol{\beta})$ as
$H_{1}(\Delta(\boldsymbol{\beta}))$ and $H_{2}(\Delta(\boldsymbol{\beta}))$.
Now, suppose that $F_{j}(\boldsymbol{\beta})$ is the distribution function of
$H_{1}(\Delta(\boldsymbol{\beta}))$ defined in $\mathcal{S}_{j}^{d-1}$, and
denote the derivative of $F_{j}(\boldsymbol{\beta})$ with respect to
$\boldsymbol{\beta}$ as $F^{\prime}(\boldsymbol{\beta})$. If the hyperplane
$\Delta(\boldsymbol{\beta})$ is a minimal hyperplane at $\mathbf{a}$, i.e.,
$m(\Delta(\boldsymbol{\beta}))=D_{F}(\mathbf{a})$ then $F_{j}$ has an extremum
in $\boldsymbol{\beta}$ and $F_{j}^{\prime}(\boldsymbol{\beta})=0$. Therefore,
under the condition that $F_{j}^{\prime}(\boldsymbol{\beta})=0$ implies
$\Delta(\boldsymbol{\beta})$ is a minimal hyperplane at $\mathbf{a}$, the
following result holds.

######  Proposition 2 (Hassairi and Regaieg (2008)).

Let $F$ be a probability distribution on $\mathbb{R}^{d}$ with connected
support, and $\mathcal{H}$ be the class of closed half-spaces $H$ in
$\mathbb{R}^{d}$. Assume that $F$ is absolutely continuous with the connected
support and that the corresponding probability density function is continuous
with the interior of the support. Then the half-space depth characterizes the
distribution.

Propositions 1 and 2 indicate that the difference of the half-space depth
values can provide a valid measure of the discrepancy, which is defined as DDD
in Equation (7). Now we propose a graphical method where we plot DDD as a
scatter plot across a horizontal axis and measure the deviation from the
horizontal axis. Needless to say, the advantage of this graphical tool-kit is
that it remains unaffected by the data dimension. This graphical approach can
be used for any distribution, but to produce a valid graphical tool-kit for
addressing the Problems 1 and 2, we continue our discussion with half-space
depth since we can exploit the characterization properties discussed in
Propositions 1 and 2 of half-space depth. As a result, we can conclude that,
if $\text{DDD}(\mathbf{x};F,G)=0$ for all $\mathbf{x}$, and graphically, if
the majority of points are concentrated on the horizontal axis, we can
conclude that the two underlying distributions are expected to be identical,
where $F$ and $G$ belong to a certain family of distribution functions. The
illustrative examples provide the different patterns of deviations from the
horizontal line for both goodness-of-fit and two-sample problems.

####  3.1.1 Illustration of graphical tool-kit: Goodness-of-fit testing
problem

For a given sample
$\mathcal{X}=\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}$ from an
unknown distribution function $F$, we are interested in plotting
$\text{DDD}(\mathbf{x};\mathcal{X},F_{0})=D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})$
with respect to indices of $\mathbf{x}$. Therefore, if two distributions are
identical, then the value of $\text{DDD}(\mathbf{x};\mathcal{X},F_{0})$ will
be a straight line with respect to the observed data. Thus, this method
provides some idea of goodness-of-fit. We consider four simulated data sets
each consisting of 100 i.i.d. observations. Four sets of samples are generated
from (1) standard bivariate normal distribution, (2) standard bivariate Cauchy
distribution, (3) bivariate t-distribution with degrees of freedom (df) 3 and
(4) standard bivariate Laplace distribution, respectively. For all of them, we
consider the bivariate normal distribution as the specified distribution
$F_{0}$ with unknown mean $\boldsymbol{\mu}$ and dispersion
$\boldsymbol{\Sigma}$, that are estimated from the sample using the sample
mean vector and the dispersion matrix, receptively. We standardize the sample
using these estimates and see the DDD-plots based on standardized data. The
data-depth discrepancy plots for the four simulated samples are shown in
Figure 1. In those plots, the dotted gray curves in the figures indicate the
two-sigma limits of DDD. It is clearly evident from the graphs that the
specified distribution fits the data in Figure 1(a) very well; moreover, the
data points in this plot are clustered around the straight line and most of
them are in the two-sigma limits. The other three plots 1(b), 1(c), 1(d)
indicate deviations from the straight line, which immediately implies that the
data are not from a normal distribution.

(a) $\mathbf{X}\sim\text{standard normal}$

(b) $\mathbf{X}\sim\text{standard Cauchy}$

(c) $\mathbf{X}\sim\text{standard Laplace}$

(d) $\mathbf{X}\sim\text{standard t with df 3}$

Figure 1: The data-depth discrepancy plot for the goodness-of-fit problem
where the specified distribution is bivariate normal. The plots of the first
and second rows are for the examples, where the distribution of the data are
form bivariate standard normal, Cauchy, Laplace and $t_{3}$ distributions
respectively. The dotted gray curves indicate the two-sigma limits of data-
depth discrepancy. Points which are outside the two-sigma limits are color
coded by red.

####  3.1.2 Illustration of graphical tool-kit: Two-sample testing problem

For given samples $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ and
$\mathcal{Y}=\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\\}$, here we plot
$\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})=D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x})$
for $\mathbf{x}\in\mathcal{X}\cup\mathcal{Y}$. If the two distributions are
identical, then the value of $\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})$
will be in a straight line with respect to observed data. Here we consider two
simulated data sets to demonstrate our proposed graphical tool-kit. In our
first problem, we simulate sample-1 consisting of 100 i.i.d. observations
generated from the trivariate normal distribution having zero mean and scatter
matrix $\boldsymbol{\Sigma}=(\sigma_{i,j})_{i,j=1,2,3,4}$ with
$\sigma_{1,2}=0.9,\sigma_{1,3}=0.2$ and $\sigma_{2,3}=0.5.$ $(F)$ and sample-2
consisting of 50 i.i.d. observations from the standard trivariate normal
distribution denoted as $(G)$. In the next problem, sample-1 consists of 100
i.i.d. observations from the standard trivariate normal distribution $(F)$,
and sample-2 consists of 50 i.i.d. observations from a trivariate skew-normal
distribution $(G)$ (Azzalini and Valle, 1996). The p.d.f. of the trivariate
skew-normal distribution is given by
$f(\mathbf{x})=2\phi_{3}(\mathbf{x};\boldsymbol{\Omega})\Phi(\boldsymbol{\alpha}^{\text{T}}\mathbf{x})$,
where,
$\boldsymbol{\alpha}^{\text{T}}=\frac{\boldsymbol{\lambda}^{\text{T}}\boldsymbol{\Psi}^{-1}\boldsymbol{\Delta}^{-1}}{\sqrt{1+\boldsymbol{\lambda}^{\text{T}}\boldsymbol{\Psi}^{-1}\boldsymbol{\lambda}}}$,
$\Delta=diag(\sqrt{1-\lambda_{1}^{2}},\sqrt{1-\lambda_{2}^{2}},\sqrt{1-\lambda_{3}^{2}})$,
$\boldsymbol{\lambda}=\left(\frac{\lambda_{1}}{\sqrt{1-\lambda_{1}^{2}}},\frac{\lambda_{2}}{\sqrt{1-\lambda_{2}^{2}}},\frac{\lambda_{3}}{\sqrt{1-\lambda_{3}^{2}}}\right)^{\text{T}}$,
and
$\boldsymbol{\Omega}=\boldsymbol{\Delta}(\boldsymbol{\Psi}+\boldsymbol{\lambda}\boldsymbol{\lambda}^{\text{T}})\boldsymbol{\Delta}$.
Here $\phi_{3}(\mathbf{x};\boldsymbol{\Omega})$ denotes the probability
density function of a trivariate normal distribution with standardized
marginals and the correlation matrix matrix $\boldsymbol{\Omega}$, and
$\boldsymbol{\Psi}$ is the distribution function of the standard univariate
normal distribution. In this example, we have considered
$\lambda_{1}=\lambda_{2}=\lambda_{3}=0.9$ and
$\boldsymbol{\Phi}=\mathbf{I}_{3}$, identity matrix of dimension $3\times 3$.
The data-depth discrepancy plots for these two toy examples are shown in
Figure LABEL:fig:two-sample. As before, the gray dotted curves in the figures
indicate the two-sigma limits of DDD.

####  3.1.3 Illustration of graphical tool-kit: High-dimension data

In the last decades or so, there have been plenty of applications involved
high-dimensional data, and we here illustrate the usefulness of proposed
graphical tool-kit on to meet the challenges of the modern data applications.
By varying the dimensions of the data sets, $d=10,20,30,40,50$ and $60$, we
generate sample-1 from $d$-variate normal distribution and sample-2 from
$d$-variate Cauchy distribution with sample size $100$ each for a two-sample
testing problem. The data-depth discrepancy plots for different dimensions are
shown in Figure 2 which illustrates the fact that differences between the
depth values of two samples are not clustered around the the horizontal axis.
This indicates that the proposed graphical tool-kit is infact useful for high-
dimensional data.

(a) $d=10$

(b) $d=20$

(c) $d=30$

(d) $d=40$

(e) $d=50$

(f) $d=60$

Figure 2: The data-depth discrepancy plot for the two-sample examples in high-
dimensional situation for different dimensions. Here sample are from
multivariate normal and Cauchy distributions. The dotted gray curves indicate
the two-sigma limits of data-depth discrepancy. Points which are outside the
two-sigma limits are color coded by red.

###  3.2 Associated tests

In Problem 1, for given sample $\mathcal{X}$ of size $n$, we are interested to
test $H_{0}:F=F_{0}$ against $H_{1}:F\neq F_{0}$, where $F_{0}$ is a pre-
specified distribution function. Due to Propositions 1 and 2, for half-space
depth function $D$, it is equivalent to test
$H^{*}_{0}:D_{F}(\mathbf{x})=D_{F_{0}}(\mathbf{x})\;\text{for
all}\;\mathbf{x}\in\mathbb{R}^{d}$ against $H^{*}_{1}:D_{F}(\mathbf{x})\neq
D_{F_{0}}(\mathbf{x})$ for some $\mathbf{x}\in\mathbb{R}^{d}$. The most well-
known test statistics for comparing two distributions are the
Kolmogorov–Smirnov (KS) test statistics and Cramér–von Mises (CvM) test
statistics. These test statistics are based on the certain differences between
the empirical distribution function $F_{n}$ and the hypothesised distribution
function $F_{0}$. Here, we propose alternate test statistics based on these
commonly known families of test statistics and DDD defined in Section 3.1 so
that we can perform the goodness-of-fit test for any arbitrary dimension of
the data.

  1. 1.

Depth-based KS test statistic:

$T^{\text{KS}}_{\mathcal{X},F_{0}}=\sqrt{n}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|\text{DDD}(\mathbf{x};\mathcal{X},F)|=\sqrt{n}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|$
(8)

  2. 2.

Depth-based CvM test statistic:

$T^{\text{CvM}}_{\mathcal{X},F_{0}}=n\int\text{DDD}^{2}(\mathbf{x};\mathcal{X},F_{0})dF_{0}(\mathbf{x})=n\int(D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x}))^{2}dF_{0}(\mathbf{x}),$
(9)

where the integral is over a closed ball in $\mathbb{R}^{d}$ with the center
at the origin.

For testing $H_{0}$ against $H_{1}$, the null hypothesis will be rejected when
the values of the test statistics are very large. The asymptotic behaviors of
the proposed test statistics are provided in Section 4. To obtain the
empirical p-value based on the proposed test statistics
$T_{\mathcal{X},F_{0}}$, where $T_{\mathcal{X},F_{0}}$ is either
$T_{\mathcal{X},F_{0}}^{\text{KS}}$ or $T_{\mathcal{X},F_{0}}^{\text{CvM}}$,
we consider a parametric bootstrap procedure (Efron, 1982) that consists of
the following steps.

  1. Step 1.1.

Generate $M$ points uniformly from the $d$-dimensional unit sphere, for $M$
sufficiently large. The set of such points is denoted as $\mathcal{U}$. The
algorithm for generating a single sample from a $d$-dimensional sphere with
radius $r_{0}$ is as follows.

    1. (a)

Generate $U$ from uniform distribution over the set $(0,r_{0})$ and define
$R=\sqrt{U}$.

    2. (b)

$\theta_{1},\cdots,\theta_{d-2}\sim Uniform(0,\pi)$, $\theta_{d-1}\sim
Uniform(0,2\pi)$ and each of $\theta_{j}$s are mutually independent for
$j=1,\cdots,d-1$.

    3. (c)

A $d$-dimensional observation from the sphere is
$\mathbf{u}=(u_{1},\cdots,u_{d})^{\text{T}}$ where,

$\displaystyle u_{1}=R\cos\theta_{1}$ $\displaystyle
u_{2}=R\sin\theta_{1}\cos\theta_{2}$ $\displaystyle\vdots$ $\displaystyle
u_{d-1}=R\sin\theta_{1}\sin\theta_{2}\cdots\sin\theta_{d-2}\cos\theta_{d-1}$
$\displaystyle
u_{d}=R\sin\theta_{1}\sin\theta_{2}\cdots\sin\theta_{d-2}\sin\theta_{d-1}.$

  2. Step 1.2.

Calculate half-space depth with respect to samples $\mathcal{X}$ and $F_{0}$,
respectively, and therefore, calculate
$\text{DDD}(\mathbf{x};\mathcal{X},F_{0})$ for all $\mathbf{x}\in\mathcal{U}$.
Construct an approximate test statistic $\widetilde{T}_{\mathcal{X},F_{0}}$
which is either
$\widetilde{T}^{\text{KS}}_{\mathcal{X},F_{0}}=\sqrt{n}\max_{\mathbf{x}\in\mathcal{U}}|\text{DDD}(\mathbf{x};\mathcal{X},F_{0})|$
or
$\widetilde{T}^{\text{CvM}}_{\mathcal{X},F_{0}}=n\times\frac{1}{M}\sum\limits_{j=1}^{M}\text{DDD}^{2}(\mathbf{x}_{j};\mathcal{X},F_{0})$
based on the choice of $T_{\mathcal{X},F_{0}}$.

  3. Step 1.3.

Let $\theta$ be the unspecified parameter of the distribution function $F_{0}$
(denoted as $F_{0}(\cdot;\theta)$) and $\widehat{\theta}_{n}$ be an estimate
of $\theta$ which could be the maximum likelihood estimate.

  4. Step 1.4.

Draw a random sample
$\mathcal{X}^{*}=\\{\mathbf{X}^{*}_{1},\cdots,\mathbf{X}^{*}_{n}\\}$ from
fitted distribution $F_{0}(\cdot;\widehat{\theta}_{n})$.

  5. Step 1.5.

Compute the half-space depth with respect to $\mathcal{X}^{*}$ and, therefore,
compute
$\text{DDD}(\mathbf{x};\mathcal{X}^{*},F_{0}(\cdot;\widehat{\theta}_{n}))=D_{\mathcal{X}^{*}}(\mathbf{x})-D_{F_{0}(\cdot;\widehat{\theta}_{n})}(\mathbf{x})$
for $\mathbf{x}\in\mathcal{U}$. Thus, calculate the test statistic
$\widetilde{T}_{\mathcal{X}^{*},F_{0}(\cdot;\widehat{\theta}_{n})}$ based on
$\text{DDD}(\mathbf{x};\mathcal{X}^{*},F_{0}(\cdot;\widehat{\theta}_{n}))$.

  6. Step 1.6.

Repeat Steps 1.4 and 1.5 for $B$ times, where $B$ is sufficiently large. Thus,
the empirical distribution of
$\\{\widetilde{T}^{(b)}_{\mathcal{X}^{*},F_{0}(\cdot;\widehat{\theta}_{n})}:b=1,\cdots,B\\}$
can be used to approximate the null distribution of ${T}_{\mathcal{X},F_{0}}$.

  7. Step 1.7.

The empirical p-value is computed using the following formula.

$\widehat{\mathbb{P}}_{H_{0}}\left\\{T_{\mathcal{X},F_{0}}>\widetilde{T}_{\mathcal{X},F_{0}}\right\\}=\frac{1}{B}\sum_{b=1}^{B}\textbf{1}\left\\{\widetilde{T}_{\mathcal{X}^{*},F_{0}(\cdot;\widehat{\theta}_{n})}^{(b)}>\widetilde{T}_{\mathcal{X},F_{0}}\right\\}$
(10)

Next, we consider the two-sample $d$-dimensional problem for two independent
samples $\mathcal{X}=\left\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\right\\}$
and $\mathcal{Y}=\left\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\right\\}$ where
$\mathbf{X}_{i}$’s are independently distributed with distribution function
$F$ and $\mathbf{Y}_{i}$’s are independently distributed with distribution
function $G$. In Problem 2, we want to test $H_{0}:F=G$ against $H_{1}:F\neq
G$, which is equivalent to test
$H^{*}_{0}:D_{F}(\mathbf{x})=D_{G}(\mathbf{x})\;\text{for
all}\;\mathbf{x}\in\mathbb{R}^{d}$ against $H^{*}_{1}:D_{F}(\mathbf{x})\neq
D_{G}(\mathbf{x})$ for some $\mathbf{x}\in\mathbb{R}^{d}$, where $D$ is the
half-space depth. Similar to goodness-of-fit test problem, the above testing
procedure is meaningful due to the characterization property of half-space
depth. We construct the KS and CvM type test statistics based on
$\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})$ and obtain the following test
statistics.

  1. 1.

Depth-based KS test statistic:

$T^{\text{KS}}_{\mathcal{X},\mathcal{Y}}=\sqrt{n+m}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})|=\sqrt{n+m}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x})|$
(11)

  2. 2.

Depth-based CvM test statistic:

$T^{\text{CvM}}_{\mathcal{X},\mathcal{Y}}=(n+m)\int\text{DDD}^{2}(\mathbf{x};\mathcal{X},\mathcal{Y})dH_{n,m}(\mathbf{x})=(n+m)\int(D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x}))^{2}dH_{n,m}(\mathbf{x}),$
(12)

where the integral is over a closed ball in $\mathbb{R}^{d}$ with the center
at origin, and $H_{n,m}(\cdot)$ is the empirical distribution function based
on the combine sample $\mathcal{X}\cup\mathcal{Y}$.

Here also, for large value of the test statistics, we reject the null
hypothesis. The asymptotic results for the proposed test statistics are
studied in Section 4. The following algorithm provides an empirical p-value
(Tibshirani and Efron, 1993) for Problem 2 based on any of the proposed test
statistics $T_{\mathcal{X},\mathcal{Y}}$ where $T_{\mathcal{X},\mathcal{Y}}$
is either $T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}$ or
$T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}$.

  1. Step 2.1.

Generate $M$ points uniformly from the $d$-dimensional unit sphere, for
sufficiently large $M$. The set of such points is denoted as $\mathcal{U}$.

  2. Step 2.2.

Calculate half-space depth with respect to samples $\mathcal{X}$ and
$\mathcal{Y}$ respectively, and therefore compute
$\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})$ for all
$\mathbf{x}\in\mathcal{U}$. Construct an approximate test statistic
$\widetilde{T}_{\mathcal{X},\mathcal{Y}}$ which is either
$\widetilde{T}^{\text{KS}}_{\mathcal{X},\mathcal{Y}}=\sqrt{n+m}\max_{\mathbf{x}\in\mathcal{U}}|\text{DDD}(\mathbf{x};\mathcal{X},\mathcal{Y})|$
or
$\widetilde{T}^{\text{CvM}}_{\mathcal{X},\mathcal{Y}}=(n+m)\times\frac{1}{M}\sum\limits_{j=1}^{M}\text{DDD}^{2}(\mathbf{x}_{j};\mathcal{X},\mathcal{Y})$.

  3. Step 2.3.

Combine the samples $\mathcal{X}$ and $\mathcal{Y}$, and the pooled sample is
denoted as
$\mathcal{Z}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}\cup\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\\}$.

  4. Step 2.4.

For the $b$-th bootstrap replicate, take a sample of size $(n+m)$ from
$\mathcal{Z}$ with replacement and treat the first $n$ as the
$\mathcal{X}_{b}^{*}$ sample and the last $m$ are $\mathcal{Y}_{b}^{*}$ for
$b=1,\cdots,B$, where $B$ is sufficiently large.

  5. Step 2.5.

Compute half-space depths based on the samples $\mathcal{X}_{b}^{*}$ and
$\mathcal{Y}_{b}^{*}$ separately for each of the replicates. Therefore,
compute the approximate test statistic
$\widetilde{T}_{\mathcal{X}^{*},\mathcal{Y}^{*}}^{(b)}$ as described in Step
2.2 by replacing $\mathcal{X}$ with $\mathcal{X}_{b}^{*}$ and $\mathcal{Y}$
with $\mathcal{Y}_{b}^{*}$. Thus, the empirical distribution of
$\widetilde{T}^{(b)}_{\mathcal{X}^{*},\mathcal{Y}^{*}}$s can be used to
approximate the null distribution of ${T}_{\mathcal{X},\mathcal{Y}}$.

  6. Step 2.6.

The empirical p-value is computed using the following formula.

$\widehat{\mathbb{P}}_{H_{0}}\left\\{T_{\mathcal{X},\mathcal{Y}}>\widetilde{T}_{\mathcal{X},\mathcal{Y}}\right\\}=\frac{1}{B}\sum_{b=1}^{B}\textbf{1}\left\\{\widetilde{T}_{\mathcal{X}^{*},\mathcal{Y}^{*}}^{(b)}>\widetilde{T}_{\mathcal{X},\mathcal{Y}}\right\\}$
(13)

######  Remark 1.

Using the Algorithms described above, it is observed that the computation of
the p-values is simple and efficient enough even when the dimension of the
data is fairly large.

##  4 Main results

In the earlier section, we studied the algorithms to compute the p-value of
the tests. However, one needs to know the distribution of the test statistics
to estimate the size and power of the tests. Observe that the derivation of
the exact distribution of the test statistics is intractable because of its
complex expression, and to overcome this issue, this section investigates the
asymptotic consistency and the asymptotic power under contiguous alternatives
of the proposed test using the asymptotic distribution of the test statistics.

###  4.1 Large sample properties

To investigate the asymptotic properties of the proposed graphical tool-kit
and tests, one first needs to assume a few technical conditions.

Suppose that $\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ are the multivariate
observations from an unknown distribution $F$, and the corresponding empirical
distribution function is denoted by $\widehat{F}_{n}$. We define
$\mathcal{L}_{n}(\mathbf{x})=\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))$
for $\mathbf{x}\in\mathbb{R}^{d}$. Therefore,
$\\{\mathcal{L}_{n}:\mathbf{x}\in\mathbb{R}^{d}\\}$ is a stochastic process
with bounded sample path which is a map into $l^{\infty}(\mathbb{R}^{d})$, a
space of bounded real valued function on $\mathbb{R}^{d}$ equipped with the
uniform norm, viz.,
$\|a\|_{\infty}=\sup\limits_{\mathbf{t}\in\mathbb{R}^{d}}|a(\mathbf{t})|$.
Moreover, define $\mathcal{V}_{n}=\sqrt{n}(\widehat{F}_{n}-F)$, and that can
be viewed as a map into space $l^{\infty}(\mathcal{H})$, where $\mathcal{H}$
is the class of closed half-spaces $H$ in $\mathbb{R}^{d}$. The technical
assumptions:

  1. (C1)

Let $\partial H$ be the topological boundary of $H$. Then the distribution
function $F$ satisfies $F(\partial H)=0$ for all $H\in\mathcal{H}$.

  2. (C2)

Let $\mathcal{F}$ be a point-wise bounded, totally bounded and permissible
subset of $L_{2}(F)$ on $\mathcal{H}$, where $L_{2}(F)$ is the set of all
square-integrable functions with respect to $F$. For each $\eta>0$ and
$\epsilon>0$, there exits a $\delta>0$ such that,
$\lim\sup\mathbb{P}\left\\{\sup_{\mathcal{D}(\delta)}|\mathcal{V}_{n}(f-g)|>\eta\right\\}<\epsilon$,
where $\mathcal{D}(\delta)=\\{(f,g):f,g\in\mathcal{F},\rho_{F}(f-g)<\delta\\}$
for the semi-norm $\rho_{F}$ on $\mathcal{H}$.

  3. (C3)

At any point with positive depth, the collection of all minimal directions
passing through that point has a finite number of elements or equals with
$\mathcal{S}^{d-1}=\\{\mathbf{u}:\|\mathbf{u}\|=1\\}$.

  4. (C4)

A minimal closed half-space at $\mathbf{x}$, viz. $H[\mathbf{x}]$ is uniquely
defined if $D_{F}(\mathbf{x})>0$.

  5. (C5)

The distribution functions are either with finite support or absolutely
continuous with continuous probability density function.

######  Remark 2.

(C1) is the smoothness condition of the underlying distribution function. This
is useful to obtain the stronger properties of the Tukey’s half-space depth in
different applications. For example, it is easy to see that if $F$ is
absolutely continuous, the condition (C1) is satisfied (Massé, 2004).
Moreover, if $F$ satisfies (C1), $D_{F}(\mathbf{x})$ can be expressed as the
probability of some closed half-space whose boundary goes through
$\mathbf{x}$. In addition to that, for the closed half-space
$H[\mathbf{x},\mathbf{u}]=\\{\mathbf{y}\in\mathbb{R}^{d}:\mathbf{u}^{\text{T}}\mathbf{x}\geq\mathbf{u}^{\text{T}}\mathbf{y}\\}$,
the function $(\mathbf{x},\mathbf{u})\mapsto F(H[\mathbf{x},\mathbf{u}])$ is
continuous on $\mathbb{R}^{d}\times\mathcal{S}^{d-1}$ and
$\mathbf{x}\rightarrow D_{F}(\mathbf{x})$ is continuous. Condition (C2) are
useful for empirical central limit theorems (see Pollard (1984). If
$D(\mathbf{x})=F(H[\mathbf{x},\mathbf{u}])$, $H[\mathbf{x},\mathbf{u}]$ is a
minimal half-space at $\mathbf{x}$ and $\mathbf{u}$ is a minimal direction at
the same point. (C3) satisfies the condition related to the multiplicity of
the minimal direction at $\mathbf{x}$. The conditions (C1) and (C3) are
referred as the “local regularity conditions” in the literature. The details
about the condition (C5) is discussed in Propositions 1 and 2.

Proposition 3 states the weak convergence of $\mathcal{V}_{n}$.

######  Proposition 3.

(Theorem 21 in section VII.5 of Pollard (1984)) Under condition (C2),
$\mathcal{V}_{n}=\sqrt{n}(\widehat{F}_{n}-F)$ converges weakly to
$\mathcal{V}_{F}$ where $\mathcal{V}_{F}$ is tight and $F$-Brownian bridge
with mean zero and covariance function $\Sigma_{F}=F(H_{1}\cap
H_{2})-F(H_{1})F(H_{2})$ for $H_{1},H_{2}\in\mathcal{H}$. Moreover,
$\mathcal{V}_{F}$ can be chosen such that each sample path is continuous with
respect to $\rho$.

Let us define
$\mathcal{J}(\mathcal{V}_{F})(\mathbf{x})=\inf\limits_{\mathbf{v}\in
V(\mathbf{x})}\mathcal{V}_{F}H[\mathbf{x},\mathbf{v}]$ for
$\mathbf{x}\in\mathbb{R}^{d}$, where $V(\mathbf{x})$ is the set of all minimal
directions passing thought that point. Proposition 4 states the asymptotic
distribution of the Tukey’s half-space depth.

######  Proposition 4.

(Massé, 2004) Under the conditions (C1)-(C4), for $A$ is a closed subset of
$\mathcal{S}_{\mathcal{F}}$, in $l^{\infty}(A)$, $\\{\mathcal{L}_{n}\\}$
converges weakly to $\mathcal{J}(\mathcal{V}_{F})$, which is tight measurable
map into $l^{\infty}(A)$. Moreover, $\mathcal{J}(\mathcal{V}_{F})$ is a random
element associated with a centered Gaussian process with covariance function
$\text{cov}\\{\mathcal{J}(\mathcal{V}_{F})(\mathbf{x}_{1}),\mathcal{J}(\mathcal{V}_{F})(\mathbf{x}_{2})\\}=F(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F(H[\mathbf{x}_{1}])F(H[\mathbf{x}_{2}])$ for
$\mathbf{x}_{1},\mathbf{x}_{2}\in A$.

We now state the results, which justifies the usefulness of the proposed
graphical tool-kits for sufficiently large sample sizes. See Theorems 1 and 2:

######  Theorem 1.

For every $\epsilon>0$, define
$\mathcal{C}_{\epsilon}(F,F_{0})=\\{(\mathbf{x},D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x})):\mathbf{x}\in\mathbb{R}^{d},|D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|<\epsilon\\}$
and for a given sample
$\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$,
$\widehat{\mathcal{C}}(\mathcal{X},F_{0})=\\{(\mathbf{x},D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})):\mathbf{x}\in\mathcal{X}\\}$.
Then, under the conditions (C1)-(C5),

$\lim_{n\rightarrow\infty}\mathbb{P}\left\\{\widehat{\mathcal{C}}(\mathcal{X},F_{0})\subset\mathcal{C}_{\epsilon}(F,F_{0})\right\\}=1.$
(14)

######  Theorem 2.

For every $\epsilon>0$, define
$\mathcal{C}_{\epsilon}(F,G)=\\{(\mathbf{x},D_{F}(\mathbf{x})-D_{G}(\mathbf{x})):\mathbf{x}\in\mathbb{R}^{d},|D_{F}(\mathbf{x})-D_{G}(\mathbf{x})|<\epsilon\\}$
and for given independent samples
$\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ and
$\mathcal{Y}=\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{m}\\}$,
$\widehat{\mathcal{C}}(\mathcal{X},\mathcal{Y})=\\{(\mathbf{x},D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x})):\mathbf{x}\in\mathcal{X}\cup\mathcal{Y}\\}$.
Then, under the conditions (C1)-(C5), for positive finite number
$\lambda=\lim\limits_{\min(n,m)\rightarrow\infty}n/(n+m)$,

$\lim_{\min(n,m)\rightarrow\infty}\mathbb{P}\left\\{\widehat{\mathcal{C}}(\mathcal{X},\mathcal{Y})\subset\mathcal{C}_{\epsilon}(F,G)\right\\}=1.$
(15)

######  Remark 3.

Theorems 1 and 2 show that for large sample size the points in the sets lie
around the horizontal axis if and only if the null hypothesis is expected to
be true for both the goodness-of-fit and two-sample testing problem under some
mild conditions.

Next, Theorems 3 and 4 assert the point-wise asymptotic properties of the
goodness-of-fit tests based on $T_{\mathcal{X},F_{0}}^{\text{KS}}$ and
$T_{\mathcal{X},F_{0}}^{\text{CvM}}$.

######  Theorem 3.

Under the conditions (C1)-(C5), the test based on test statistic
$T_{\mathcal{X},F_{0}}^{\text{KS}}$ for testing $H_{0}:F=F_{0}$ against
$H_{1}:F\neq F_{0}$ is point-wise consistent in power, i.e.,
$\mathbb{P}_{F}\\{T_{\mathcal{X},F_{0}}^{KM}>s_{1-\alpha}^{(1)}\\}\rightarrow
1$ as $n\rightarrow\infty$ under $H_{1}$, where $s_{1-\alpha}^{(1)}$ is such
that
$\lim\limits_{n\rightarrow\infty}\mathbb{P}_{H_{0}}\left\\{T_{\mathcal{X},F_{0}}^{\text{KS}}>s_{1-\alpha}^{(1)}\right\\}=\alpha$.

######  Theorem 4.

Under the conditions (C1)-(C5), the test based on test statistic
$T_{\mathcal{X},F_{0}}^{\text{CvM}}$ for testing $H_{0}:F=F_{0}$ against
$H_{1}:F\neq F_{0}$ is point-wise consistent in power, i.e.,
$\mathbb{P}\left\\{T_{\mathcal{X},F_{0}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}\rightarrow
1$ as $n\rightarrow\infty$, where $s_{1-\alpha}^{(2)}$ is such that
$\lim\limits_{n\rightarrow\infty}\mathbb{P}_{H_{0}}\left\\{T_{\mathcal{X},F_{0}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}=\alpha$.

Now, Theorems 5 and 6 assert the asymptotically uniform power properties of
the goodness-of-fit test based on $T_{\mathcal{X},F_{0}}^{\text{KS}}$ and
$T_{\mathcal{X},F_{0}}^{\text{CvM}}$.

######  Theorem 5.

Suppose that $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ is a
collection of i.i.d. random variables with distribution function $F$. For
testing $F=F_{0}$ against $H_{1}:F\neq F_{0}$, under the conditions (C1)-(C5),
the power of the test based on $T_{\mathcal{X},F_{0}}^{\text{KS}}$ tends to 1
uniformly over all alternatives $F$ satisfying
$\sqrt{n}d_{K}(D_{F},D_{F_{0}})\geq\Delta_{n}$, if
$\Delta_{n}\rightarrow\infty$ as $n\rightarrow\infty$.

######  Theorem 6.

Suppose that $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ is a
collection of i.i.d. random variables with distribution function $F$. For
testing $H_{0}:F=F_{0}$ against $H_{1}:F\neq F_{0}$, under the conditions
(C1)-(C5), the power of test based on $T_{\mathcal{X},F_{0}}^{\text{CvM}}$
tends to 1 uniformly over all alternatives $F_{n}$ satisfying
$\sqrt{n}d_{K}(D_{F_{n}},D_{F_{0}})\geq\Delta_{n}$ if
$\Delta_{n}\rightarrow\infty$ as $n\rightarrow\infty$.

We now describe the results related to point-wise and uniform consistency of
the proposed two-sample testing procedure.

######  Theorem 7.

Let us denote
$T^{(1)}_{\mathcal{X},\mathcal{Y}}:=T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}$
and
$T^{(2)}_{\mathcal{X},\mathcal{Y}}:=T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}$,
and $n$ and $m$ are such that
$\lim\limits_{\min(n,m)\rightarrow\infty}\frac{n}{n+m}=\lambda\in(0,1)$. For
$i=1$ and $2$, let $t_{1-\alpha}^{(i)}$ be such that
$\lim\limits_{\min(n,m)\rightarrow\infty}\mathbb{P}_{H_{0}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{(i)}>t_{1-\alpha}^{(i)}\right\\}=\alpha$.
Further, suppose that $H(\mathbf{x})=\lambda
F(\mathbf{x})+(1-\lambda)G(\mathbf{x})$ and
$\int(D_{F}(\mathbf{x})-D_{G}(\mathbf{x}))dH(\mathbf{x})<\infty$. Then, under
the conditions (C1)-(C5),
$\mathbb{P}_{H_{1}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{(i)}>t_{1-\alpha}^{(i)}\right\\}\rightarrow
1$, as $\min(n,m)\rightarrow\infty$. Moreover, for $i=1$ and $2$, the power of
the test based on $T^{(i)}_{\mathcal{X},\mathcal{Y}}$ tends to one, uniformly
over all alternatives satisfying
$\sqrt{n+m}d_{k}(D_{F_{n}},D_{G_{m}})\geq\Delta_{m,n}$ if
$\Delta_{m,n}\rightarrow\infty$ as $\min(m,n)\rightarrow\infty$.

###  4.2 Asymptotic local power study

In Section 4.1, we have established that data-depth based KS and CvM tests are
all asymptotically consistent, and therefore, a natural question is how is the
asymptotic power of the tests under local/contiguous alternatives (see Sidak
et al. (1999)). Let $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$ be the sequences of
the probability measures defined on sequence of probability spaces
$(\Omega_{n},\mathcal{A}_{n})$. Then, $\mathbb{Q}_{n}$ is said to be
contiguous with respect to $\mathbb{P}_{n}$ when
$\mathbb{P}_{n}\\{A_{n}\\}\rightarrow 0$ implies that
$\mathbb{Q}_{n}\\{A_{n}\\}\rightarrow 0$ for every sequence of measurable sets
$A_{n}$. It is important to note that the sequence of set $A_{n}$ changes with
$n$ along with the $\sigma-$field $\mathcal{A}_{n}$, and hence, it does not
directly follow from the definition of contiguity that any distribution
function $\mathbb{Q}_{n}$ is contiguous with respect to $\mathbb{P}_{n}$. In
order to characterize the contiguity in terms of the asymptotic behavior of
the likelihood ratio between $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$, Le Cam
proposed some results which are popularly known as “Le Cam’s Lemma”. A
consequence of Le Cam’s first lemma is that the sequence $\mathbb{Q}_{n}$ will
be contiguous with respect to the sequence $\mathbb{P}_{n}$ if
$\log(\mathbb{Q}_{n}/\mathbb{P}_{n})$ is an asymptotically normal random
variable with mean $-\sigma^{2}/2$ and variance $\sigma^{2}$, where $\sigma$
is a positive constant. Moreover, the consequence of Le Cam’s third lemma is
that for $\mathbf{X}_{n}\in\mathbb{R}^{d}$, the joint distribution of
$\mathbf{X}_{n}$ and $\log(\mathbb{Q}_{n}/\mathbb{P}_{n})$ is distributed as
multivariate normal with mean vector
$(\boldsymbol{\mu},-\sigma^{2}/2)^{\text{T}}$ and covariance
$\begin{pmatrix}\boldsymbol{\Sigma}&\boldsymbol{\tau}\\\
\boldsymbol{\tau}^{\text{T}}&\sigma^{2}\end{pmatrix}$ under $\mathbb{P}_{n}$,
then under the alternative distribution $\mathbb{Q}_{n}$, the asymptotic
distribution of $\mathbf{X}$ is also a normal with mean
$\boldsymbol{\mu}+\boldsymbol{\tau}$ and covariance $\boldsymbol{\Sigma}$.

In order to test $H_{0}:F=F_{0}$, we consider the sequence of alternatives

$H_{n}:F_{n}=\left(1-\frac{\gamma}{\sqrt{n}}\right)F_{0}+\frac{\gamma}{\sqrt{n}}H$
(16)

for a fixed $\gamma>0$ and $n=1,2,\cdots$. In terms of depth function, this
testing of hypothesis problem can equivalently be written as
$H_{0}^{*}:D_{F}=D_{F_{0}}$, and the sequence of alternatives
$H_{n}^{*}:D_{F_{n}}=(1-\frac{\gamma}{\sqrt{n}})D_{F_{0}}+\frac{\gamma}{\sqrt{n}}D_{H}$.
It follows from the Propositions 1 and 2 that the aforesaid hypothesis
statement is valid for half-space depth function. Theorems 8 and 9 state the
asymptotic power properties of the proposed test.

######  Theorem 8.

Assume that $F_{0}$ and $H$ (see (16)) have continuous and positive densities
$f_{0}$ and $h$, respectively on $\mathbb{R}^{d}\,(d\geq 2)$ such that
$\mathbb{E}_{H_{0}}\left\\{\frac{h(\mathbf{x})}{f_{0}(\mathbf{x})}-1\right\\}^{4}<\infty$,
and suppose that the optimal half-space depth associated to $\mathbf{x}$ is
unique. In addition, conditions (C1)-(C5) hold. Then the sequence of
alternatives is a contiguous sequence. Further, assume that
$\mathcal{G}_{1}(\mathbf{x})$ is a random element associated with a Gaussian
process with mean function zero and covariance kernel
$F_{0}(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F_{0}(H[\mathbf{x}_{1}])F_{0}(H[\mathbf{x}_{2}])$, and
$\mathcal{G}_{1}^{\prime}(\mathbf{x})$ is a random element associated with a
Gaussian process with mean function $-\gamma\mathbb{E}_{\mathbf{x}\sim
h}\\{D_{F_{0}}(\mathbf{x})\\}$ and covariance kernel
$F_{0}(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F_{0}(H[\mathbf{x}_{1}])F_{0}(H[\mathbf{x}_{2}])$. Under
the alternatives described in (16), the asymptotic power of the test based on
$T_{\mathcal{X},F_{0}}^{\text{KS}}$ is
$\mathbb{P}_{\gamma}\left\\{\sup\limits_{\mathbf{x}}|\mathcal{G}_{1}^{\prime}(\mathbf{x})|>s_{1-\alpha}^{(1)}\right\\}$,
where $s_{1-\alpha}^{(1)}$ is such that
$\mathbb{P}_{\gamma=0}\left\\{\sup\limits_{\mathbf{x}}|\mathcal{G}_{1}(\mathbf{x})|>s_{1-\alpha}^{(1)}\right\\}=\alpha$.
Moreover, under the alternative hypothesis described in (16), the asymptotic
power of the test based on $T_{\mathcal{X},F_{0}}^{\text{CvM}}$ is
$\mathbb{P}_{\gamma}\left\\{\int\limits|\mathcal{G}_{1}^{\prime}(\mathbf{x})|^{2}dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}\right\\}$,
where $s_{1-\alpha}^{(2)}$ is such that
$\mathbb{P}_{\gamma=0}\left\\{\int\limits|\mathcal{G}_{1}(\mathbf{x})|^{2}dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}\right\\}=\alpha$.

Now consider the scenario of a two-sample problem. The null hypothesis is
given by $H_{0}:F=G$ against the sequences of alternatives

$H_{n,m}:G=\left(1-\frac{\gamma}{\sqrt{n+m}}\right)F+\frac{\gamma}{\sqrt{n+m}}H$
(17)

for a fixed $\gamma>0$ and $n,m=1,2,\cdots$, which is equivalent to test
$H_{0}^{*}:D_{F}=D_{G}$ against the sequence of alternatives
$H_{n,m}^{*}:D_{G}=\left(1-\frac{\gamma}{\sqrt{n+m}}\right)D_{F}+\frac{\gamma}{\sqrt{n+m}}D_{H}$.

######  Theorem 9.

Assume $F$ and $H$ (see (17)) have continuous and positive densities $f$ and
$h$, respectively on $\mathbb{R}^{d}\,(d\geq 2)$ such that
$\mathbb{E}_{F}\left\\{\frac{h(\mathbf{x})}{f(\mathbf{x})}-1\right\\}^{4}<\infty$.
Suppose that the optimal half-space depth associated to $\mathbf{x}$ is unique
and $\lim\limits_{\min(n,m)\rightarrow\infty}\frac{n}{m+n}=\lambda\,\in(0,1)$,
and in addition, conditions (C1)-(C5) hold. Then the sequence of alternatives
is a contiguous sequence. Furthermore, assume that
$\mathcal{G}_{2}(\mathbf{x})$ is a random element associated with a Gaussian
process with mean function zero and covariance kernel
$\left\\{F(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F(H[\mathbf{x}_{1}])F(H[\mathbf{x}_{2}])\right\\}/\lambda(1-\lambda)$
and $\mathcal{G}_{2}^{\prime}(\mathbf{x})$ is a random element associated with
a Gaussian process with mean function
$\gamma\sqrt{\lambda/(1-\lambda)}\mathbb{E}_{\mathbf{u}\sim
h}\left\\{D_{F}(\mathbf{x})\right\\}$ and covariance kernel
$\left\\{F(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F(H[\mathbf{x}_{1}])F(H[\mathbf{x}_{2}])\right\\}/\lambda(1-\lambda)$.
Under the alternatives described in (17), the asymptotic power of the test
based on KS is
$\mathbb{P}_{\gamma}\left\\{\sup_{\mathbf{x}}|\mathcal{G}_{2}^{\prime}(\mathbf{x})>t_{1-\alpha}\right\\}$
where $t_{1-\alpha}^{(1)}$ is such that
$\mathbb{P}_{\gamma=0}\\{\sup_{\mathbf{x}}|\mathcal{G}_{2}(\mathbf{x})|>t_{1-\alpha}^{(1)}\\}=\alpha$.
Moreover, under the alternative hypothesis described in (17), the asymptotic
power of the test based on CvM is
$\mathbb{P}_{\gamma}\left\\{\int|\mathcal{G}_{2}^{\prime}(\mathbf{x})|^{2}d\mathbf{x}>t^{(2)}_{1-\alpha}\right\\}$
where
$\mathbb{P}_{\gamma=0}\left\\{\int|\mathcal{G}_{2}^{\prime}(\mathbf{x})|^{2}d\mathbf{x}>t^{(2)}_{1-\alpha}\right\\}=\alpha$.

The results in Theorems 8 and 9 provide us the asymptotic power of the
proposed tests under the local alternatives described in (16) and (17),
respectively. Using those results, we compute the asymptotic powers of the
proposed tests for various choices of $\gamma$. In this study, we consider
that $F_{0}$ is the standard bivariate normal distribution, and $H$ is the
standard bivariate Laplace distribution. Figure 3 illustrate the summarised
results, and it is clearly indicated by those diagram that our proposed tests
perform well for this example. For the sake of concise presentation, we are
not reporting here the results for other choices of $F_{0}$ and $H$, and some
other choices of dimension. Nevertheless, our preliminary investigation
suggests that the tests perform well for alternative choices also.

(a) Goodness-of-fit test

(b) Two-sample test

Figure 3: The asymptotic power for (a) goodness-of-fit test and (b) two-sample
test with $\lambda=0.5$ under contiguous alternative described in (16) and
(17), respectively, where $F_{0}$ is standard bivariate normal distribution
and $H$ is the standard Laplace distribution. The black solid line indicates
the results based on proposed depth-based KS test statistics and the red
dotted line indicates the results based on proposed depth-based CvM test
statistics.

##  5 Simulation studies

In this section, we demonstrate the finite sample performance by reporting the
size and power of the proposed tests obtained from multiple studies compared
to other existing methods. In order to implement the proposed tests, we
prepared R codes which are available at https://github.com/pratimguhaniyogi.
Precisely speaking, we use mvtnorm and LaplacesDemon to generate the data from
multivariate normal, t and Laplace distributions, respectively. Next, in order
to compute the half-space depth, the package, viz., ddalpha is used. For other
details, the users are referred to the aforementioned web-link. Next, in order
to implement the other tests, mvnTest is used to run goodness-of-fit tests.
All the simulation results in the following examples are based on 500
replications, and the critical value of the test is estimated by 500 bootstrap
samples in each of the simulation runs. We also study the simulation results
for cases with different dimensions, choosing $d\in\\{2,6,10\\}$. Furthermore,
we denote $\boldsymbol{\mu}\in\mathbb{R}^{d}$ as a location parameter and
$\boldsymbol{\Sigma}$ is a $d\times d$ positive definite matrix as a scale
parameter in those examples.

In the goodness-of-fit problem (see Problem 1), we assume that $F_{0}$ is a
standard $d$-dimensional normal distribution. The choices of unknown
distributions described in Problem 1 are listed below.

  1. Model A.1.

Standard $d$-variate normal distribution
$N_{d}(\mathbf{0}_{d},\mathbf{I}_{d})$, where
$\mathbf{0}_{d}=(0,\cdots,0)^{\text{T}}$ is the $d$-dimensional null vector,
and $\mathbf{I}_{d}$ refers to the $d\times d$ identity matrix.

  2. Model A.2.

Mixture of the $d$-variate normal distribution with mixing probability $0.8$,
where the samples are taken from
$0.8N_{d}(\mathbf{0}_{d},\mathbf{I}_{d})+0.2N_{d}(5\times\textbf{1}_{d},\mathbf{I}_{d})$.

  3. Model A.3.

Mixture of the $d$-variate normal distribution with mixing probability $0.8$,
where the samples are taken as a from
$0.8N_{d}(\mathbf{0}_{d},\mathbf{I}_{d})+0.2N_{d}(\textbf{0}_{d},\boldsymbol{\Sigma})$.
Here, $\boldsymbol{\Sigma}$ be a compound symmetric matrix with off-diagonal
elements as $0.5$ and diagonal elements as $1$.

  4. Model A.4.

$d$-variate t-distribution with degrees of freedom 3,
$t(\mathbf{0}_{d},\mathbf{I}_{d},3)$.

  5. Model A.5.

Standard $d$-variate Cauchy distribution $C(\mathbf{0}_{d},\mathbf{I}_{d})$ or
$t(\mathbf{0}_{d},\mathbf{I}_{d},1)$ with p.d.f.
$f(\mathbf{x})=[\Gamma((d+1)/2)/(\sqrt{\pi}\Gamma(d/2))]\times(1+\|\mathbf{x}\|^{2})^{-(d+1)/2}$.

  6. Model A.6.

Standard $d$-variate Laplace distribution $L(\mathbf{0}_{d},\mathbf{I}_{d})$
with p.d.f.
$f(\mathbf{x})=[\Gamma(d/2)/(2\Gamma(d)\pi^{d/2})]\times\exp(-\|\mathbf{x}\|)$.

We fix the sample size $n\in\\{10,25,50,100,300,500\\}$ and the nominal
significance level $\alpha$ is assumed to be $0.05$. We compare the proposed
method denoted as KS.depth and CvM.depth respectively with Anderson-Darling
(AD), Cramér-von Mises (CvM), Doornik-Hansen (DH), Henze-Zirkler (HZ), Royston
(R) tests and a series of $\chi^{2}$-type tests such as McCulloch (M),
Nikulin-Rao-Robson (NRR), Dzhaparidze-Nikulin (DN) tests.

Table 1: Goodness-of-fit test: Observed relative frequency for rejecting the null, $d=2$. $n$ | KS.depth | CvM.depth | AD | CM | DH | HZ | R | McCulloch | NRR | DN
---|---|---|---|---|---|---|---|---|---|---
Model A.1. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$
10 | 0.074 | 0.040 | 0.060 | 0.050 | 0.052 | 0.020 | 0.040 | 0.032 | 0.050 | 0.062
25 | 0.040 | 0.028 | 0.034 | 0.030 | 0.044 | 0.052 | 0.044 | 0.030 | 0.044 | 0.054
50 | 0.038 | 0.056 | 0.060 | 0.058 | 0.060 | 0.042 | 0.056 | 0.044 | 0.048 | 0.048
100 | 0.052 | 0.058 | 0.042 | 0.056 | 0.048 | 0.046 | 0.044 | 0.044 | 0.060 | 0.062
300 | 0.055 | 0.090 | 0.060 | 0.055 | 0.065 | 0.050 | 0.055 | 0.080 | 0.055 | 0.050
500 | 0.070 | 0.066 | 0.052 | 0.050 | 0.072 | 0.038 | 0.060 | 0.038 | 0.064 | 0.054
Model A.2. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim
0.8N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})+0.2N_{2}(5\textbf{1}_{2},\mathbf{I}_{2})$
10 | 0.600 | 0.800 | 0.000 | 0.000 | 0.000 | 0.000 | 0.600 | 0.000 | 0.200 | 0.200
25 | 0.751 | 0.788 | 0.168 | 0.157 | 0.129 | 0.980 | 0.996 | 0.105 | 0.122 | 0.090
50 | 0.886 | 0.930 | 0.212 | 0.214 | 0.198 | 1.000 | 1.000 | 0.144 | 0.158 | 0.122
100 | 0.994 | 1.000 | 0.328 | 0.334 | 0.488 | 1.000 | 1.000 | 0.118 | 0.250 | 0.224
300 | 1.000 | 1.000 | 0.656 | 0.662 | 0.982 | 1.000 | 1.000 | 0.134 | 0.668 | 0.684
500 | 1.000 | 1.000 | 0.810 | 0.818 | 1.000 | 1.000 | 1.000 | 0.132 | 0.832 | 0.822
Model A.3. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim
0.8N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})+0.2N_{2}(\mathbf{0}_{2},\boldsymbol{\Sigma})$
10 | 0.600 | 0.600 | 0.400 | 0.200 | 0.200 | 0.000 | 0.200 | 0.200 | 0.200 | 0.200
25 | 0.469 | 0.611 | 0.050 | 0.035 | 0.048 | 0.046 | 0.050 | 0.046 | 0.046 | 0.066
50 | 0.532 | 0.650 | 0.052 | 0.042 | 0.056 | 0.056 | 0.054 | 0.036 | 0.048 | 0.040
100 | 0.588 | 0.732 | 0.046 | 0.050 | 0.054 | 0.048 | 0.040 | 0.048 | 0.050 | 0.056
300 | 0.644 | 0.808 | 0.082 | 0.068 | 0.064 | 0.078 | 0.060 | 0.064 | 0.068 | 0.056
500 | 0.690 | 0.844 | 0.054 | 0.060 | 0.062 | 0.052 | 0.064 | 0.084 | 0.048 | 0.040
Model A.4. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim t(\mathbf{0}_{2},\mathbf{I}_{2},3)$
10 | 0.296 | 0.300 | 0.104 | 0.072 | 0.172 | 0.176 | 0.178 | 0.026 | 0.040 | 0.054
25 | 0.412 | 0.600 | 0.478 | 0.444 | 0.524 | 0.462 | 0.456 | 0.396 | 0.356 | 0.174
50 | 0.702 | 0.900 | 0.822 | 0.818 | 0.822 | 0.772 | 0.758 | 0.788 | 0.756 | 0.430
100 | 0.956 | 0.994 | 0.994 | 0.992 | 0.986 | 0.978 | 0.972 | 0.988 | 0.986 | 0.810
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.5. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim C(\mathbf{0}_{2},\mathbf{I}_{2})$
10 | 0.728 | 0.790 | 0.476 | 0.438 | 0.588 | 0.708 | 0.550 | 0.190 | 0.258 | 0.234
25 | 0.974 | 0.988 | 0.982 | 0.974 | 0.954 | 0.988 | 0.960 | 0.922 | 0.948 | 0.860
50 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.998 | 1.000 | 0.996
100 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.6. $H_{0}:\mathbf{X}\sim N_{2}(\mathbf{0}_{2},\mathbf{I}_{2})$ against
$H_{1}:\mathbf{X}\sim L(\mathbf{0}_{2},\mathbf{I}_{2})$
10 | 0.272 | 0.252 | 0.116 | 0.078 | 0.194 | 0.158 | 0.168 | 0.022 | 0.040 | 0.066
25 | 0.422 | 0.670 | 0.550 | 0.474 | 0.520 | 0.520 | 0.458 | 0.260 | 0.302 | 0.224
50 | 0.720 | 0.960 | 0.908 | 0.882 | 0.786 | 0.842 | 0.752 | 0.622 | 0.738 | 0.514
100 | 0.946 | 0.998 | 0.992 | 0.992 | 0.950 | 0.992 | 0.934 | 0.914 | 0.980 | 0.878
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Table 2: Goodness-of-fit test: Observed relative frequency for rejecting the null, $d=6$. $n$ | KS.depth | CvM.depth | AD | CM | DH | HZ | R | McCulloch | NRR | DN
---|---|---|---|---|---|---|---|---|---|---
Model A.1. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$
10 | 0.072 | 0.042 | 0.066 | 0.062 | 0.054 | 0.004 | 0.060 | 0.008 | 0.376 | 0.610
25 | 0.066 | 0.052 | 0.050 | 0.070 | 0.050 | 0.030 | 0.034 | 0.046 | 0.084 | 0.092
50 | 0.064 | 0.052 | 0.060 | 0.052 | 0.034 | 0.038 | 0.034 | 0.048 | 0.060 | 0.066
100 | 0.046 | 0.050 | 0.056 | 0.066 | 0.044 | 0.046 | 0.036 | 0.064 | 0.070 | 0.068
300 | 0.039 | 0.072 | 0.083 | 0.072 | 0.056 | 0.028 | 0.050 | 0.050 | 0.061 | 0.078
500 | 0.038 | 0.068 | 0.038 | 0.038 | 0.050 | 0.060 | 0.062 | 0.048 | 0.046 | 0.044
Model A.2. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim
0.8N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})+0.2N_{6}(5\textbf{1}_{6},\mathbf{I}_{6})$
10 | 0.941 | 1.000 | 0.059 | 0.000 | 0.059 | 0.000 | 0.706 | 0.000 | 0.412 | 0.529
25 | 1.000 | 1.000 | 0.083 | 0.083 | 0.042 | 0.292 | 1.000 | 0.042 | 0.083 | 0.125
50 | 0.996 | 1.000 | 0.048 | 0.040 | 0.068 | 0.920 | 1.000 | 0.064 | 0.054 | 0.064
100 | 1.000 | 1.000 | 0.094 | 0.086 | 0.088 | 1.000 | 1.000 | 0.084 | 0.096 | 0.082
300 | 1.000 | 1.000 | 0.150 | 0.136 | 0.130 | 1.000 | 1.000 | 0.054 | 0.104 | 0.106
500 | 1.000 | 1.000 | 0.234 | 0.186 | 0.184 | 1.000 | 1.000 | 0.050 | 0.140 | 0.150
Model A.3. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim
0.8N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})+0.2N_{6}(\mathbf{0}_{6},\boldsymbol{\Sigma})$
10 | 0.882 | 1.000 | 0.059 | 0.059 | 0.059 | 0.000 | 0.059 | 0.000 | 0.176 | 0.529
25 | 0.958 | 1.000 | 0.042 | 0.042 | 0.000 | 0.000 | 0.000 | 0.042 | 0.000 | 0.083
50 | 0.910 | 0.962 | 0.032 | 0.036 | 0.064 | 0.046 | 0.044 | 0.038 | 0.060 | 0.060
100 | 0.930 | 0.970 | 0.054 | 0.054 | 0.040 | 0.088 | 0.044 | 0.044 | 0.056 | 0.052
300 | 0.952 | 0.986 | 0.152 | 0.124 | 0.046 | 0.092 | 0.044 | 0.090 | 0.078 | 0.068
500 | 0.976 | 0.992 | 0.300 | 0.240 | 0.074 | 0.152 | 0.050 | 0.158 | 0.168 | 0.124
Model A.4. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim t(\mathbf{0}_{6},\mathbf{I}_{6},3)$
10 | 0.116 | 0.154 | 0.018 | 0.016 | 0.080 | 0.054 | 0.066 | 0.056 | 0.240 | 0.354
25 | 0.640 | 0.836 | 0.646 | 0.492 | 0.572 | 0.742 | 0.470 | 0.272 | 0.404 | 0.284
50 | 0.904 | 0.994 | 0.996 | 0.994 | 0.966 | 0.994 | 0.908 | 0.924 | 0.980 | 0.860
100 | 0.994 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.998 | 1.000 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.5. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim C(\mathbf{0}_{6},\mathbf{I}_{6})$
10 | 0.308 | 0.498 | 0.024 | 0.004 | 0.146 | 0.314 | 0.124 | 0.308 | 0.162 | 0.078
25 | 0.998 | 0.998 | 0.996 | 0.992 | 0.990 | 0.998 | 0.968 | 0.880 | 0.976 | 0.946
50 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.998 | 1.000 | 1.000
100 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.6. $H_{0}:\mathbf{X}\sim N_{6}(\mathbf{0}_{6},\mathbf{I}_{6})$ against
$H_{1}:\mathbf{X}\sim L(\mathbf{0}_{6},\mathbf{I}_{6})$
10 | 0.114 | 0.168 | 0.024 | 0.014 | 0.084 | 0.100 | 0.076 | 0.084 | 0.294 | 0.344
25 | 0.656 | 0.912 | 0.828 | 0.590 | 0.568 | 0.896 | 0.438 | 0.080 | 0.400 | 0.418
50 | 0.910 | 0.996 | 1.000 | 0.990 | 0.942 | 1.000 | 0.854 | 0.682 | 0.966 | 0.930
100 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.982 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Table 3: Goodness-of-fit test: Observed relative frequency for rejecting the null, $d=10$. $n$ | KS.depth | CvM.depth | AD | CM | DH | HZ | R | McCulloch | NRR | DN
---|---|---|---|---|---|---|---|---|---|---
Model A.1. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
25 | 0.050 | 0.050 | 0.038 | 0.042 | 0.038 | 0.028 | 0.050 | 0.034 | 0.158 | 0.168
50 | 0.052 | 0.042 | 0.036 | 0.040 | 0.066 | 0.040 | 0.052 | 0.060 | 0.088 | 0.092
100 | 0.054 | 0.062 | 0.050 | 0.052 | 0.052 | 0.052 | 0.052 | 0.056 | 0.080 | 0.078
300 | 0.040 | 0.023 | 0.051 | 0.040 | 0.069 | 0.069 | 0.057 | 0.057 | 0.040 | 0.040
500 | 0.050 | 0.060 | 0.058 | 0.052 | 0.056 | 0.054 | 0.060 | 0.042 | 0.048 | 0.050
Model A.2. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim
0.8N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})+0.2N_{10}(5\textbf{1}_{10},\mathbf{I}_{10})$
25 | 0.995 | 0.997 | 0.032 | 0.037 | 0.037 | 0.099 | 1.000 | 0.048 | 0.136 | 0.171
50 | 0.998 | 1.000 | 0.042 | 0.042 | 0.062 | 0.330 | 1.000 | 0.072 | 0.096 | 0.094
100 | 1.000 | 1.000 | 0.051 | 0.042 | 0.051 | 0.850 | 1.000 | 0.060 | 0.066 | 0.060
300 | 1.000 | 1.000 | 0.050 | 0.054 | 0.082 | 1.000 | 1.000 | 0.046 | 0.070 | 0.080
500 | 1.000 | 1.000 | 0.093 | 0.086 | 0.064 | 1.000 | 1.000 | 0.061 | 0.068 | 0.068
Model A.3. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim
0.8N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})+0.2N_{10}(\mathbf{0}_{10},\boldsymbol{\Sigma})$
25 | 0.979 | 1.000 | 0.019 | 0.024 | 0.051 | 0.043 | 0.043 | 0.061 | 0.131 | 0.120
50 | 0.976 | 1.000 | 0.014 | 0.010 | 0.060 | 0.052 | 0.060 | 0.032 | 0.032 | 0.034
100 | 0.978 | 0.996 | 0.040 | 0.028 | 0.068 | 0.078 | 0.054 | 0.068 | 0.060 | 0.056
300 | 0.934 | 0.986 | 0.288 | 0.212 | 0.074 | 0.160 | 0.056 | 0.084 | 0.120 | 0.128
500 | 0.894 | 0.982 | 0.548 | 0.430 | 0.068 | 0.310 | 0.054 | 0.122 | 0.260 | 0.238
Model A.4. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim t(\mathbf{0}_{10},\mathbf{I}_{10},3)$
25 | 0.330 | 0.708 | 0.296 | 0.098 | 0.384 | 0.676 | 0.272 | 0.054 | 0.160 | 0.198
50 | 0.830 | 0.996 | 1.000 | 0.990 | 0.978 | 1.000 | 0.912 | 0.886 | 0.974 | 0.922
100 | 0.998 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.5. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim C(\mathbf{0}_{10},\mathbf{I}_{10})$
25 | 0.930 | 0.996 | 0.970 | 0.904 | 0.908 | 0.998 | 0.844 | 0.092 | 0.848 | 0.854
50 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.998 | 1.000 | 1.000
100 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
Model A.6. $H_{0}:\mathbf{X}\sim N_{10}(\mathbf{0}_{10},\mathbf{I}_{10})$
against $H_{1}:\mathbf{X}\sim L(\mathbf{0}_{10},\mathbf{I}_{10})$
25 | 0.436 | 0.900 | 0.740 | 0.312 | 0.446 | 0.948 | 0.328 | 0.224 | 0.470 | 0.394
50 | 0.894 | 1.000 | 1.000 | 0.998 | 0.950 | 1.000 | 0.890 | 0.430 | 0.968 | 0.966
100 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.998 | 0.978 | 1.000 | 1.000
300 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
500 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000

Tables 1, 2 and 3 describe Monte Carlo results for goodness-of-fit tests that
represents the observed relative frequency for rejecting $H_{0}$ for $d=2,6$
and $10$ respectively. Here, the observed relative frequency for rejecting the
null for Model A.1 indicates the size of the test and for Model A.2-A.6
represent the power of the test for different alternatives. If we increase the
sample size $n$ the power of the test increases for all the methods and the
size of the test are more or less stabilized at the nominal significance
level. If we increase the dimension of the data, the aforementioned statement
holds true. The power of the proposed tests are consistently better for all
the choices of alternative distributions. Moreover, for mixture distributions
like Model A.2 and A.3, the competing testing procedures are worse, whereas
the proposed methods KS.depth and CvM.depth perform satisfactorily. The Figure
4 represents receiver operating characteristics (ROC) curve that illustrates
the performance of testing procedure. The ROC curve is created by plotting the
power as a function of the Type I error of the testing procedure. For brevity,
here we represent the ROC curve for the goodness-of-fit test when the
alternative hypotheses are (a) Cauchy and (d) Laplace distributions for small
sample size. Here, this figure tells us that our proposed testing procedures
(especially the data-depth based CvM method) outperform the existing methods.

(a) $H_{1}:\mathbf{X}\sim C(\mathbf{0}_{d},\mathbf{I}_{d})$

(b) $H_{1}:\mathbf{X}\sim L(\mathbf{0}_{d},\mathbf{I}_{d})$

Figure 4: ROC plot for the goodness-of-fit test when the alternative
distributions are standard Cauchy (a) and standard Laplace distribution (b),
respectively, where the null distribution is standard normal distribution.

In the two-sample test scenario (see Problem 2), we use sample sizes $(n,m)$
such that $\rho=n/(n+m)\in\\{0.3,0.5,0.8\\}$ and fix $n\in\\{50,100,300\\}$.
We generate samples from the following distribution family.

  1. Model B.

The first samples are taken from the standard multivariate distribution
$N_{d}(\mathbf{0}_{d},\mathbf{I}_{d})$ and the second sample is taken from the
$d$ variate normal $N_{d}(\mu\times\textbf{1}_{d},\mathbf{I}_{d})$ where we
chose $\mu\in\\{0,0.5,01\\}$.

Table 4 represents the observed relative frequency for rejecting $H_{0}$ based
on the above model. The values on the table corresponding to $\mu=0$ is the
size of the test and that to $\mu=\mu_{0}\in\\{0.5,1\\}$ represent the power
of the test at $\mu_{0}$. If $\mu$ creases, so do the powers of the tests.
Even if the dimension of the data is increased, the performance of the
proposed test is not affected.

Table 4: Two-sample test: Observed relative frequency for rejecting the null based on Model B. |  |  | $\mu=0$ | $\mu=0.5$ | $\mu=1$
---|---|---|---|---|---
| n | m | KS.depth | CvM.depth | KS.depth | CvM.depth | KS.depth | CvM.depth
$d=2$
$\rho=0.3$ | 50 | 117 | 0.054 | 0.052 | 0.908 | 0.940 | 1.000 | 1.000
100 | 234 | 0.050 | 0.058 | 0.994 | 0.998 | 1.000 | 1.000
300 | 700 | 0.061 | 0.046 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.5$ | 50 | 50 | 0.064 | 0.064 | 0.790 | 0.836 | 0.998 | 1.000
100 | 100 | 0.066 | 0.068 | 0.974 | 0.986 | 1.000 | 1.000
300 | 300 | 0.048 | 0.058 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.8$ | 50 | 13 | 0.092 | 0.076 | 0.416 | 0.452 | 0.904 | 0.926
100 | 25 | 0.064 | 0.060 | 0.676 | 0.780 | 0.996 | 1.000
300 | 75 | 0.058 | 0.054 | 0.996 | 0.998 | 1.000 | 1.000
$d=6$
$\rho=0.3$ | 50 | 117 | 0.072 | 0.070 | 0.996 | 0.986 | 1.000 | 1.000
100 | 234 | 0.074 | 0.064 | 1.000 | 1.000 | 1.000 | 1.000
300 | 700 | 0.056 | 0.060 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.5$ | 50 | 50 | 0.132 | 0.058 | 0.998 | 0.996 | 1.000 | 1.000
100 | 100 | 0.052 | 0.052 | 1.000 | 1.000 | 1.000 | 1.000
300 | 300 | 0.048 | 0.056 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.8$ | 50 | 13 | 0.290 | 0.206 | 0.914 | 0.892 | 0.986 | 0.998
100 | 25 | 0.132 | 0.102 | 0.976 | 0.988 | 1.000 | 1.000
300 | 75 | 0.068 | 0.064 | 1.000 | 1.000 | 1.000 | 1.000
$d=10$
$\rho=0.3$ | 50 | 117 | 0.086 | 0.078 | 1.000 | 0.997 | 1.000 | 1.000
100 | 234 | 0.066 | 0.070 | 1.000 | 1.000 | 1.000 | 1.000
300 | 700 | 0.031 | 0.073 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.5$ | 50 | 50 | 0.106 | 0.048 | 1.000 | 1.000 | 1.000 | 1.000
100 | 100 | 0.068 | 0.062 | 1.000 | 1.000 | 1.000 | 1.000
300 | 300 | 0.052 | 0.050 | 1.000 | 1.000 | 1.000 | 1.000
$\rho=0.8$ | 50 | 13 | 0.486 | 0.264 | 0.976 | 0.976 | 0.990 | 0.996
100 | 25 | 0.140 | 0.116 | 1.000 | 0.998 | 1.000 | 1.000
300 | 75 | 0.066 | 0.064 | 1.000 | 1.000 | 1.000 | 1.000

##  6 Real data analysis

To understand the practicability of the proposed tests, two well-known data
sets are analyzed here.

##### Fisher’s Iris data

This data set, available in the base R, consists of three multivariate samples
corresponding to three different species of Iris, namely Iris setosa, Iris
virginica and Iris versicolor with sample size 50 for each species. In each
species, the length and width of the sepals are measured in centimeters. We
would like to determine how close the distribution of each sample associated
with sepal length and sepal width is to a bivariate normal sample. This can be
formulated as a goodness-of-fit testing problem, with $F_{0}$ being a
bivariate normal distribution with unknown mean $\boldsymbol{\mu}$ and unknown
dispersion matrix $\boldsymbol{\Sigma}$. For each species, we estimate these
unknown parameters using the corresponding mean vector and covariance matrix.
The data-depth discrepancy is graphically represented for three different
species in Figure 5, where we observe that most points are tightly clustered
around the straight line, indicating that the standardized data are from a
standard bivariate normal distribution. Moreover, our goodness-of-fit test for
testing $H_{0}:F=F_{0}$ against $H_{1}:F\neq F_{0}$ led to very high empirical
p-values for all types of spices (see Table 5). This indicates that $H_{0}$ is
to be accepted, and therefore, bivariate normal distributions show signs of
good fit for the data for the sepal length and sepal width of three Iris
species.

(a) Iris setosa

(b) Iris virginica

(c) Iris versicolor

Figure 5: Goodness-of-fit test for Iris data: The data-depth discrepancy plot for Iris setosa, Iris virginica and Iris versicolor respectively for testing normality. The dotted gray curves indicate the two-sigma limits of data-depth discrepancy Table 5: Empirical p-values of the proposed goodness-of-fit tests obtained over 1000 replications. | p-value of depth based tests
---|---
Species | KS.depth | CvM.depth
Iris setosa | 0.286 | 0.418
Iris virginica | 0.43 | 0.39
Iris versicolor | 0.19 | 0.118

#####  gilgais data

This data set is available in the MASS package in R. It is collected on a line
transect survey in the gilgai territory in New South Wales, Australia
(Webster, 1977). On a 4-meter spaced linear grid, 365 sampling locations are
selected. Electrical conductivity (in mS / cm), pH, and chloride contained (in
ppm) are collected at three different depths below the surface, 0-10 cm, 30-40
cm, and 80-90 cm. We would like to investigate how close the joint
distribution of three variables at each level is to a trivariate normal
distribution. As the previous example, this can be formulated as a goodness-
of-fit problem, where $F_{0}$ is specified as a trivariate normal distribution
with unknown mean $\mu$ and unknown dispersion matrix $\Sigma$. We estimate
these unknown parameters using their respective maximum likelihood estimates.
We then standardize the data in each sample using the corresponding mean
vector and covariance matrix. The proposed data-depth dispersion plot is shown
in Figure 6 where we observe a majority of the data cloud is far from the
straight line. This phenomenon indicates that the data are not from a
trivariate normal distribution. Moreover, the test of $H_{0}:F=F_{0}$ against
$H_{1}:F\neq F_{0}$ led to zero empirical p-values for all height levels below
the surface. Moreover, we are interested to study whether distribution of any
height is same as the other. Therefore, this can be viewed as two-sample test
where $F$ is the distribution of chemical parameters art height at $0-10$ cm
and G is that of either at height $30–40$ cm or $80–90$ cm, respectively. The
Figure 7(b) shows that the data cloud is far from the straight line, this
phenomenon indicates that the distributions are different. Moreover, the
empirical p-value for each of situations are close to zero based on the
proposed test which determines that the distributions are significantly
different.

(a) height 0–10 cm

(b) height 30–40 cm

(c) height 80–90 cm

Figure 6: Goodness-of-fit test for gilgais data: The data-depth discrepancy
plot for tuple at depths below the surface of height 0-10 cm, 30-40 cm and
80-90 cm respectively for testing normality. The dotted gray curves indicate
the two-sigma limits of data-depth discrepancy

(a) heights 0–10 cm against 30-40 cm

(b) height 0–10 cm against 30-40 cm

Figure 7: Two-sample test for gilgias data: The data-depth discrepancy plot
for heights 0–10 cm against 30-40 cm and height 0–10 cm against 30-40 cm are
shown. The dotted gray curves indicate the two-sigma limits of data-depth
discrepancy

##  7 Conclusion

In this paper, we have proposed a data-depth discrepancy and a graphical tool-
kit based on Tukey’s half-space depth. This device is used to test the
equality of two distribution functions and is applicable to any dimension of
the distributions. The motivations behind the graphical device are shown
through simulation examples. Influenced by the graphical device, we have
proposed test statistics based on the Kolmogorov-Smirnov and Cramér-von Mises
tests. For goodness-of-fit and two-sample testing problems, the proposed test
statistics can test the equality of two distribution functions, which are
common in practice. We have shown that the test statistics are point-wise and
also uniformly constant. Moreover, we have studied the asymptotic power under
contiguous alternatives. The applicability of the proposed method is
illustrated by simulation studies and real data analysis. The graphical tool
and test statistics developed in this article can be used regardless of the
dimension of the data and therefore, significantly contributes to statistics
and machine learning literature.

Acknowledgement : Subhra Sankar Dhar is thankful to the research grants
MTR/2019/000039 and CRG/2022/001489, Government of India.

## Appendix

###### Proof of Theorem 1.

Note that, for every $\epsilon>0$,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left\\{\widehat{\mathcal{C}}(\mathcal{X},F_{0})\subset\mathcal{C}_{\epsilon}(F,F_{0})\right\\}\geq\mathbb{P}\left\\{\sup_{\mathbf{x}}|D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|<\epsilon\right\\}\rightarrow
1$ (18)

due to the Proposition 4.4 of Massé (2004). ∎

###### Proof of Theorem 2.

By the similar argument of Theorem 2, using the fact that
$\lim\limits_{\min(n,m)\rightarrow\infty}n/(n+m)=\lambda\in(0,1)$, observe
that,

$\displaystyle\lim_{\min(n,m)\rightarrow\infty}\mathbb{P}\\{\widehat{\mathcal{C}}(\mathcal{X},\mathcal{Y})\subset\mathcal{C}_{\epsilon}(F,G)\\}\geq\mathbb{P}\left\\{\sup_{\mathbf{x}}|D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x})|<\epsilon\right\\}\rightarrow
1$

∎

###### Proof of Theorem 3.

Let us define the Kolmogorv-Smirnov distance between empirical and theoretical
distribution functions $\widehat{F}_{n}$ and $F_{0}$ based on Tukey’s half-
space depth $D$ as
$d_{K}(\mathcal{X},F_{0})=\sup\limits_{\mathbf{x}\in\mathbb{R}^{d}}|D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|=\sup\limits_{\mathbf{x},\mathbf{u}}|\widehat{F}_{n}(H[\mathbf{x},\mathbf{u}])-F_{0}(H[\mathbf{x},\mathbf{u}])|$
where $\mathcal{X}=\\{\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\\}$ and $H$ is the
half-space with
$H[\mathbf{x},\mathbf{u}]=\\{\mathbf{y}\in\mathbb{R}^{d}:\mathbf{u}^{\text{T}}\mathbf{x}\geq\mathbf{u}^{\text{T}}\mathbf{y}\\}$.
It is important to note that $d_{K}(\cdot,\cdot)$ has the Glivenko-Cantalli
property (Pollard, 1984; Donoho and Gasko, 1992), i.e., if the data samples
are i.i.d. from $F_{0}$, $d_{K}(\mathcal{X},F_{0})\rightarrow 0$ as
$n\rightarrow\infty$. Therefore, under an alternative $F$, we have

$\sup_{\mathbf{x}}|D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|\rightarrow
d_{K}(D_{F},D_{F_{0}})>0$ (19)

almost surely. Now fix any $F$ with $d_{K}(D_{F},D_{F_{0}})>0$, then there
exists some $\mathbf{x}$ with $D_{\mathcal{X}}(\mathbf{x})\neq
D_{F}(\mathbf{x})$. Without loss of generality, assume that
$D_{\mathcal{X}}(\mathbf{x})>D_{F_{0}}(\mathbf{x})$. Then

$\displaystyle\mathbb{P}_{F}\left\\{T_{\mathcal{X},F_{0}}^{\text{KS}}>s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F}\left\\{\sqrt{n}|D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|>s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle=\mathbb{P}_{F}\left\\{\sqrt{n}[D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})+D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x})]>s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F}\left\\{\sqrt{n}[D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})]>s_{1-\alpha}^{(1)}-\sqrt{n}[D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x})]\right\\}$
$\displaystyle\rightarrow 1\text{ as }n\rightarrow\infty.$ (20)

The last implication follows from the fact that
$\left\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\right\\}$ is
uniformly bounded in probability in view of Proposition 4 and an application
of Prokhorov’s theorem, and for finite $s_{1-\alpha}^{(1)}$, we have
$s_{1-\alpha}^{(1)}-\sqrt{n}|D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x})|\rightarrow-\infty$
as $n\rightarrow\infty$. Hence the limiting power is one when
$D_{F}(\mathbf{x})>D_{F_{0}}(\mathbf{x})$. Similar argument can show that the
limiting power is one when $D_{F}(\mathbf{x})<D_{F_{0}}(\mathbf{x})$ for some
$\mathbf{x}$. Hence, the proof is complete. ∎

###### Proof of Theorem 4.

Due to the Glivenko-Cantalli property of $d_{K}(\cdot,\cdot)$, under an
alternative $F$, 19 holds almost surely. Thus fix any $F$with
$d_{K}(D_{F},D_{F_{0}})>0$, then there exits some $\mathbf{x}$ with
$D_{\mathcal{X}}(\mathbf{x})\neq D_{F}(\mathbf{x})$. Therefore, observe that

$\displaystyle\mathbb{P}_{F}\left\\{T_{\mathcal{X},F_{0}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle=\mathbb{P}_{F}\left\\{n\int(D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x}))^{2}dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F}\left\\{n\int(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))^{2}dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}-a_{n}+2b_{n}\right\\}$
$\displaystyle\rightarrow 1\text{ as }n\rightarrow\infty,$ (21)

where
$a_{n}=n\int(D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x}))^{2}dF_{0}(\mathbf{x})$
and
$b_{n}=n\int(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))(D_{F}(\mathbf{x})-D_{F_{0}}(\mathbf{x}))dF_{0}(\mathbf{x})$.
The last implication follows from the following facts: (i)
$\left\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\right\\}$ is
uniformly bounded in probability in view of Proposition 4, (ii)
$s_{1-\alpha}^{(2)}$ is positive finite, and (iii)
$\sup\limits_{\mathbf{x}}|D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})|\rightarrow
0$ as $n\rightarrow\infty$ under $F$ almost surely (see Donoho and Gasko
(1992) and Proposition 4.4 of Massé (2004)). An application of Prokhorov’s
theorem shows that the limiting power is one when
$D_{F}(\mathbf{x})>D_{F_{0}}(\mathbf{x})$. Similar argument can show that the
limiting power is one when $D_{F}(\mathbf{x})<D_{F_{0}}(\mathbf{x})$ for some
$\mathbf{x}$. Hence, the proof is complete. ∎

###### Proof of Theorem 5.

Let $F_{n}$ be any distribution function that satisfies
$\sqrt{n}d_{K}(D_{F_{n}},D_{F_{0}})\geq\Delta_{n}$. By the triangle
inequality, we have

$d_{K}(D_{F_{n}},D_{F_{0}})\leq
d_{K}(D_{F_{n}},D_{\mathcal{X}})+d_{K}(D_{\mathcal{X}},D_{F_{0}}),$ (22)

which implies
$T_{\mathcal{X},F_{0}}^{\text{KS}}\geq\Delta_{n}-\sqrt{n}d_{K}(D_{F_{n}},D_{F_{0}})$.
Therefore,

$\displaystyle\mathbb{P}_{F_{n}}\left\\{T_{\mathcal{X},F_{0}}^{\text{KS}}>s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n}}\left\\{\Delta_{n}-\sqrt{n}d_{K}(D_{F_{n}},D_{\mathcal{X}})>s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n}}\left\\{\sqrt{n}d_{K}(D_{F_{n}},D_{\mathcal{X}})\leq\Delta_{n}-s_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\rightarrow 1\text{ as }n\rightarrow\infty.$ (23)

The last implication follows from the following facts: (i)
$\left\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\right\\}$ is
uniformly bounded under $F_{n}$ in probability in view of Proposition 4 (ii)
$\Delta_{n}\rightarrow\infty$, and (iii) $s_{1-\alpha}^{(1)}$ is positive and
finite. An application of Prokhorov’s theorem, we get
$\mathbb{P}_{F_{n}}\left\\{T_{\mathcal{X},F_{0}}^{\text{KS}}>s_{1-\alpha}^{(1)}\right\\}\rightarrow
1$. ∎

###### Proof of Theorem 6.

Let $F_{n}$ be any distribution function satisfies
$\sqrt{n}d_{K}(D_{F_{n}},D_{F_{0}})\geq\Delta_{n}$. Therefore, using the
triangle inequality 22, we have

$\displaystyle\mathbb{P}_{F_{n}}\left\\{T_{\mathcal{X},F_{0}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n}}\left\\{\int\Delta_{n}^{2}dF_{0}(\mathbf{x})-n\int
d_{K}^{2}(D_{F_{n}},D_{\mathcal{X}})dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n}}\left\\{n\int
d_{K}^{2}(D_{F_{n}},D_{\mathcal{X}})dF_{0}(\mathbf{x})\leq\int\Delta_{n}^{2}dF_{0}(\mathbf{x})-s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\rightarrow 1\text{ as }n\rightarrow\infty.$ (24)

The last implication follows from the the following fact: (i)
$\left\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\right\\}$ is
uniformly bounded under $F_{n}$ in probability in view of Proposition 4 (ii)
$\Delta_{n}\rightarrow\infty$, and (iii) $s_{1-\alpha}^{(2)}$ is positive and
finite. An application of Prokhorov’s theorem, we get
$\mathbb{P}_{F_{n}}\left\\{T_{\mathcal{X},F_{0}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}\rightarrow
1$. ∎

###### Proof of Theorem 7.

Define $a_{n,m}=\sqrt{n+m}\times
d_{K}(D_{F},D_{G})=\sqrt{n+m}|D_{F}(\mathbf{x})-D_{G}(\mathbf{x})|$.

$\displaystyle\mathbb{P}_{H_{1}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}>t_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{H_{1}}\left\\{\sqrt{n+m}|D_{\mathcal{X}}(\mathbf{x})-D_{\mathcal{Y}}(\mathbf{x})|>t_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{H_{1}}\left\\{\sqrt{n+m}|D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})|+\sqrt{n+m}|D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x})|>t_{1-\alpha}^{(1)}-a_{n,m}\right\\}$
$\displaystyle\geq\mathbb{P}_{H_{1}}\left\\{\lambda^{-1/2}\sqrt{n}|D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})|+(1-\lambda)^{-1/2}\sqrt{m}|D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x})|>t_{1-\alpha}^{(1)}-a_{n,m}\right\\}$
$\displaystyle\rightarrow 1\text{ as }\min(n,m)\rightarrow\infty.$ (25)

The last implication follows from the fact that
$\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x})\\}$ and
$\\{\sqrt{m}(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x})\\}$ are uniformly
bounded in probability in view of Proposition 4, for finite
$t_{\alpha}^{(1)}$, we have
$\sqrt{n+m}|D_{F}(\mathbf{x})-D_{G}(\mathbf{x})|\rightarrow\infty$ as
$\min(n,m)\rightarrow\infty$. By an application of Prokhorov’s theorem, we get
$\mathbb{P}_{H_{1}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}>t_{1-\alpha}^{(1)}\right\\}\rightarrow
1$ as $\min(n,m)\rightarrow\infty$.

Furthermore,

$\displaystyle\mathbb{P}_{H_{1}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}>t_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}\left\\{\lambda^{-1}n\int(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))^{2}dH_{n,m}(\mathbf{x})\right.$
$\displaystyle\qquad+\left.(1-\lambda)^{-1}n\int(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x}))^{2}dH_{n,m}(\mathbf{x})>t_{1-\alpha}^{(2)}-a^{*}_{n,m}\right\\}$
$\displaystyle\rightarrow 1\text{ as }\min(n,m)\rightarrow\infty$ (26)

where $a^{*}_{n,m}=a_{n,m}-2b_{n,m}+2c_{n,m}-2d_{n,m}$ with

$\displaystyle a_{n,m}$
$\displaystyle=(n+m)\int(D_{F}(\mathbf{x})-D_{G}(\mathbf{x}))^{2}dH_{n,,m}(\mathbf{x})$
$\displaystyle b_{n,m}$
$\displaystyle=(n+m)\int(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x}))dH_{n,m}(\mathbf{x})$
$\displaystyle c_{n,m}$
$\displaystyle=(n+m)\int(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))(D_{F}(\mathbf{x})-D_{G}(\mathbf{x}))dH_{n,m}(\mathbf{x})$
$\displaystyle d_{n,m}$
$\displaystyle=(n+m)\int(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x}))(D_{F}(\mathbf{x})-D_{G}(\mathbf{x}))dH_{n,m}(\mathbf{x})$

The last implication follows from the following facts: (i)
$\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\\}$ and
$\sqrt{m}(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x}))\\}$ are uniformly
bounded in probability in view of Proposition 4, (ii) $t_{1-\alpha}^{(1)}$ is
positive finite, (iii)
$\lambda=\lim\limits_{\min(n,m)\rightarrow\infty}\frac{n}{n+m}=\lambda\in(0,1)$,
(iv) $H_{n,m}\rightarrow H$ almost surely, due to Glivenko-Cantalli’s theorem,
(iv) $a_{n,m}\rightarrow\infty$ $b_{n,m},c_{n,m}$ and $d_{n,m}$ are finite
(due to (i)-(iv)). Thus, with an application of Prokhorov’s theorem, we get
$\mathbb{P}_{H_{1}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}>t_{1-\alpha}^{(2)}\right\\}\rightarrow
1$ as $\min(n,m)\rightarrow\infty$.

Let $F_{n}$ and $G_{m}$ be any distribution functions that satisfies
$\sqrt{n+m}d_{K}(D_{F_{n}},D_{G_{m}})\geq\Delta_{n,m}$. By triangle
inequality, we have

$d_{K}(D_{F_{n}},D_{G_{m}})\leq
d_{K}(D_{F_{n}},D_{\mathcal{X}})+d_{K}(D_{\mathcal{X}},D_{\mathcal{Y}})+d_{K}(D_{\mathcal{Y}},D_{D_{m}})$
(27)

which implies
$T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}\geq\Delta_{n,m}-\sqrt{n+m}d_{K}(D_{\mathcal{X}},D_{F_{n}})-\sqrt{n+m}d_{K}(D_{\mathcal{Y}},D_{G_{m}})$.
Therefore,

$\displaystyle\mathbb{P}_{F_{n},G_{m}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}>t_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n},G_{m}}\left\\{\Delta_{n,m}-\sqrt{n+m}d_{K}(D_{\mathcal{X}},D_{F_{n}})-\sqrt{n+m}d_{K}(D_{\mathcal{Y}},D_{G_{m}})>t_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n},G_{m}}\left\\{\sqrt{n+m}d_{K}(D_{\mathcal{X}},D_{F_{n}})+\sqrt{n+m}d_{K}(D_{\mathcal{Y}},D_{G_{m}})\leq\Delta_{n,m}-t_{1-\alpha}^{(1)}\right\\}$
$\displaystyle\rightarrow 1\text{ as }\min(n,m)\rightarrow\infty.$ (28)

The last implication follows from the following facts: (i)
$\\{\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F}(\mathbf{x}))\\}$ and
$\\{\sqrt{m}(D_{\mathcal{Y}}(\mathbf{x})-D_{G}(\mathbf{x}))\\}$ are uniformly
bounded in probability under $F_{n}$ and $G_{m}$ respectively in the view of
Proposition 4, (ii) $\lambda\in(0,1)$ and $t_{1-\alpha}^{(1)}$ is positive
finite. (iii) $\Delta_{n,m}\rightarrow\infty$. An application of Prokhorov’s
theorem, we get
$\mathbb{P}_{F_{n},G_{m}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}>t_{1-\alpha}^{(1)}\right\\}\rightarrow
1$ as $\min(n,m)\rightarrow\infty$.

Moreover, by the inequality 27,

$\displaystyle\mathbb{P}_{F_{n},G_{m}}\left\\{T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}>s_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n},G_{m}}\left\\{\int\Delta_{n,m}^{2}dH_{n,m}(\mathbf{x})-(n+m)\int
d_{K}^{2}(D_{\mathcal{X}},D_{F_{n}})dH_{n,m}(\mathbf{x})\right.$
$\displaystyle\qquad\left.-(n+m)\int
d_{K}^{2}(D_{\mathcal{Y}},D_{G_{m}})dH_{n,m}(\mathbf{x})>t_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\geq\mathbb{P}_{F_{n},G_{m}}\left\\{(n+m)\int
d_{K}^{2}(D_{\mathcal{X}},D_{F_{n}})dH_{n,m}(\mathbf{x})\right.$
$\displaystyle\qquad\left.+(n+m)\int
d_{K}^{2}(D_{\mathcal{Y}},D_{G_{m}})dH_{n,m}(\mathbf{x})\leq\int\Delta_{n,m}^{2}d\mathbf{x}-t_{1-\alpha}^{(2)}\right\\}$
$\displaystyle\rightarrow 1\text{ as }\min(n,m)\rightarrow\infty.$ (29)

The last implication follows from the same facts that are used in Equation
(Proof of Theorem 7.). ∎

###### Proof of Theorem 8.

Observe that the log-likelihood ratio for testing $H_{0}:F=F_{0}$ against
$H_{n}$ described in (16),

$\displaystyle\mathcal{L}_{n}$
$\displaystyle=\sum_{i=1}^{n}\log\frac{(1-\gamma/\sqrt{n})f_{0}(\mathbf{x}_{i})+(\gamma/\sqrt{n})h(\mathbf{x}_{i})}{f_{0}(\mathbf{x}_{i})}$
$\displaystyle=\sum_{i=1}^{n}\log\left\\{1+(\gamma/\sqrt{n})\left[\frac{h(\mathbf{x}_{i})}{f_{0}(\mathbf{x}_{i})}-1\right]\right\\}$
$\displaystyle=\frac{\gamma}{\sqrt{n}}\sum_{i=1}^{n}\left\\{\frac{h(\mathbf{x}_{i})}{f_{0}(\mathbf{x}_{i}}-1\right\\}-\frac{\gamma^{2}}{2n}\sum_{i=1}^{n}\left\\{\frac{h(\mathbf{x}_{i})}{f_{0}(\mathbf{x}_{i}}-1\right\\}^{2}+\mathcal{R}_{n}$
$\displaystyle=\frac{\gamma}{\sqrt{n}}\sum_{i=1}^{n}\mathcal{K}_{i}-\frac{\gamma^{2}}{2n}\sum_{i=1}^{n}\mathcal{K}_{i}^{2}+\mathcal{R}_{n}$
(30)

where $\mathcal{K}_{i}=\frac{h({\mathbf{x}}_{i})}{f_{0}({\mathbf{x}}_{i})}-1$.
Since
$\sigma^{2}=\mathbb{E}_{F_{0}}\left\\{\frac{h(\mathbf{x})}{f_{0}(\mathbf{x})}-1\right\\}^{2}$
is finite, $\mathcal{R}_{n}\xrightarrow{\mathbb{P}}0$ as $n\rightarrow\infty$.
Contiguity of the sequence $H_{n}$ directly follows from Dhar et al. (2014)
(see Theorem 6.1) since the first term of Equation (Proof of Theorem 8.) is
asymptotically normal with mean zero and variance $\gamma^{2}\sigma^{2}$ due
to central limit theorem and second term converges in probability to
$\gamma^{2}\sigma^{2}/2$ due to weak law of large numbers. Therefore, by
Slutsky’s theorem, $\mathcal{L}_{n}$ is asymptotically normal with mean
$-\gamma^{2}\sigma^{2}/2$ and variance $\gamma^{2}\sigma^{2}$. An application
of Le Cam’s first lemma, we can deduce the first part of the theorem.

Now, to apply Le Cam’s third lemma, we need to calculate the covariance
between $D_{\mathcal{X}}-D_{F_{0}}$ and $\mathcal{L}_{n}$ under null.

Let $\mathcal{T}(F)$ be a functional defined for all distributions in a stable
class, then the influence function of $\mathcal{T}$ at $F$ is defined as
$\text{IF}(\mathbf{z};\mathcal{T}(F))=\lim\limits_{\epsilon\rightarrow
0^{+}}\frac{\mathcal{T}((1-\epsilon)F+\epsilon\delta_{\mathbf{z}})-\mathcal{T}(F)}{\epsilon}$
where $\delta_{\mathbf{z}}$ is the point mass probability measure at
$\mathbf{z}\in\mathbb{R}^{d}$ and $\epsilon\in[0,1]$. For any
$\mathbf{z}\in\mathbb{R}^{d}$, partition the set closed half-space
$\mathcal{H}$ as $\mathcal{H}_{\mathbf{z}}=\\{H\in\mathcal{H}:\mathbf{z}\in
H\\}$ and
$\mathcal{H}_{\overline{\mathbf{z}}}=\\{H\in\mathcal{H}:\mathbf{z}\notin
H\\}$. Thus, corresponding depths are defined as
$D_{F}^{\mathbf{z}}(\mathbf{x})=\inf_{\mathcal{H}_{\mathbf{z}}}F(H)$ and
$D_{F}^{\overline{\mathbf{z}}}(\mathbf{x})=\inf_{\mathcal{H}_{\overline{\mathbf{z}}}}F(H)$
respectively. Then the influence function of half-space depth becomes
$\text{IF}(\mathbf{z};D_{F}(\mathbf{x}))=-D_{F}(\mathbf{x})\textbf{1}\\{D_{F}^{(\overline{\mathbf{z}})}(\mathbf{x})<D_{F}^{(\mathbf{z})}(\mathbf{x})\\}+(1-D_{F}(\mathbf{x}))\textbf{1}\\{D_{F}^{(\overline{\mathbf{z}})}(\mathbf{x})\geq
D_{F}^{(\mathbf{z})}(\mathbf{x})\\}$ (Romanazzi, 2001), where
$\textbf{1}\\{a\in A\\}$ takes value $1$ if $a\in A$ and zero otherwise. If
$\mathbf{z}=\mathbf{x}$, then, $\mathcal{H}_{\mathbf{z}}=\mathcal{H}$ and
$\mathcal{H}_{\overline{\mathbf{z}}}$ becomes null-set. Hence
$\text{IF}(\mathbf{x};D_{F}(\mathbf{x}))=1-D_{F}(\mathbf{x})$. Note that, IF
is bounded, and is a step function. For all $\mathbf{x}$ belonging to optimal
half-space $\text{IF}(\mathbf{z};D_{F_{0}}(\mathbf{x}))$ is constant and equal
to $1-D_{F}(\mathbf{x})$ and for all $\mathbf{z}$ belonging to non-optimal
half-spaces, $\text{IF}(\mathbf{z};D_{F_{0}}(\mathbf{x}))$ is constant and
equal to $D_{F}(\mathbf{x})$. Since by the assumption optimal half-space depth
associated to $\mathbf{x}$ is unique, then for suitable regularity conditions
on $F_{0}$ (Serfling, 2009) we can write the von-Mises expansion (Romanazzi,
2001) as

$D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}\text{IF}(\mathbf{X}_{i};D_{F_{0}}(\mathbf{x}))+\mathcal{E}_{n}$
(31)

where $\mathcal{E}_{n}$ is the remaining term. Due to central limit theorem,
it can be shown that $\sqrt{n}\mathcal{E}_{n}\xrightarrow[]{\mathbb{P}}0$ (see
Appendix 2 of Romanazzi (2001)). Therefore, the asymptotic covariance function
between $\sqrt{n}(D_{\mathcal{X}}(\mathbf{x})-D_{F_{0}}(\mathbf{x}))$ and
$\mathcal{L}_{n}$ is

$\displaystyle\text{cov}_{H_{0}}\left\\{\sqrt{n}(D_{\mathcal{X}}({\mathbf{x}})-D_{F_{0}}({\mathbf{x}})),\mathcal{L}_{n}\right\\}$
$\displaystyle=\frac{\gamma}{n}\text{cov}_{H_{0}}\left\\{\sum_{i=1}^{n}\text{IF}(\mathbf{X}_{i};D_{F_{0}}({\mathbf{x}})),\sum_{i=1}^{n}\mathcal{K}_{i}\right\\}-\frac{\gamma^{2}}{2n^{3/2}}\text{cov}_{H_{0}}\left\\{\sum_{i=1}^{n}\text{IF}({\mathbf{X}_{i}};D_{F_{0}}({\mathbf{x}})),\sum_{i=1}^{n}\mathcal{K}_{i}^{2}\right\\}$
$\displaystyle=-\gamma\text{cov}_{H_{0}}\left\\{D_{F_{0}}({\mathbf{x}}),\frac{h({\mathbf{x}})}{f_{0}({\mathbf{x}})}\right\\}+\frac{\gamma^{2}}{2\sqrt{n}}\text{cov}_{H_{0}}\left\\{D_{F_{0}}(\mathbf{x}),\left(\frac{h(\mathbf{x})}{f_{0}(\mathbf{x})}-1\right)^{2}\right\\}$
$\displaystyle=-\gamma\int D_{F_{0}}(\mathbf{x})h(\mathbf{x})d\mathbf{x}+o(1)$
(32)

The last implication follows from the fact that
$\text{cov}_{H_{0}}\left\\{D_{F_{0}}(\mathbf{x}),\left(\frac{h(\mathbf{x})}{f_{0}(\mathbf{x})}-1\right)^{2}\right\\}$
is finite that can be shown by applying Cauchy-Schwarz inequality with the
fact that $\text{Var}_{H_{0}}\\{D_{F_{0}(\mathbf{x})}\\}$ and
$\mathbb{E}_{H_{0}}\left\\{\frac{h(\mathbf{x})}{f_{0}(\mathbf{x})}-1\right\\}^{4}$
are finite. Thus, by Le Cam’s third lemma, under contiguous alternatives the
empirical Tukey’s half-space depth process i.e.
$\sqrt{n}(D_{\mathcal{X}}({\mathbf{x}})-D_{F_{0}}({\mathbf{x}}))$ converges to
$\mathcal{G}_{1}^{\prime}$ which is Gaussian process with mean
$-\gamma\mathbb{E}_{\mathbf{x}\sim h}\left\\{D_{F_{0}}(\mathbf{x})\right\\}$
and the covariance kernel $F_{0}(H[\mathbf{x}_{1}]\cap
H[\mathbf{x}_{2}])-F_{0}(H[\mathbf{x}_{1}])F_{0}(H[\mathbf{x}_{2}])$.
Moreover, $\sqrt{n}(D_{\mathcal{X}}({\mathbf{x}})-D_{F_{0}}({\mathbf{x}}))$
satisfies the tightness condition under contiguous alternatives since it is
tight under null. The tightness under null follows from the weak convergence
of the empirical Tukey’s depth process (see Proposition 4). Therefore, by Le
Cam’s third lemma, under the contiguous alternatives alternatives $H_{n}$, the
asymptotic power of the test based on $T_{\mathcal{X},F_{0}}^{\text{KS}}$ and
$T_{\mathcal{X},F_{0}}^{\text{CvM}}$ are
$\mathbb{P}_{\gamma}\left\\{\sup\limits_{\mathbf{x}}|\mathcal{G}_{1}^{\prime}(\mathbf{x})|>s_{1-\alpha}^{(1)}\right\\}$
and
$\mathbb{P}_{\gamma}\left\\{\int\limits|\mathcal{G}_{1}^{\prime}(\mathbf{x})|^{2}dF_{0}(\mathbf{x})>s_{1-\alpha}^{(2)}\right\\}$,
respectively. ∎

###### Proof of Theorem 8.

Observe that the likelihood ratio for testing $H_{0}:F=G$ against $H_{n,m}$
described in (17),

$\displaystyle\mathcal{L}_{n,m}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{m}\log\frac{f(\mathbf{x}_{i})\left\\{(1-\gamma/\sqrt{n+m})f(y_{j})+\gamma
h(\mathbf{y}_{j})/\sqrt{n+m}\right\\}}{f(\mathbf{x}_{i})f(\mathbf{y}_{i})}$
$\displaystyle=\sum_{j=1}^{m}\log\left\\{1+\frac{\gamma}{\sqrt{n+m}}\left(\frac{h(y_{j})}{f(y_{j})}-1\right)\right\\}$
$\displaystyle=\frac{\gamma}{\sqrt{n+m}}\sum_{j=1}^{m}\left\\{\frac{h(\mathbf{y}_{j})}{f(\mathbf{y}_{j})}-1\right\\}-\frac{\gamma^{2}}{2(n+m)}\sum_{j=1}^{m}\left\\{\frac{h(\mathbf{y}_{j})}{f(\mathbf{y}_{j})}-1\right\\}^{2}+\mathcal{R}_{n,m}$
$\displaystyle=\frac{\gamma}{\sqrt{n+m}}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime}-\frac{\gamma^{2}}{2(n+m)}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime
2}+\mathcal{R}_{n,m}$ (33)

where
$\mathcal{K}_{i}^{\prime}=\frac{h(\mathbf{y}_{j})}{f(\mathbf{y}_{j})}-1$.
Since
$\sigma^{2}=\mathbb{E}\left\\{\frac{h(\mathbf{y})}{f(\mathbf{y})}-1\right\\}^{2}$
is finite, $\mathcal{R}_{n,m}\xrightarrow[]{\mathbb{P}}0$ as
$\min(n,m)\rightarrow\infty$. Contiguity of the sequence $H_{n,m}$ directly
follows from Dhar et al. (2014) (see Theorem 6.2) since the first term of
Equation (Proof of Theorem 8.) is asymptotically normal with mean zero and
variance $\gamma^{2}\sigma^{2}(1-\lambda)$ due to central limit theorem and
second term converges in probability to $-\gamma^{2}\sigma^{2}(1-\lambda)/2$.
Therefore, by Slutsky’s theorem, $\mathcal{L}_{n,m}$ is asymptotically normal
with mean $-\gamma^{2}\sigma^{2}(1-\lambda)/2$ and variance
$\gamma^{2}\sigma^{2}(1-\lambda)$. An application of Le Cam’s first lemma, we
can deduce the first part of the theorem.

Now, to apply Le Cam’s third lemma, we need to calculate the covariance
between $D_{\mathcal{X}}-D_{\mathcal{Y}}$ and $\mathcal{L}_{n,m}$ under null.
As in the proof of Theorem 8, using von-Mises expansion described in Equation
(31) and the fact that $\mathcal{X}$ and $\mathcal{Y}$ are independent, we
have

$\displaystyle\text{cov}_{H_{0}}\left\\{\sqrt{n+m}(D_{\mathcal{X}}(\mathbf{u})-D_{\mathcal{Y}}(\mathbf{u})),\mathcal{L}_{n,m}\right\\}$
$\displaystyle=\gamma\text{cov}_{H_{0}}\left\\{\sqrt{n+m}(D_{\mathcal{X}}(\mathbf{u})-D_{\mathcal{Y}}(\mathbf{u})),\frac{1}{\sqrt{n+m}}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime}\right\\}$
$\displaystyle\qquad-\frac{\gamma^{2}}{2}\text{cov}_{H_{0}}\left\\{\sqrt{n+m}(D_{\mathcal{X}}(\mathbf{u})-D_{\mathcal{Y}}(\mathbf{u})),\frac{1}{(n+m)}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime
2}\right\\}:=\gamma\mathcal{A}_{1}-\frac{\gamma^{2}}{2}\mathcal{A}_{2}$ (34)

where

$\displaystyle\mathcal{A}_{1}$
$\displaystyle=\text{cov}_{H_{0}}\left\\{\sqrt{n+m}(D_{\mathcal{X}}(\mathbf{u})-D_{\mathcal{Y}}(\mathbf{u})),\frac{1}{\sqrt{n+m}}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime}\right\\}$
$\displaystyle=-\text{cov}_{H_{0}}\left\\{(1-\lambda)^{-1/2}\sqrt{m}(D_{\mathcal{Y}}(\mathbf{u})-D_{F}(\mathbf{u})),\frac{1}{\sqrt{n+m}}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime}\right\\}$
$\displaystyle=-(1-\lambda)^{-1/2}\text{cov}_{H_{0}}\left\\{\frac{1}{\sqrt{m}}\sum_{j=1}^{m}\text{IF}(\mathbf{y}_{j};D_{F}(\mathbf{u})),\frac{1}{\sqrt{n+m}}\sum_{j=1}^{m}\mathcal{K}_{j}^{\prime}\right\\}$
$\displaystyle=\sqrt{\frac{\lambda}{1-\lambda}}\text{cov}_{H_{0}}\left\\{D_{F}(\mathbf{u}),\frac{h(\mathbf{u})}{f(\mathbf{u})}\right\\}=\sqrt{\frac{\lambda}{1-\lambda}}\int
D_{F}(\mathbf{u})h(\mathbf{u})d\mathbf{u}$ (35)

and similar to the second term of Equation (Proof of Theorem 8.), under the
condition
$\mathbb{E}\left\\{\frac{h(\mathbf{x})}{f(\mathbf{x})}-1\right\\}^{4}<\infty$,
by applying Cauchy-Schwarz inequality, $\mathcal{A}_{2}=o(1)$ as
$\min(n,m)\rightarrow\infty$. Thus, the covariance between
$\sqrt{n+m}(D_{\mathcal{X}}(\mathbf{u})-D_{\mathcal{Y}}(\mathbf{u}))$ and
$\mathcal{L}_{n,m}$ is
$\gamma\sqrt{\lambda/(1-\lambda)}\mathbb{E}_{\mathbf{u}\sim
h}\left\\{D_{F}(\mathbf{u})\right\\}$.

Thus, by Le Cam’s third lemma, under contiguous alternatives the empirical
Tukey’s half-space depth process i.e.
$\sqrt{n+m}(D_{\mathcal{X}}({\mathbf{u}})-D_{\mathcal{Y}}({\mathbf{u}}))$
converges to $\mathcal{G}_{2}^{\prime}$ which is Gaussian process with mean
$\gamma\sqrt{\lambda/(1-\lambda)}\mathbb{E}_{\mathbf{u}\sim
h}\left\\{D_{F}(\mathbf{u})\right\\}$ and the covariance kernel
$\\{F(H[\mathbf{u}_{1}]\cap
H[\mathbf{u}_{2}])-F(H[\mathbf{u}_{1}])F_{0}(H[\mathbf{u}_{2}])\\}/\lambda(1-\lambda)$.
Moreover,
$\sqrt{n+m}(D_{\mathcal{X}}({\mathbf{u}})-D_{\mathcal{Y}}({\mathbf{u}}))$
satisfies the tightness condition under contiguous alternatives since it is
tight under null. The tightness under null follows from the weak convergence
of the empirical Tukey’s depth process (under the independence of
$\mathcal{X}$ and $\mathcal{Y}$, for $\lambda\in(0,1)$, see Proposition 4).
Therefore, by Le Cam’s third lemma, under the contiguous alternatives
alternatives $H_{n}$, the asymptotic power of the test based on
$T_{\mathcal{X},\mathcal{Y}}^{\text{KS}}$ and
$T_{\mathcal{X},\mathcal{Y}}^{\text{CvM}}$ are
$\mathbb{P}_{\gamma}\left\\{\sup\limits_{\mathbf{u}}|\mathcal{G}_{2}^{\prime}(\mathbf{u})|>t_{1-\alpha}^{(1)}\right\\}$
and
$\mathbb{P}_{\gamma}\left\\{\int\limits|\mathcal{G}_{2}^{\prime}(\mathbf{u})|^{2}dH_{n,m}(\mathbf{u})>t_{1-\alpha}^{(2)}\right\\}$,
respectively. ∎

## References

  * Anderson (1962) Anderson, T. W. (1962).  On the distribution of the two-sample cramer-von mises criterion.  The Annals of Mathematical Statistics, 1148–1159.
  * Anderson and Darling (1952) Anderson, T. W. and D. A. Darling (1952).  Asymptotic theory of certain” goodness of fit” criteria based on stochastic processes.  The annals of mathematical statistics, 193–212.
  * Azzalini and Valle (1996) Azzalini, A. and A. D. Valle (1996).  The multivariate skew-normal distribution.  Biometrika 83(4), 715–726.
  * Brown and Hettmansperger (1989) Brown, B. M. and T. P. Hettmansperger (1989).  An affine invariant bivariate version of the sign test.  Journal of the Royal Statistical Society: Series B (Methodological) 51(1), 117–125.
  * Chambers et al. (2018) Chambers, J. M., W. S. Cleveland, B. Kleiner, and P. A. Tukey (2018).  Graphical methods for data analysis.  Chapman and Hall/CRC.
  * Chen and Friedman (2017) Chen, H. and J. H. Friedman (2017).  A new graph-based two-sample test for multivariate and object data.  Journal of the American statistical association 112(517), 397–409.
  * Cuesta-Albertos and Nieto-Reyes (2008) Cuesta-Albertos, J. and A. Nieto-Reyes (2008).  The tukey and the random tukey depths characterize discrete distributions.  Journal of Multivariate Analysis 99(10), 2304 – 2311\.
  * Daniel (1990) Daniel, W. W. (1990).  Kolmogorov–smirnov one-sample test.  Applied nonparametric statistics 2.
  * Dhar et al. (2014) Dhar, S. S., B. Chakraborty, and P. Chaudhuri (2014).  Comparison of multivariate distributions using quantile–quantile plots and related tests.  Bernoulli 20(3), 1484–1506.
  * Donoho and Gasko (1992) Donoho, D. L. and M. Gasko (1992).  Breakdown properties of location estimates based on halfspace depth and projected outlyingness.  The Annals of Statistics, 1803–1827.
  * Doornik and Hansen (2008) Doornik, J. A. and H. Hansen (2008).  An omnibus test for univariate and multivariate normality.  Oxford bulletin of economics and statistics 70, 927–939.
  * Efron (1982) Efron, B. (1982).  The jackknife, the bootstrap and other resampling plans.  SIAM.
  * Friedman and Rafsky (1981) Friedman, J. H. and L. C. Rafsky (1981).  Graphics for the multivariate two-sample problem.  Journal of the American Statistical Association 76(374), 277–287.
  * Gnanadesikan (2011) Gnanadesikan, R. (2011).  Methods for statistical data analysis of multivariate observations.  John Wiley & Sons.
  * Gower (1974) Gower, J. (1974).  Algorithm as 78: The mediancentre.  Journal of the Royal Statistical Society. Series C (Applied Statistics) 23(3), 466–470.
  * Hassairi and Regaieg (2008) Hassairi, A. and O. Regaieg (2008).  On the tukey depth of a continuous probability distribution.  Statistics & probability letters 78(15), 2308–2313.
  * Henze and Zirkler (1990) Henze, N. and B. Zirkler (1990).  A class of invariant consistent tests for multivariate normality.  Communications in statistics-Theory and Methods 19(10), 3595–3617.
  * Hodges (1955) Hodges, J. L. (1955).  A bivariate sign test.  The Annals of Mathematical Statistics 26(3), 523–527.
  * Kong and Zuo (2010) Kong, L. and Y. Zuo (2010).  Smooth depth contours characterize the underlying distribution.  Journal of Multivariate Analysis 101(9), 2222 – 2226\.
  * Koshevoy (2002) Koshevoy, G. A. (2002).  The tukey depth characterizes the atomic measure.  Journal of Multivariate Analysis 83(2), 360–364.
  * Koshevoy (2003) Koshevoy, G. A. (2003).  Lift-zonoid and multivariate depths.  In Developments in Robust Statistics, pp. 194–202. Springer.
  * Koziol (1982) Koziol, J. A. (1982).  A class of invariant procedures for assessing multivariate normality.  Biometrika 69(2), 423–427.
  * Liu (1990) Liu, R. (1990).  Liu, on a notion of data depth based on random simplices.  Ann. Statist 18, 405–414.
  * Liu et al. (1999) Liu, R. Y., J. M. Parelius, and K. Singh (1999).  Multivariate analysis by data depth: descriptive statistics, graphics and inference, (with discussion and a rejoinder by liu and singh).  Ann. Statist. 27(3), 783–858.
  * Lorenz (1905) Lorenz, M. O. (1905).  Methods of measuring the concentration of wealth.  Publications of the American statistical association 9(70), 209–219.
  * Mahalanobis (1936) Mahalanobis, P. C. (1936).  On the generalized distance in statistics.  National Institute of Science of India.
  * Marden (1998) Marden, J. I. (1998).  Bivariate qq-plots and spider web plots.  Statistica Sinica, 813–826.
  * Marden (2004) Marden, J. I. (2004).  Positions and qq plots.  Statistical Science, 606–614.
  * Massé (2004) Massé, J.-C. (2004).  Asymptotics for the tukey depth process, with an application to a multivariate trimmed mean.  Bernoulli 10(3), 397–419.
  * Michael (1983) Michael, J. R. (1983).  The stabilized probability plot.  Biometrika 70(1), 11–17.
  * Mosler (2013) Mosler, K. (2013).  Depth statistics.  Robustness and complex data structures, 17–34.
  * Newland et al. (2016) Newland, R. P., S. H. Parade, S. Dickstein, and R. Seifer (2016).  Goodness of fit between prenatal maternal sleep and infant sleep: Associations with maternal depression and attachment security.  Infant behavior and development 44, 179–188.
  * Oja (1983) Oja, H. (1983).  Descriptive statistics for multivariate distributions.  Statistics & Probability Letters 1(6), 327–332.
  * Paraschiv-Ionescu et al. (2013) Paraschiv-Ionescu, A., E. Buchser, and K. Aminian (2013).  Unraveling dynamics of human physical activity patterns in chronic pain conditions.  Scientific reports 3(1), 1–10.
  * Paulson et al. (1987) Paulson, A., P. Roohan, and P. Sullo (1987).  Some empirical distribution function tests for multivariate normality.  Journal of Statistical Computation and Simulation 28(1), 15–30.
  * Pollard (1984) Pollard, D. (1984).  Convergence of stochastic processes.  Springer Science & Business Media.
  * Romanazzi (2001) Romanazzi, M. (2001).  Influence function of halfspace depth.  Journal of Multivariate Analysis 77(1), 138 – 161.
  * Rousseeuw and Ruts (1996) Rousseeuw, P. J. and I. Ruts (1996).  Algorithm as 307: Bivariate location depth.  Journal of the Royal Statistical Society. Series C (Applied Statistics) 45(4), 516–526.
  * Rousseeuw and Struyf (1998) Rousseeuw, P. J. and A. Struyf (1998).  Computing location depth and regression depth in higher dimensions.  Statistics and Computing 8(3), 193–203.
  * Royston (1992) Royston, P. (1992).  Approximating the shapiro-wilk w-test for non-normality.  Statistics and computing 2(3), 117–119.
  * Ruts and Rousseeuw (1996) Ruts, I. and P. J. Rousseeuw (1996).  Computing depth contours of bivariate point clouds.  Computational Statistics & Data Analysis 23(1), 153 – 168.  Classification.
  * Serfling (2009) Serfling, R. J. (2009).  Approximation theorems of mathematical statistics, Volume 162.  John Wiley & Sons.
  * Sidak et al. (1999) Sidak, Z., P. K. Sen, and J. Hajek (1999).  Theory of rank tests.  Academic press.
  * Tibshirani and Efron (1993) Tibshirani, R. J. and B. Efron (1993).  An introduction to the bootstrap.  Monographs on statistics and applied probability 57, 1–436.
  * Tukey (1975) Tukey, J. W. (1975).  Mathematics and the picturing of data.  In Proceedings of the international congress of mathematicians, Volume 2, pp. 523–531.
  * Voinov et al. (2016) Voinov, V., N. Pya, R. Makarov, and Y. Voinov (2016).  New invariant and consistent chi-squared type goodness-of-fit tests for multivariate normality and a related comparative simulation study.  Communications in Statistics-Theory and Methods 45(11), 3249–3263.
  * Webster (1977) Webster, R. (1977).  Spectral analysis of gilgai soil.  Soil Research 15(3), 191–204.
  * Wilk and Gnanadesikan (1968) Wilk, M. B. and R. Gnanadesikan (1968).  Probability plotting methods for the analysis for the analysis of data.  Biometrika 55(1), 1–17.
  * Zuo and Serfling (2000) Zuo, Y. and R. Serfling (2000).  General notions of statistical depth function.  Annals of statistics, 461–482.