# Living on the Edge:
An Unified Approach to Antithetic Sampling

Roberto Casarin22footnotemark: 2 Radu V. Craiu33footnotemark: 3
Lorenzo Frattarolo44footnotemark: 4 Christian P. Robert 55footnotemark: 5

22footnotemark: 2 University Ca’ Foscari of Venice
33footnotemark: 3 University of Toronto
44footnotemark: 4 European Commission Joint Research Centre
55footnotemark: 5 CEREMADE, Université Paris-Dauphine and University of
Warwick

###### Abstract

We identify recurrent ingredients in the antithetic sampling literature
leading to a unified sampling framework. We introduce a new class of
antithetic schemes that includes the most used antithetic proposals. This
perspective enables the derivation of new properties of the sampling schemes:
i) optimality in the Kullback-Leibler sense; ii) closed-form multivariate
Kendall’s $\tau$ and Spearman’s $\rho$; iii)ranking in concordance order and
iv) a central limit theorem that characterizes stochastic behavior of Monte
Carlo estimators when the sample size tends to infinity. Finally, we provide
applications to Monte Carlo integration and Markov Chain Monte Carlo Bayesian
estimation.

Keywords: Antithetic variables, Countermonotonicity, Monte Carlo, Negative
dependence, Variance reduction

##  1 Introduction

The Monte Carlo method is at the core of model-based scientific exploration.
In its simplest form it relies on approximating an integral $\mathfrak{I}=\int
f(\mathbf{x})\pi(d\mathbf{x})$ with
$\hat{\mathfrak{I}}_{d}=\frac{1}{d}\sum_{i=1}^{d}f(\mathbf{X}_{i})$ when $\pi$
is a probability measure, $f:\mathbb{R}^{p}\mapsto\mathbb{R}$ is a integrable
function with respect to $\pi$, $d$ is the Monte Carlo sample size and
$\mathbf{X}_{1},\ldots,\mathbf{X}_{d}$ are independent, identically
distributed (henceforth, iid) samples from $\pi$.

In modern computational problems, sampling from the distribution $\pi$ may be
expensive, in terms of either computational effort or time, so techniques
needed to reduce the Monte Carlo sample size $d$, while maintaining the
desired precision in estimation, are essential. A relevant class is
represented by the variance reduction techniques that use statistical
properties induced by the sampling design to reduce
$Var(\hat{\mathfrak{I}}_{d})$. For instance, in the case $p=1$, if the
independence condition between samples $X_{1},\ldots,X_{d}$ is dropped then

$\hbox{Var}(\hat{\mathfrak{I}}_{d})=\dfrac{1}{d}\sum^{d}_{i=1}\hbox{Var}\left(f\left(X_{i}\right)\right)+\dfrac{1}{d^{2}}\sum_{i\neq
j}\hbox{Cov}\left(f\left(X_{i}\right),f\left(X_{j}\right)\right),$ (1)

and the variance is reduced, compared to independent sampling, if the average
covariance is negative.

Antithetic sampling designs aim at minimizing the covariances between samples
while preserving their marginal distribution. A historical perspective on the
strategies for antithetic sampling (e.g., Hammersley and Mauldon,, 1956;
Hammersley and Morton,, 1956) allows us to better understand the rationale
behind various constructions and to establish useful relationships with the
results available from related fields, such as stochastic orders (e.g., Barlow
and Proschan,, 1975), optimal transport (e.g., Gaffke and Rüschendorf,, 1981),
Fréchet classes (e.g. Whitt,, 1976), and group transformation (e.g.
Andréasson,, 1972). Our historical review identifies some key recurrent
ingredients used to propose a unified framework for antithetic sampling. We
introduce a new class of antithetic constructions that also includes some of
the, to our knowledge, most used antithetic proposals, which are reviewed
later in this section. Moreover, this new perspective enables the derivation
of new properties of the sampling schemes for $p$ stochastically independent
replications ($p\geq 1$): i) optimality in the Kullback-Leibler sense; ii)
closed-form multivariate Kendall’s $\tau$ and Spearman’s $\rho$; iii) ranking
in concordance order, and iv) a central limit theorem that characterizes
stochastic behavior when $d$ tends to infinity.

The pairwise antithetic coupling introduced by Hammersley and Morton, (1956)
achieves variance reduction by generating $d/2$ (we assume $d$ is even in (1))
iid pairs of negatively correlated random variables. The joint distribution of
any pair, say $(X_{1},X_{2})$, achieves the lower Fréchet bound (see Fréchet,,
1935). This is achieved by sampling $X_{1}\sim\pi$ and sampling the second
variable $X_{2}$ via

$X_{2}=F_{\pi}^{-1}\left(1-F_{\pi}\left(X_{1}\right)\right),$

where $F_{\pi}$ is the cumulative distribution function (CDF) for the
distribution $\pi$. This procedure minimizes the correlation for any monotonic
$f$. The explanation of this reduction, as correctly pointed out by Whitt,
(1976), can be found in Hoeffding, (1940) and Fréchet, (1951) and it is due to
the rearrangement inequality (see chapter X of Hardy et al.,, 1934).
Interestingly, the construction cannot be unambiguously extended beyond pairs
because the lower Fréchet bound of all $k$-variate distributions is itself a
distribution only when $k=2$.

In the following, we establish a relationship between variance reduction and
stochastic orders. This relationship allows generalizing variance reduction to
higher dimensions $d>2$.

######  Definition 1.1 (Concordance Order (Joe,, 1990)).

Let $\mathbf{X}$ and $\mathbf{Y}$ be random vectors with CDFs $F$ and $G$ and
survival function $\bar{F}$ and $\bar{G}$, respectively. Then $\mathbf{Y}$ is
more concordant than $\mathbf{X}$ (written $\mathbf{X}\prec_{C}\mathbf{Y}$) if

$\begin{array}[]{ccc}F\leq G&\hbox{ and }&\bar{F}\leq\bar{G}\end{array}$ (2)

The pairwise antithetic construction being based on the lower Fréchet bound is
minimal in concordance order. If the univariate marginals are the same in the
two-dimensional case, then the concordance order is equivalent to the
correlation order, which is a stochastic order induced by the covariance among
coordinate-wise monotonic functions. The correlation order was introduced in
Barlow and Proschan, (1975) an studied further in Dhaene and Goovaerts,
(1996). In this paper, we follow the multivariate definition and notation from
Lu and Yi, (2004).

######  Definition 1.2 (Correlation Order (Lu and Yi,, 2004)).

Denote with $\mathcal{G}^{d}$ the set of real-valued bounded coordinate-wise
non-decreasing Borel functions on $\mathbb{R}^{d}$, $d\in\mathbb{N}$ and with
$\mathcal{D}=\left\\{1,\ldots,d\right\\}$ the set of coordinate indexes. The
d-variate random vector $\mathbf{X}$ is less correlated than the d-variate
random vector $\mathbf{Y}$, written $\mathbf{X}\prec_{corr}\mathbf{Y}$, if
they have the same univariate marginal CDFs $F_{i}$, $i=1,\ldots,d$ and if for
every pair of disjoint finite subsets
$\mathcal{I},\mathcal{J}\subset\mathcal{D}$ and all functions
$f\in\mathcal{G}^{\left|\mathcal{I}\right|},g\in\mathcal{G}^{\left|\mathcal{J}\right|}$,
for which the covariances exist:

$\mathbb{C}\mathrm{ov}\left(f\left(X_{\mathcal{I}}\right),g\left(X_{\mathcal{J}}\right)\right)\leq\mathbb{C}\mathrm{ov}\left(f\left(Y_{\mathcal{I}}\right),g\left(Y_{\mathcal{J}}\right)\right)$

Since concordance and correlation orders are equivalent in dimension two (see,
for example, Dhaene and Goovaerts,, 1996), antithetic coupling minimizes the
covariance among monotone functions.

If we drop the monotonicity assumption about $f$ in (1), the discussion and
derivation of lower bounds for the variance are more complex and less general.
If $p=1$ and $d=2$, for non-monotonic, bounded $f$, Hammersley and Mauldon,
(1956) prove that the lower bound of the variance can be attained only by a
multivariate transformation of a single standard uniform random variable
which, almost surely, is coordinate-wise monotonic. The proof relies on two
main ingredients. First, the monotonic transformation introduces an
approximate representation of the class of bivariate distributions with
uniform marginals. The candidate member of the class is approximated by
partitioning the unit square in sub-squares of side $1/n$. The approximation
is a doubly stochastic matrix in which each element corresponds to a sub-
square and has a value equal to the mass assigned to the corresponding sub-
square. This construction relies on the bijective rearrangement (Puccetti and
Wang,, 2015) also known as measure-preserving transformation (e.g. Brown,,
1966; Vitale,, 1990). In the interpretation of Vitale, (1990), for every
random vector $\left(U_{1},U_{2}\right)$ on the unit square, with standard
uniform marginals, there is a sequence of bijective maps $f_{n}$ such that
$\left(U_{1},f_{n}\left(U_{1}\right)\right)$ weakly converges to
$\left(U_{1},U_{2}\right)$. Bijective rearrangement and the induced stochastic
dependence (Durante and Sanchez,, 2012) are relevant to our discussion of the
antithetic constructions for $d>2$.

The second ingredient is the Birkhoff-Von Neumman’s decomposition (Birkhoff,,
1946; Von Neumann,, 1953) of doubly-stochastic matrices in which they are
represented as convex combinations of permutation matrices. Handscomb, (1958)
extends the latter result to characterize the extremal points of multi-
stochastic arrays as higher-dimension permutation arrays and to provide a
generalization for $d>2$ of the results in Hammersley and Mauldon, (1956).

Those early results based on discretization give sufficient conditions to
characterize the transformations needed to obtain a minimal variance.
Unfortunately, the characterizations are not constructive and do not provide
feasible random sampling algorithms. Moreover, the existence of the optimal
transformation minimizing the variance is not guaranteed.

This led earlier researchers to propose feasible, yet sub-optimal sampling
solutions, including Hammersley and Morton’s (Hammersley and Morton,, 1956)
proposal for $d>2$. Andréasson, (1972), Andréasson and Dahlquist, (1972) and
Roach and Wright, (1977) follow a group theoretic approach. In particular,
Roach and Wright, (1977), build on Andréasson, (1972), Andréasson and
Dahlquist, (1972) and Tukey, (1957), draws a parallel with systematic
sampling. Its solutions for $d=2$, in the case of non-monotonic $f$s, relies
on discretization, optimal transport and a branch and bound algorithm. The
group theoretic approach was also used in Fishman and Huang, (1983) to obtain
a reinterpretation of the original Hammersley and Morton, (1956) proposal for
$d>2$. Their construction was named rotation sampling and is described next.

######  Example 1.1 (Rotation sampling).

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle U\sim\mathcal{U}[0,1]$
$\displaystyle U_{l}$ $\displaystyle=$
$\displaystyle\left(\dfrac{l-1}{d}+U\right)\mathrm{mod}\,1,\quad
l\in\left\\{2,\ldots,d\right\\}$ (3)

This proposal is a particular case of our stochastic representation.

The extension to unbounded functions of the theorems in Hammersley and
Mauldon, (1956) and Handscomb, (1958) can be found in Wilson, (1979) for
$d\geq 2$ and $p=1$ and in Wilson, (1983) for $d\geq 2$, $p>=1$. The latter
paper combines discretization and bijective rearrangement with the optimal
transport assignment problem to prove the results. Bijective rearrangement and
Monge-Kantorowitch transportation problem are used in Gaffke and Rüschendorf,
(1981) to obtain minimum variance constructions for $f$ equal to the identity
function. The authors are the first to realize that the Hammersley and Morton,
(1956) bivariate antithetic vector has an almost sure constant sum, which is
one of the main ingredients of our unified approach.

The relationship between constant sum and variance reduction is trivial.
Random vectors of dimension $d\geq 2$ with constant sum achieve the smallest
variance for the sum of their components. Then they minimize the variance in
the case $p=1$ and $f$ equal to the identity function. Beyond that, a recent
stream of papers proves that the constant sum vectors are among the minimal
vectors with respect to the concordance order. In particular, one possible
generalization of the constant sum constraint is the following one.

######  Definition 1.3 ($l$-countermonotonic).

A $d$-dimensional random vector $\mathbf{U}$ with uniform marginals, is said
to be $l$-countermonotonic ($l$-CTM), if there exist some
$\mathcal{L}\subseteq\mathcal{D}$ with $\left|\mathcal{L}\right|=l$, a family
$\left\\{g_{l}\right\\}_{l\in\mathcal{L}}$ of strictly increasing continuous
functions $\left[0,1\right]\mapsto\mathbb{R}$ and some $k\in\mathbb{R}$ such
that :

$\sum_{l\in\mathcal{L}}g_{l}\left(U_{l}\right)=k\hbox{ a.s.}$ (4)

Theorem 2 and Proposition 1 in Lee and Ahn, (2014) show that the antithetic
vector is the only element of the $2$-CTM class and it is minimal in the
concordance and correlation orders. The equivalence among orders does not hold
in the multivariate case, but several relationships have been identified. For
instance, Lu and Yi, (2004) have showed that
$\mathbf{X}\prec_{corr}\mathbf{Y}$ implies $\mathbf{X}\prec_{C}\mathbf{Y}$.

In addition, conditions for achieving minimality in the concordance order were
linked to $l$-CTM random vectors. If $\mathbf{V}$ is $d$-CTM or $(d-1)$-CTM
and $\mathbf{U}\prec_{C}\mathbf{V}$ then $\mathbf{U}$ has the same
distribution as $\mathbf{V}$ (Lee et al.,, 2017). The same result has been
shown to hold for $\prec_{corr}$ order (see Ahn and Fuchs,, 2020, remark 3.1).
Since $d$-CTM and $(d-1)$-CTM are among the minimal elements in the
correlation order, they are promising candidates for constructing variance
reduction techniques in any dimension $d$. In this paper, we study some of the
existing sampling methods and propose new constructions for $d$-CTM vectors of
Uniform(0,1) random variables with a.s. constant sum, that is
$g_{l}\left(U_{l}\right)=U_{l}$, $l\in\mathcal{D}$ in Eq (4). This subclass is
known in the literature as strict $d$-CTM (Lee and Ahn,, 2014). Gaffke and
Rüschendorf, (1981) recognize that Hammersley and Morton, (1956) is strict
$2$-CTM and provide the first strict $3$-CTM construction given below.

######  Example 1.2 (Gaffke and Rüschendorf, (1981) strict $3$-CTM).

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle U,\qquad U\sim
U\left[0,1\right]$ $\displaystyle U_{2}$ $\displaystyle=$ $\displaystyle
U+\frac{1}{2}\mathbb{I}_{\left[0,1/2\right]}\left(U\right)-\frac{1}{2}\mathbb{I}_{\left[1/2,1\right]}\left(U\right)$
(5) $\displaystyle U_{3}$ $\displaystyle=$
$\displaystyle-2U+\mathbb{I}_{\left[0,1/2\right]}\left(U\right)+2\mathbb{I}_{\left[1/2,1\right]}\left(U\right)$

For the case $d>3$, the authors propose to generate a sequence of independent
random vectors using their representation in (1.2) and the bivariate
antithetic vector of Hammersley and Morton, (1956).

######  Example 1.3 (Gaffke and Rüschendorf, (1981) strict $d$-CTM).

Let $V_{i}$, $i=1,\ldots,\left\lfloor(d-2)/2\right\rfloor+1$ independent
random variables, and

$\displaystyle U_{2i-1}=V_{i},\quad U_{2i}=1-V_{i},\quad
i=1,\ldots,\left\lfloor(d-2)/2\right\rfloor$

with

$\displaystyle U_{d-1}=V_{\left\lfloor(d-2)/2\right\rfloor+1},\quad
U_{d}=1-V_{\left\lfloor(d-2)/2\right\rfloor}$ (6)

if $d$ even, and

$\displaystyle U_{d-2}$ $\displaystyle=$ $\displaystyle
V_{\left\lfloor(d-2)/2\right\rfloor+1}$ $\displaystyle U_{d-1}$
$\displaystyle=$ $\displaystyle
V_{\left\lfloor(d-2)/2\right\rfloor+1}+1/2\mathbb{I}_{\left[0,1/2\right]}\left(V_{\left\lfloor(d-2)/2\right\rfloor+1}\right)-1/2\mathbb{I}_{\left[1/2,1\right]}\left(V_{\left\lfloor(d-2)/2\right\rfloor+1}\right)$
$\displaystyle U_{d}$ $\displaystyle=$
$\displaystyle-2U+\mathbb{I}_{\left[0,1/2\right]}\left(V_{\left\lfloor(d-2)/2\right\rfloor+1}\right)+2\mathbb{I}_{\left[1/2,1\right]}\left(V_{\left\lfloor(d-2)/2\right\rfloor+1}\right)$
(7)

if $d$ is odd.

Almost contemporaneously, Arvidsen and Johnsson, (1982) put forward the
apparently different proposal given in the following.

######  Example 1.4 (Arvidsen and Johnsson, (1982) strict $d$-CTM).

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle U\qquad
U\sim\mathcal{U}[0,1]$ $\displaystyle U_{i}$ $\displaystyle=$
$\displaystyle\left(2^{i-2}U_{1}+1/2\right)\mathrm{mod}\,1,\quad
i=\left\\{2,\ldots,d-1\right\\}$ (8) $\displaystyle U_{d}$ $\displaystyle=$
$\displaystyle 1-\left(2^{d-2}U_{1}\right)\mathrm{mod}\,1.$

We will show that for $d=3$, the two proposals in Examples 1.2 and 1.4
coincide and are special cases of our general stochastic representation. Both
constructions yield vectors with a constant sum, but Gaffke and Rüschendorf,
(1981) proposal’s use of independent random variates for $d>3$ made us wonder
about its efficiency, especially since the results in Hammersley and Mauldon,
(1956) and Handscomb, (1958) suggest that combinations of independent vectors
could be sub-optimal. We compare different strict $d$-CTM constructions using
the concordance order. According to Ahn and Fuchs, (2020), all strict $d$-CTM
have minimal multivariate Kendall’s $\tau$, but they can have different
multivariate Spearman’s $\rho$ values. For example, Gaffke and Rüschendorf,
(1981), and Arvidsen and Johnsson, (1982) proposals have the same values for
Kendall’s $\tau$, but different ones for Spearman’s $\rho$.

Other examples of strict $d$-CTM vectors, partially covered by our
representation, can be found in Knott and Smith, (2006) , Lee and Ahn, (2014),
and after a linear transformation also the construction in Bubenik and
Holbrook, (2007) and references therein, can be seen as strict $d$-CTM
vectors.

The range of application of CTM constructions has been extended to other
marginal distributions. For instance, Rüschendorf and Uckelmann, (2002),
expands the work in Gaffke and Rüschendorf, (1981) to random variables
$\\{Y_{1},\ldots,Y_{d}\\}$ with unimodal distributions using the Levy-Shepp
form of the Khinchine representation theorem (Lévy,, 1962; Shepp,, 1962), as
$Y_{i}=XV_{i}$ where $V_{i}\sim U(-1,1)$ for all $1\leq i\leq d$. Hence, a CTM
construction for $\\{V_{1},\ldots,V_{d}\\}$ implies constant sum for for
$\\{Y_{1},\ldots,Y_{d}\\}$. Trivially, such a vector will achieve the smallest
variance for the sum of its components. In those cases, the literature refers
to these vectors as complete or joint mix, differentiating between having
identical or different marginals (see Puccetti and Wang,, 2015, and references
therein). Our general construction can be extend to non-uniform marginals
following Rüschendorf and Uckelmann, (2002).

In Rubinstein and Samorodnitsky, (1987), a different extension of Handscomb,
(1958) theorem was proposed by dropping the bijective condition for the
rearrangement. They prove the existence of antithetic solutions that minimize
the variance. According to the authors, optimal antithetic solutions should be
a function of only one uniform random variable, without restriction on the
functional dependence. Unfortunately, dropping the bijective condition result
in a tautological statement because, as shown in Brown, (1966), Whitt, (1976),
Vitale, (1990) and recently reformulated in theorem 1 of Puccetti and Wang,
(2015), every random vector can be expressed as a function of only one uniform
random variable. This difficulty of narrowing down conditions for the
existence of optimal antithetic variables is linked to the challenge of
extending Birkoff-von Neumman representation to the continuous case. It is, in
fact, well known that bijective rearrangements are only a sub-class of the
extremal transformations. For example, the $d=2$ case is known as Birkhoff’s
problem 111 (Isbell,, 1955), and even if there exists a characterization
(Lindenstrauss,, 1965), the necessary and sufficient conditions in their most
recognizable form (Moameni,, 2016) are of limited practical relevance. For a
discussion and an example of an extremal non-bijective class in the
multivariate case, refer to Durante et al., (2014). Since a general
characterization is out of reach, we solve the optimal transport problem for
transformations in the extremal class and produce a stochastic representation
that depends on a single standard uniform.

Historically, given the impossibility of obtaining optimal and feasible
antithetic plans, by mid 80’s the literature shifted the focus to negative
dependence. In particular, a procedure considered close to antithetic
sampling, but applicable to the general $d\geq 2$, $p\geq 1$ is the Latin
Hypercube sampling introduced in McKay et al., (1979).

######  Example 1.5 (McKay et al., (1979) Latin Hypercube).

Given a standard uniform $d$-dimensional random vector $\mathbf{U}$ and let
$\mathcal{D}^{\pi}_{t}=\left(\pi_{t}\left(0\right),\ldots,\pi_{t}\left(d-1\right)\right)^{T}$
be a permutation of $\\{0,1,\ldots,d-1\\}$ independent of
$\mathbf{U},\ldots,\mathbf{U}_{1}$ and

$\displaystyle\mathbf{U}_{1}=\dfrac{1}{d}\left(\mathcal{D}^{\pi}_{t}+\mathbf{U}\right)$
(9)

The simplicity of the method, the guarantee of asymptotic variance reduction
(Stein,, 1987) and the availability of a central limit theorem (Owen,, 1992)
made it one of the most common variance reduction strategies. The relationship
with antithetic variates was studied in Craiu and Meng, (2005), where through
the introduction of an iterative version of the method, the Iterated Latin
Hypercube (ILH), it is shown that, in the iteration limit, the resulting
random vector has an almost-sure constant sum. Our new representation allows
comparing ILH and its combination with other antithetic proposals. Finally, we
extend the central limit theorem in Owen, (1992), showing the irrelevance of
the starting distribution when $d$ goes to infinity, and the number of
iterations is fixed.

The paper is structured as follows: Section 2 contains the description of the
unified representation, its interpretation, and its invariant transformations.
Section 3 discusses distributional properties and concordance measures. New
and old illustrations of the unified representation and their ranking are
presented in Section 4 followed by the derivation of the general central limit
theorem for Latin Hypercube in Section 5. Numerical illustrations are
presented in Section 6 and the paper ends with a discussion of future
directions for research in Section 7.

##  2 Sampling on Line Segments

We introduce a general method for constructing antithetic vectors whose
components have a standard uniform, Uniform(0,1), as marginal distribution.
Non-uniform variables can be obtained using various transformations, e.g. the
inverse CDF method or the Levy-Shepp form of the Khinchine representation
theorem (Lévy,, 1962; Shepp,, 1962). We study conditions for achieving $d$-CTM
and show that several known countermonotonic random vectors used in variance
reduction, can be obtained as special cases of our general construction.

###  2.1 Standard Antithetic Construction

Our method relies on sampling with equal probability on a collection
$\mathcal{S}$ of line segments in the $d$-dimensional Euclidean space. Since
each segment is uniquely characterized by its endpoints or vertexes, the
collection $\mathcal{S}$ can be equivalently represented by the set of vertex
pairs that define the segments and their coordinates. This representation is
efficient in large dimensions even when the segments share some of their
vertexes.

More formally, let us define a vertex set
$\mathcal{V}=\left\\{1,\ldots,n\right\\}$ as a set of points in the
$d$-dimensional hypercube, the coordinates of the $k$-th vertex as the column
vector
$\mathbf{x}_{k}\equiv\left(x_{1k},\ldots,x_{dk}\right)^{T}\in\left[0,1\right]^{d}$
and the coordinate matrix
$\mathbf{X}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right)$ as the
collection of vertex coordinates. We assume there is an edge $e=(i,j)$ between
$i$ and $j$, with $i<j$, if there is a segment joining the two vertices $i$
and $j$, and define the collection of segments by the edge set
$\mathcal{E}=\left\\{\left(i,j\right)\in\mathcal{V}\times\mathcal{V}\right\\}$.
Then $\mathcal{G}=\left\\{\mathcal{V},\mathcal{E}\right\\}$ is an undirected
graph and $\mathcal{S}=\left\\{\mathcal{G},\mathbf{X}\right\\}$ is the
collection of segments. The lexicographic order on vertex indexing induces an
order on the edge set, such that the $k$-th element $e_{k}\in\mathcal{E}$ is
uniquely associated to its couple of vertices, defining the map
$\\{1,\ldots,\left|\mathcal{E}\right|\\}\mapsto\mathcal{E}$,
$k\rightarrow\left(i\left(k\right),j\left(k\right)\right)$.

Our stochastic construction relies on the graph representation $\mathcal{G}$
and requires a properly chosen vertex matrix $\mathbf{X}$ and two standard
uniform random numbers666As in the antithetic couple case it is also possible
to use only one standard uniform random number $W$ by setting
$V=\left\\{\right|\mathcal{E}\left|W\right\\}$.

  1. 1.

draw $V\sim\mathcal{U}[0,1]$ and $W\sim\mathcal{U}[0,1]$

  2. 2.

choose with uniform probability on the edge set $\mathcal{E}$ the edge $e_{K}$
by computing $K=\left\lfloor\left|\mathcal{E}\right|W\right\rfloor+1$

  3. 3.

obtain the random vertexes pair
$\left(I,J\right)=\left(i\left(K\right),j\left(K\right)\right)$

  4. 4.

obtain a random point on the segment joining vertexes $I$ and $J$ with uniform
probability

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle
x_{1I}V+\left(1-V\right)x_{1J}$ (10) $\displaystyle\vdots$ $\displaystyle
U_{d}$ $\displaystyle=$ $\displaystyle x_{dI}V+\left(1-V\right)x_{dJ}$

The following example shows that the standard antithetic method is a special
case of our general sampling construction.

######  Example 2.1.

Let us consider $d=2$ and sampling in the unit square on the diagonal joining
the vertex $1$ of coordinates
$\mathbf{x}_{1}\equiv\left(x_{11}=1,x_{21}=0\right)^{T}$ and the vertex $2$ of
coordinates $\mathbf{x}_{2}\equiv\left(x_{12}=0,x_{22}=1\right)^{T}$ (see
Figure 1). Let $V\sim\mathcal{U}[0,1]$ and compute :

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle
x_{11}V+\left(1-V\right)x_{12}=V$ $\displaystyle U_{2}$ $\displaystyle=$
$\displaystyle x_{21}V+\left(1-V\right)x_{22}=1-V$

$x_{2}$$x_{1}$$1$$1$$\mathbf{x}_{2}$$\mathbf{x}_{1}$$\mathbf{U}$$1-V$$V$ | $2$$e_{1}$$1$
---|---
Figure 1: Left: support set, i.e. the segment joining vertexes $1$ and $2$ of
coordinates $x_{1}$ and $x_{2}$, for the antithetic sampling. Right:
dependence graph $\mathcal{G}$ with $\mathcal{V}=\\{1,2\\}$ and
$\mathcal{E}=\\{e_{1}\\}$.

Thus, sampling one uniform antithetic couple $V$ and $1-V$ is equivalent to
sampling on a segment (see left plot in Fig. 1) and the support of the samples
can be summarized by the vertex coordinate matrix

$\displaystyle\mathbf{X}=\left(\begin{array}[]{cc}x_{11}&x_{12}\\\
x_{21}&x_{22}\end{array}\right)$ (13)

and the couple of vertexes $e_{1}=\left(1,2\right)$ of the segment we are
sampling on (right plot).

The marginal uniformity of the samples on the segment follows from the convex
combination representation and the standard uniform assumption for $V$

$\displaystyle\begin{array}[]{cc}U_{l}\sim\mathcal{U}\left[\min\left(x_{l1},x_{l2}\right),\max\left(x_{l1},x_{l2}\right)\right]&l=1,2\end{array}$
(15)

In addition, the standard uniformity of $U_{l}$ follows from the assumption
$\max\left(x_{l1},x_{l2}\right)=1$ and $\min\left(x_{l1},x_{l2}\right)=0$ for
all $l=1,2$, and the $d$-CTM property from conditions on the vertex
coordinates: $x_{11}+x_{21}=x_{12}+x_{22}=1$.

Example 2.1 illustrates the fact, which is easy to prove in full generality,
that our construction generates samples with uniform probability on
$\mathcal{S}$. However, it does not guarantee that all the marginal
distributions of the components $U_{l}$, for $l\in\mathcal{D}$, are
$\mbox{Uniform}(0,1)$. In the following, we study the conditions on
$\mathcal{G}$ and $\mathbf{X}$, such that the variables $U_{l}$,
$l\in\mathcal{D}$ are conditionally uniform, given the choice of the segment,
and marginally $\mbox{Uniform}(0,1)$.

###  2.2 Uniformity

We provide conditions on the collection of segments $\mathcal{S}$ in
$d$-dimensional hypercube to achieve standard marginal uniformity and $d$-CTM
when using our construction. Before introducing the general result we
introduce some notation and discuss the main assumptions.

For the general case, the following condition rules out atomic and mixed
probability measures.

######  Assumption 1 (Admissibility).

The set $\mathcal{S}=\\{\mathcal{G},\mathbf{X}\\}$ is admissible if all
segments in $\mathcal{S}$ are not contained in any of the $(d-1)$-hyperplanes
that are parallel to a $(d-1)$-dimensional hyperface of the unit hypercube
$[0,1]^{d}$.

We provide some intuition for Assumption 1 through the following
$2$-dimensional example.

######  Example 2.2.

Consider the collection segments in the left plot of Figure 2 with coordinate
matrix:

$\displaystyle\mathbf{X}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\alpha&\beta&\gamma&\alpha\\\
\beta&\alpha&\beta&\alpha\end{array}\right)$ (18)

where $\alpha<\beta\leq\gamma\in\mathbb{R}$. According to the lexicographic
order map $k\rightarrow\left(i\left(k\right),j\left(k\right)\right)$ the edge
set $\mathcal{E}^{*}$ is:

$\displaystyle\left\\{e^{*}_{1}=\left(1,2\right),e^{*}_{2}=\left(1,3\right),e^{*}_{3}=\left(1,4\right),e^{*}_{4}=\left(2,3\right),e^{*}_{5}=\left(2,4\right),e^{*}_{6}=\left(3,4\right)\right\\}$

(see right plot).

$x_{2}$$x_{1}$$1$$1$$\mathbf{x}_{2}$$\mathbf{x}_{1}$$\mathbf{x}_{3}$$\mathbf{x}_{4}$ |  $2$$e^{*}_{1}$$1$$4$$e^{*}_{6}$$3$$e^{*}_{4}$$e^{*}_{5}$$e^{*}_{2}\quad$$e^{*}_{3}$
---|---
Figure 2: Left: support set of the sampling scheme. Right: admissible (solid)
and not-admissible (dashed) edges. Edge labeling follows the vertexes
lexicographic order.

Since the sampling method can concentrate the probability mass at some points
or along some directions of the hypercube we need to impose some admissibility
conditions to have non-degenerate distributions. The first set of conditions
excludes degenerate segments, that are segments with equal end-points. Thus,
we rule out self-loops from the graph, i.e. edges from one vertex to itself.
The second set of admissibility conditions excludes edges where the
distribution concentrates along some coordinates. The edge $e^{*}_{2}$ joining
vertex 1 to 3 is not admissible since:

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle
x_{11}V+x_{13}\left(1-V\right)=\gamma-\left(\gamma-\alpha\right)V$
$\displaystyle U_{2}$ $\displaystyle=$ $\displaystyle
x_{21}V+x_{23}\left(1-V\right)=\beta$

and $U_{2}$ is almost surely constant conditionally on being on that edge. A
similar remark applies to the edges $e^{*}_{3}$ and $e^{*}_{5}$, whereas
$e^{*}_{4}$ (red dashed in Figure 2) is only admissible if $\beta<\gamma$.

In summary, the admissible edge set is:

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ccc}\left\\{e_{1}=e^{*}_{1}=\left(1,2\right),e_{2}=e^{*}_{4}=\left(2,3\right),e_{3}=e^{*}_{6}=\left(3,4\right)\right\\},&&\mbox{
if }\beta<\gamma\\\ &&\\\
\left\\{e_{1}=e^{*}_{1}=\left(1,2\right),e_{2}=e^{*}_{6}=\left(3,4\right)\right\\},&&\mbox{
if }\beta=\gamma\end{array}\right.$

In our construction, conditional uniformity is a necessary condition for
standard marginal uniformity and the admissibility assumption implies the
conditional uniformity. To enhance the paper’s readability, all proofs are
deferred to the Appendix.

######  Lemma 2.1 (Conditional Uniformity).

Let $\mathcal{S}$ satisfy Assumption 1. Conditionally on being on the $k$-th
segment of edge $e_{k}=\left(i\left(k\right),j\left(k\right)\right)$, for each
$l\in\mathcal{D}$, the random variable $U_{l}$ in (10) is uniform on
$\left[\alpha_{l,k},\beta_{l,k}\right]$ with
$\alpha_{l,k}=\min\left(x_{l,i\left(k\right)},x_{l,j\left(k\right)}\right)$
and
$\beta_{l,k}=\max\left(x_{l,i\left(k\right)},x_{l,j\left(k\right)}\right)$.

Another requirement for our construction is that sampling points are in the
unit hypercube which is guaranteed by the following range condition.

######  Assumption 2 (Range).

The range requirement is satisfied if $\max\\{x_{lk},k=1,\ldots,n\\}=1$ and
$\min\\{x_{lk},k=1,\ldots,n\\}=0$.

Since $U_{l}$ is a convex combination of $x_{l,i}$,
$i\in\left\\{1,\ldots,n\right\\}$, Assumption 2 is needed in light of the
requirement that $U_{l}\in\left[0,1\right]$ for each $l\in\mathcal{D}$.

In order to introduce the third assumption we need some notation. The set of
admissible edges depends on the partitions induced by the distinct elements,
sorted in ascending order, in the rows of $\mathbf{X}$. For each row
$\left(x_{l,1},\ldots,x_{l,n}\right)$, of $\mathbf{X}$ with $l=1,\ldots,d$ we
define
$\mathbf{a}_{l}=\left(a_{l,1},a_{l,2},\ldots,a_{l,n_{l}-1},a_{l,n_{l}}\right)$
the sequence of $n_{l}\leq n$ sorted distinct elements. The unique values
define a partition, say $\left\\{A_{l,m}\right\\}^{n_{l}}_{m=2}$ of the unit
interval with elements:

$A_{l,m}=\left[a_{l,m-1},a_{l,m}\right),\quad
m\in\left\\{2,\ldots,n_{l}\right\\}.$ (20)

For each unique value $a_{lk}$, the position set
$\mathcal{M}_{l,k}=\left\\{i\in\left\\{1,\ldots,n\right\\}:x_{l,i}=a_{l,k}\right\\}$,
$k=1,\ldots,n_{l}$ denotes the set of points that provides the unique
projected values of $\mathbf{X}$’s columns into the $l$-th dimension. For each
row, $\mathcal{M}_{l,k}$ must therefore satisfy:

$\mathcal{M}_{l,m}\cap\mathcal{M}_{l,m^{\prime}}=\emptyset,\quad\bigcup^{n_{l}}_{m=1}\mathcal{M}_{l,m}=\left\\{1,\ldots,n\right\\},$

and one can represent the coordinates of the vertex $\mathbf{x}_{k}$ by using
the sets of positions of the unique values:

$\displaystyle
x_{l,k}=\sum^{n_{l}}_{m=1}a_{l,m}\mathbb{I}_{\mathcal{M}_{l,m}}\left(k\right),\,l\in\mathcal{D},$

which will be used in the following to write the CM constraint.

For illustration purposes, let us focus on the first coordinate of Example 2.2
in the case $\beta<\gamma$. The sequence of $n_{l}=3$ sorted distinct elements
of the first row is
$\mathbf{a}_{1}=\left(a_{1,1}=\alpha,a_{1,2}=\beta,a_{1,3}=\gamma\right)$ and
the sets of positions of the unique values are
$\mathcal{M}_{1,1}=\left\\{1,4\right\\}$,
$\mathcal{M}_{1,2}=\left\\{2\right\\}$, $\mathcal{M}_{1,3}=\left\\{3\right\\}$
which is a partition of $\\{1,2,3,4\\}$.

We denote with
$\mathcal{G}_{l}\equiv\left\\{\mathcal{V}_{l},\mathcal{E}_{l}\right\\}$ the
projection of $\mathcal{G}$ on the $l$-coordinate induced by the unique values
$\mathbf{a}_{l}$. $\mathcal{G}_{l}$ is graph obtained by assigning a node to
each of the $n_{l}$ components of $\mathbf{a}_{l}$ and defining the edge set
$\mathcal{E}_{l}=\\{e_{l,1},\ldots,e_{l,|\mathcal{E}|}\\}$ through the map
$\\{1,\ldots,\left|\mathcal{E}\right|\\}\mapsto\mathcal{E}_{l}$,
$k\in\\{1,\ldots,\left|\mathcal{E}\right|\\}\rightarrow
e_{l,k}=\left(m\left(k\right),m^{\prime}\left(k\right)\right)\in\mathcal{E}_{l}$
where

$\displaystyle m\left(k\right)$ $\displaystyle=$
$\displaystyle\left\\{m\in\left\\{1,\ldots,n_{l}\right\\}:i\left(k\right)\in\mathcal{M}_{l,m}\right\\},$
$\displaystyle m^{\prime}\left(k\right)$ $\displaystyle=$
$\displaystyle\left\\{m^{\prime}\in\left\\{1,\ldots,n_{l}\right\\}:j\left(k\right)\in\mathcal{M}_{l,m^{\prime}}\right\\}.$

We call $n^{l}_{\left(m,m^{\prime}\right)}$ the edges multiplicity between the
projected nodes $m$ and $m^{\prime}$, by convention if
$\left(m,m^{\prime}\right)$ is not in $\mathcal{E}_{l}$ then
$n^{l}_{\left(m,m^{\prime}\right)}=0$.

######  Remark 2.1.

In the proposed notation, Assumption 1 corresponds to $m\left(k\right)\neq
m^{\prime}\left(k\right)$ when
$e_{l,k}=\left(m\left(k\right),m^{\prime}\left(k\right)\right)$ for all
$k=1,\ldots,|\mathcal{E}|$, $l=1,\ldots,d$.

In Figure 3 we report the projection of the graph of Figure 2 onto the $l$-th
coordinates given in Equation (18) of Example 2.2.

$2$$e_{1,1}$$3$$e_{1,3}$$e_{1,2}$$1$ | $2$$e_{2,2}$$e_{2,1}$$e_{2,3}$$1$
---|---
$2$$e_{1,1}$$e_{1,2}$$1$ | $2$$e_{2,1}$$e_{2,2}$$1$
Figure 3:  Graph projection for the graph in Example 2.2 in the case
$\beta<\gamma$ (top) and $\beta=\gamma$ (bottom), considering the first (left)
or the second (right) coordinate.

Set $\alpha_{l,k}=a_{l,m_{\alpha}\left(k\right)}$ and
$\beta_{l,k}=a_{l,m_{\beta}\left(k\right)}$, where $m_{\alpha}\left(k\right)$
and $m_{\beta}\left(k\right)$ are the minimum and the maximum value between
$m\left(k\right)$ and $m^{\prime}\left(k\right)$, and define

$\mathcal{K}_{l,m}\equiv\left\\{k\in\left\\{1,\ldots,\left|\mathcal{E}\right|\right\\}:m_{\alpha}\left(k\right)+1\leq
m\leq m_{\beta}\left(k\right)\right\\}$

We state now the condition on the coordinates of the sampling construction for
marginal standard uniformity.

######  Assumption 3 (Coordinate).

The following set of $n_{l}-2$ equations in the variables $a_{l,m}$,
$m=2,\ldots,n_{l}-1$ are satisfied

$\displaystyle F_{l,m}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{l,m}}\dfrac{1}{a_{l,m_{\beta}\left(k\right)}-a_{l,m_{\alpha}\left(k\right)}}-1=0$
(21)

$m=2,\ldots,n_{l}-1$.

The following example clarifies the relationship between Assumptions 2-3 and
standard uniformity.

######  Example 2.3.

We discuss separately the two cases: $\beta=\gamma$ and $\beta<\gamma$ since
they correspond to two different admissible edge sets. If $\beta=\gamma$ the
first component of $\mathbf{U}$ in the stochastic construction of Equation
(10) is:

$\displaystyle
U_{1}=\left\\{\begin{array}[]{ccc}\alpha+\left(\beta-\alpha\right)V&\hbox{if}&K=1\\\
\beta-\left(\beta-\alpha\right)V&\hbox{if}&K=2\end{array}\right.$

Since we are sampling uniformly on the edge set $\mathcal{E}$ it follows that:

$\displaystyle\mathbb{P}\left(K=1\right)=\mathbb{P}\left(K=2\right)=\dfrac{1}{\left|\mathcal{E}\right|}=\dfrac{1}{2}$

and $U_{1}$ has marginal PDF

$\displaystyle f\left(u_{1}\right)$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\dfrac{1}{\beta-\alpha}\mathbb{I}_{\left[\alpha,\beta\right]}\left(u_{1}\right)+\dfrac{1}{2}\dfrac{1}{\beta-\alpha}\mathbb{I}_{\left[\alpha,\beta\right]}\left(u_{1}\right)=\dfrac{1}{\beta-\alpha}\mathbb{I}_{\left[\alpha,\beta\right]}\left(u_{1}\right).$
(23)

We get a standard uniform random variable if $\alpha=0$ and $\beta=1$. By the
same argument, $U_{2}$ is a standard uniform random variable. In conclusion,
by sampling on the two diagonals of the unit square with a mixture of an
antithetic couple $\left(V,1-V\right)$ and a comonotonic couple
$\left(V,V\right)$, the method can attain marginal standard uniformity (left
plot in Figure 4).

$x_{2}$$x_{1}$$1$$1$$\mathbf{x}_{2}$$\mathbf{x}_{1}$$\mathbf{x}_{3}$$\mathbf{x}_{4}$ |  $x_{2}$$x_{1}$$1$$1$$\mathbf{x}_{1}$$\mathbf{x}_{2}$$\mathbf{x}_{3}$$\mathbf{x}_{4}$
---|---
Figure 4: Sampling on line segments. The support set of the 2 dimensional
example with 4 vertexes. Cases: $\beta=\gamma$ (left) and $\beta<\gamma$
(right).

If $\beta<\gamma$ the PDF of $U_{1}$ is:

$\displaystyle f\left(u_{1}\right)$ $\displaystyle=$ $\displaystyle
1/3\dfrac{1}{\beta-\alpha}\mathbb{I}_{\left[\alpha,\beta\right]}\left(u_{1}\right)+1/3\dfrac{1}{\gamma-\beta}\mathbb{I}_{\left[\beta,\gamma\right]}\left(u_{1}\right)+1/3\dfrac{1}{\gamma-\alpha}\mathbb{I}_{\left[\alpha,\gamma\right]}\left(u_{1}\right).$

A necessary and sufficient condition for $u_{1}\in\left[0,1\right]$ is
$\alpha=\min\\{x_{1k},k=1,\ldots,4\\}=0$ and
$\gamma=\max\\{x_{1k},k=1,\ldots,4\\}=1$ which implies the PDF is piece-wise
constant on the partition $\left[0,\beta\right)\subset[0,1]$ and
$\left[\beta,1\right]\subset[0,1]$, induced by the unique values in the first
row of $\mathbf{X}$, and is null on
$\left[0,\beta\right)\cap\left[\beta,1\right]$. The value of $\beta$ such that
the PDF takes the same value on the elements of the partition, i.e.:

$\displaystyle\left\\{\begin{array}[]{ccc}\dfrac{1}{\beta}+1=3&\hbox{if}&u_{1}\in\left[0,\beta\right)\\\
\dfrac{1}{1-\beta}+1=3&\hbox{if}&u_{1}\in\left[\beta,1\right]\\\
\end{array}\right.$ (26)

is $\beta=1/2$ which implies $U_{1}\sim\mathcal{U}[0,1]$ and
$U_{2}\sim\mathcal{U}[0,1/2]$. Thus, for $\beta<\gamma$ our construction is
not standard uniform along all coordinates of the vector. The right plot of
Figure 4 shows that the range of $U_{1}$ is $[0,1]$ (horizontal axis) and the
range of $U_{2}$ is $[0,1/2]$ (vertical axis) 777Similar arguments can be
applied to show that imposing standard uniformity for $U_{2}$ requires
$\min\\{x_{1k},k=1,\ldots,4\\}=\alpha=0$ and
$\max\\{x_{1k},k=1,\ldots,4\\}=\beta=1$ which is not satisfied since
$\beta<\gamma$.

We are ready to state the main result of this section which proves the
marginal standard uniformity of our construction. In addition, we provide a
method to find $\mathbf{a}_{l}$ satisfying the condition in Assumption 3.

######  Theorem 2.1 (Marginal Standard Uniformity).

Under Assumptions 1-3, each coordinate of the random vector
$\mathbf{U}=(U_{1},\ldots,U_{d})$ in the stochastic representation (10) has a
$\mbox{Uniform}(0,1)$ marginal distribution.

######  Theorem 2.2.

The system in (21) is equivalent to the following convex minimization problem

$\displaystyle\displaystyle\min_{\left(a_{l,2},\ldots,a_{l,n_{l}-1}\right)\in\left[0,1\right]^{n_{l}-2}}\Psi_{l}\left(\mathbf{a}_{l}\right),$
(27)

with

$\displaystyle\Psi_{l}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle-\dfrac{1}{2\left|\mathcal{E}\right|}\displaystyle\sum^{n_{l}}_{m,m^{\prime}=1}n^{l}_{\left(m,m^{\prime}\right)}\log\left|a_{l,m^{\prime}}-a_{l,m}\right|.$

Since the optimization problem is convex, if a solution exists it is a global
minimum. In addition, since $\Psi_{l}\left(\mathbf{a}_{l}\right)$ is a sum of
lower semi-continuous functions, it is lower semi-continuous. This, together
with the compactness of the unit hypercube, guarantees the existence of a
solution, by the lower version of Weierstrass theorem (see for example Theorem
2.43 in Aliprantis and Border, (2007)).

In the next theorem we show that our construction satisfies an optimality
criteria involving the Kullback–-Leibler (KL) divergence from the uniform
distribution. We remind the reader that the KL divergence of the probability
measure $\mathbb{P}$ with respect to the probability measure $\mathbb{Q}$ is

$\displaystyle D_{KL}\left(\mathbb{P}||\mathbb{Q}\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\mathbb{P}}\left[\log\left(\dfrac{d\mathbb{P}}{d\mathbb{Q}}\right)\right].$

In our case $\mathbb{P}$ is the joint measure of $U_{l}$ and $K$ given by the
stochastic representation (10):

$\displaystyle d\mathbb{P}\left(u_{l},k\right)$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}f\left(u_{l}\right|\left.K=k\right)du_{l}$

$\mathbb{Q}$ is the joint uniform independent measure in $U_{l}$ and $K$:

$\displaystyle d\mathbb{Q}\left(u_{l},k\right)$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\mathbb{I}_{\left[0,1\right]}\left(u_{l}\right)du_{l}$

######  Theorem 2.3.

The minimization problem in (27) is equivalent to minimizing the Kullback-
Leibler divergence of the joint measure of $U_{l}$ and $K$ given by the
stochastic representation (10) from a joint uniform independent measure in
$U_{l}$ and $K$.

The results in Theorems 2.1-2.3 allow us to find the $l$-th row
$\mathbf{a}_{l}$, of the coordinates matrix $X$ such that $U_{l}$ in
representation (10) is standard uniform. Thus, we have $d$ minimization
problems each with a different number of variables, $n_{l}$, where the
dependence structure in $\mathbf{U}$ is given by the graph $\mathcal{G}$ of
the segments in $\mathcal{S}$ and marginally encoded by $\mathbf{a}_{l}$ and a
projected graphs $\mathcal{G}_{l}$.

######  Definition 2.1.

Let $\mathcal{S}=(\mathcal{G},\mathbf{X})$ be the segment set with
$\mathcal{G}=\left\\{\mathcal{V},\mathcal{E}\right\\}$. The convex
minimization problem

$\displaystyle\displaystyle\min_{\left\\{a_{l,m},m=1,\ldots,n_{l},\,l=1,\ldots,d\right\\}\in\left[0,1\right]^{n_{1}+\ldots+n_{d}}}\sum^{d}_{l=1}\Psi_{l}\left(\mathbf{a}_{l}\right),$

is called standard uniform on $\mathcal{S}$ problem.

###  2.3 Strict Countermonotonicity on Segments

Limiting the study to strict $d$-CTM leads to an unique value for the constant
$k$ in (4):

$\displaystyle
k=\mathbb{E}\left[\sum^{d}_{j=1}U_{j}\right]=\sum^{d}_{j=1}\mathbb{E}\left[U_{j}\right]=\dfrac{d}{2}$

The constant sum condition can be written as a linear restriction on the
coordinates of the vertexes $\mathbf{x}_{k}$, that is

$\displaystyle\sum^{d}_{l=1}U_{l}$ $\displaystyle=$
$\displaystyle\sum^{d}_{l=1}x_{lJ}+V\left[\sum^{d}_{l=1}x_{lI}-\sum^{d}_{l=1}x_{lJ}\right]=\dfrac{d}{2}$

and since the previous relationship should be valid for all $V$ and
$\left(I,J\right)$ (i.e. for all $W$ in our setting) we obtain the condition
that all vertexes should be in the hyperplane of constant sum, i.e.

$\displaystyle\begin{array}[]{ccc}\displaystyle\sum^{d}_{l=1}x_{lk}=\sum^{d}_{l=1}\sum^{n_{l}}_{m=1}a_{l,m}\mathbb{I}_{\mathcal{M}_{l,m}}\left(k\right)=\dfrac{d}{2}&&k=1,\ldots,n.\end{array}$
(29)

The convexity of the minimization problem in Eq. (2.1) is not altered by the
inclusion of a linear constraint, nevertheless the constraint couples the
coordinates and yields the following non-separable optimization problem in
$n_{d}=\sum^{d}_{l=1}n_{l}$ variables.

######  Definition 2.2.

Let $\mathcal{S}=(\mathcal{G},\mathbf{X})$ be the segment set with
$\mathcal{G}=\left\\{\mathcal{V},\mathcal{E}\right\\}$. The convex
minimization problem:

$\displaystyle\displaystyle\min_{\left\\{a_{l,m},m=1,\ldots,n_{l},\,l=1,\ldots,d\right\\}\in\left[0,1\right]^{n_{1}+\ldots+n_{d}}}\sum^{d}_{l=1}\Psi_{l}\left(\mathbf{a}_{l}\right)$

subject to:

$\displaystyle\sum^{d}_{l=1}\sum^{n_{l}}_{m=1}a_{l,m}\mathbb{I}_{\mathcal{M}_{l,m}}\left(k\right)=\dfrac{d}{2},$
$\displaystyle k=1,\ldots,n$

will be referred as the strict $d$-CTM on $\mathcal{S}$ problem.

Since all the constraints are affine the problem in the above definition
represents an ordinary convex problem in the terminology of Rockafellar,
(1970) and can be solved using the method of Lagrange multipliers. Local
minima are also global if they exist. Finally, the existence of a solution is
guaranteed by lower semi-continuity of
$\sum^{d}_{l=1}\Psi_{l}\left(\mathbf{a}_{l}\right)$ and the fact that the
intersection of the unit hypercube with the hyperplane of constant sum is
compact.

###  2.4 Antithetic Vector Transformations

We introduce some transformations of the random vector of representation (10)
which are useful to study the properties of our sampling on line segments.
These transformations preserve uniformity and provide a general framework to
the construction of antithetic vectors. The use of transformations to span the
set of antithetic vectors was pioneered by Andréasson, (1972) and Andréasson
and Dahlquist, (1972) and discussed further by Roach and Wright, (1977).
First, we define deterministic compositions and present reflections as a
particular case. We provide sufficient conditions for preserving strict
countermonotoniticy under deterministic compositions. Secondly, we introduce
stochastic compositions as particular instances of probabilistic mixtures that
preserve uniformity and the strict countermonotonicity property. In addition,
we show that permutations preserve strict countermonotonicity. Finally, we
generalize of our stochastic representation (10) inspired by the deterministic
composition.

The notion of deterministic composition relies on the interpretation of the
components of the stochastic representation (10) as a deterministic map on the
unit interval. Let $\mathbf{X}$ be a solution to the standard uniform on
$\mathcal{S}^{(x)}=(\mathcal{G}^{(x)},\mathbf{X})$ with edge set
$\mathcal{E}^{(x)}$ and let
$T_{x,l,k}:\left[0,1\right]\mapsto\left[0,1\right]$

$T_{x,l,k}\left(u\right)=x_{li(k)}u+x_{lj(k)}\left(1-u\right).$ (30)

Then each element of (10) can be written as

$U_{l}=T_{x,l,K}\left(V\right)$ (31)

with $U_{l}$, $K$, and $V$ as in (10). Let $\mathbf{Y}$ be a different
solution to the standard uniform
$\mathcal{S}^{(y)}=(\mathcal{G}^{(y)},\mathbf{Y})$ with edge set
$\mathcal{E}^{(y)}$. It induces a map
$T_{y,l,k}:\left[0,1\right]\mapsto\left[0,1\right]$

$T_{y,l,k}\left(u\right)=y_{li(k)}u+y_{lj(k)}\left(1-u\right).$ (32)

The deterministic composition is defined as $T_{y,l,k^{\prime}}\circ
T_{x,l,k}\left(u\right)=T_{y,l,k}\left(T_{x,l,k}(u)\right)$. The following
lemma states that the random vector with components $T_{y,l,k^{\prime}}\circ
T_{x,l,k}\left(U_{l}\right)$ can be represented as in (10) and provides the
set of segments associated with the composition.

######  Lemma 2.2.

The deterministic composition $T_{z,l,k}\left(V\right)=T_{y,l,k^{\prime}}\circ
T_{x,l,k}\left(V\right)$, $l\in\mathcal{D}$ is uniform on the segments
$\mathcal{S}^{\prime}=(\mathcal{G}^{\prime},\mathbf{Z})$ with graph
$\mathcal{G}^{\prime}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime})$, vertex set
$\mathcal{V}^{\prime}=\mathcal{E}^{(y)}\times\mathcal{V}^{(x)}$, and edge set
$\mathcal{E}=\mathcal{E}^{(y)}\times\mathcal{E}^{(x)}$. The coordinate matrix
$\mathbf{Z}$ has elements

$\displaystyle
z_{l,m}=\left(y_{li\left(k\right)}x_{li^{\prime}}+y_{lj\left(k\right)}\left(1-x_{li^{\prime}}\right)\right),$

where the index $m\in\mathcal{V}^{\prime}$ corresponds to the pair
$\left(k,i^{\prime}\right)$ with
$k\in\left\\{1,\ldots,\left|\mathcal{E}^{(y)}\right|\right\\}$ and
$i^{\prime}\in\\{1,\ldots,|\mathcal{V}^{(x)}|\\}$ sorted in lexicographic
order.

We provide a sufficient condition for the deterministic composition to satisfy
the constant sum constraint and leave the necessary conditions for future
research.

######  Lemma 2.3.

Let $\mathbf{X}$ be a solution of the strict countermonotonic problem on
$\mathcal{S}^{(x)}=(\mathcal{G}^{(x)},\mathbf{X})$ and $\mathbf{Y}$ a solution
of the standard uniform problem on
$\mathcal{S}^{(y)}=(\mathcal{G}^{(y)},\mathbf{Y})$ with the additional
constraints that for every $l\in\mathcal{D}$ and
$k\in\left\\{1,\ldots,\mathcal{E}^{(y)}\right\\}$

$\displaystyle|y_{li\left(k\right)}-y_{lj\left(k\right)}|=c_{1},\quad\sum^{d}_{l=1}y_{lj\left(k\right)}=c_{2}.$

Let $W_{l}=T_{y,l,k^{\prime}}\circ T_{x,l,k}\left(V\right)$, $l\in\mathcal{D}$
then $\sum^{d}_{l=1}W_{l}=d/2$ if and only if

$\left(1-c_{1}\right)\dfrac{d}{2}=c_{2}.$ (33)

######  Remark 2.2.

$\mathbf{Y}$ in Lemma 2.3 cannot be a strict countermonotonic solution on
$\mathcal{S}^{(y)}=(\mathcal{G}^{(y)},\mathbf{Y})$, because this implies
$c_{2}=\dfrac{d}{2}$ and $c_{1}=2$ while the maximum value of $c_{1}$ is $1$
for all solutions of standard uniform problems by the Assumption 2.

Let $I_{d}$ denote the identity matrix and $\mathbf{e}_{l}$ its $l$-th column.
Reflections are defined in the following.

######  Definition 2.3 (Reflections).

Let $\mathbf{U}\in[0,1]^{d}$ be a random vector, the transformed vector
$\mathbf{U}\mapsto\mathbf{W}=R_{l,\frac{1}{2}}\left(\mathbf{U}\right)$ with
$R_{l,\frac{1}{2}}\left(\mathbf{U}\right)=\left[I_{d}-2\left(\mathbf{e}_{l}\mathbf{e}^{T}_{l}\right)\right]\mathbf{U}+\mathbf{e}_{l}$,
defines a reflection of $\mathbf{U}$ through an hyperplane with $l$-coordinate
equal to 1/2. Given a index subset $\mathcal{L}\subseteq\mathcal{D}$, the
sequential reflection transformation
$R_{\mathcal{L},\frac{1}{2}}(\mathbf{U})=[I_{d}-2(\sum_{l\in\mathcal{L}}\mathbf{e}_{l}\mathbf{e}^{T}_{l})]\mathbf{U}-\sum_{l\in\mathcal{L}}\mathbf{e}_{l}$
is the transformation obtained by sequentially reflecting $\mathbf{U}$ for a
sequence of indexes $l\in\mathcal{L}$ through a hyperplane with $l$-coordinate
equal to 1/2.
$R_{\mathcal{D},\frac{1}{2}}(\mathbf{U})=(I_{d}-2I_{d})\mathbf{U}+\mathbf{1}_{d}=(1-U_{1},\ldots,1-U_{l},\ldots,1-U_{d})$
is the central inversion through the center of the unit hypercube.

If $\mathbf{U}$ is uniform on segments, the reflection can be interpreted as a
deterministic composition. Let $l\in\mathcal{D}$ and

$\displaystyle
y_{l,1}=\left\\{\begin{array}[]{ccc}1&\hbox{if}&l\in\mathcal{L}\\\
0&\hbox{if}&\hbox{otherwise}\end{array}\right.$ $\displaystyle
y_{l,2}=\left\\{\begin{array}[]{ccc}0&\hbox{if}&l\in\mathcal{L}\\\
1&\hbox{if}&\hbox{otherwise}\end{array}\right.$

Using the definition of reflection we can write the elements of
$\mathbf{W}=R_{\mathcal{L},\frac{1}{2}}(\mathbf{U})$ as:

$\displaystyle W_{l}=y_{l,1}U_{l}+y_{l,2}\left(1-U_{l}\right).$

If $\mathbf{U}$ is strict $d$-CTM on segments and since $c_{1}=1$ and
$c_{2}=d-\left|\mathcal{L}\right|$, according to Lemma 2.3 $\mathbf{W}$ is
$d$-CTM only if $\mathcal{L}=\mathcal{D}$.

The following notion of stochastic composition augments the class of
deterministic compositions and it is useful for constructing $d$-CTM vectors
on segments.

######  Lemma 2.4 (Stochastic Composition).

Let $\mathbf{U}^{t}$ , $t=1,2$ be strict $d$-CTM vectors on
$\mathcal{S}^{t}=(\mathcal{G}^{t},\mathbf{X}^{t})$, with $\mathcal{E}^{t}$ the
edge set of $\mathcal{G}_{t}$, the random vector

$\mathbf{W}=\left\\{\begin{array}[]{ccc}\mathbf{U}^{1},&\hbox{with
probability}&\dfrac{\left|\mathcal{E}^{1}\right|}{\left|\mathcal{E}^{2}\sqcup\mathcal{E}^{1}\right|}\\\
\mathbf{U}^{2},&\hbox{with
probability}&\dfrac{\left|\mathcal{E}^{2}\right|}{\left|\mathcal{E}^{1}\sqcup\mathcal{E}^{2}\right|}\end{array}\right.,$

defines a stochastic composition of $U^{1}$ and $U^{2}$ and is strict $d$-CTM
on
$\mathcal{S}^{\prime}=(\mathcal{G}^{1}\sqcup\mathcal{G}^{2},\left(\mathbf{X}^{1},\mathbf{X}^{2}\right))$,
where $\sqcup$ denotes the disjoint union of graphs.

Another transformation used later on in this paper to obtain exchangeable
random vectors is the following.

######  Lemma 2.5 (Permutation).

Let $\mathbf{U}\in[0,1]^{d}$ be a random vector and
$\left\\{\pi\left(1\right)\right.,$ $\ldots,$
$\left.\pi\left(d\right)\right\\}$ a permutation of $\mathcal{D}$ with
$P_{\pi}=\left(\mathbf{e}_{\pi\left(1\right)},\ldots,\mathbf{e}_{\pi\left(d\right)}\right)^{T}$
the corresponding row-permutation matrix. $\mathbf{W}$ is a permutation of
$\mathbf{U}$ if $\mathbf{W}=P_{\pi}\mathbf{U}$. If $\mathbf{U}$ is strict
$d$-CTM on $\mathcal{S}=(\mathcal{G},\mathbf{X})$ then $\mathbf{W}$ is strict
$d$-CTM on $\mathcal{S}^{\prime}=(\mathcal{G},P_{\pi}\mathbf{X})$.

######  Remark 2.3.

From the properties of transformations in Lemma 2.4 and 2.5 it follows that
the stochastic combination of all the permutations of a strict
countermonotonic vector is strict countermonotonic and exchangeable.

Using the inner composition we can rephrase standard uniform on segments as
the deterministic composition of the maps induced by collection of segments
$\mathcal{S}$ with the comonotonic random vector $U_{l}=V$ for each
$l\in\mathcal{D}$ which is a uniform vector on the principal diagonal of the
hypercube. In particular, it can be shown that our construction includes the
multivariate shuffles of min (Mikusiński and Taylor,, 2010). This new
interpretation of Eq. (10) allows for a generalization to any random vectors
with marginal standard uniforms. Given a $d$-dimensional $\mathbf{V}\sim F$
with $V_{l}\sim\mathcal{U}[0,1]$ for $l=1,\ldots,d$ and the stochastic
representation

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle
x_{1I}V_{1}+x_{1J}\left(1-V_{1}\right)$ (36) $\displaystyle\vdots$
$\displaystyle U_{d}$ $\displaystyle=$ $\displaystyle
x_{dI}V_{d}+x_{dJ}\left(1-V_{d}\right),$

we have the following corollary to Theorem 2.1.

######  Corollary 1.

Under Assumptions 1-3 the $l$-th component of the random vector
$\mathbf{U}=(U_{1},\ldots,U_{d})$ in the stochastic representation (36) is a
standard uniform random variable.

##  3 Distributional Properties

###  3.1 Distribution

The joint cumulative distribution function of the stochastic representation
$\mathbf{U}$ of Eq. (36) can be written using the distributions of the
reflections of $\mathbf{V}$. Given a index subset $\mathcal{L}$, we denote the
distribution of the reflection
$R_{\mathcal{L},\frac{1}{2}}\left(\mathbf{V}\right)$ with:

$\displaystyle
F_{\mathbf{V},\mathcal{L}}\left(\mathbf{u}\right)=\mathbb{P}\left(R_{\mathcal{L},\frac{1}{2}}\left(\mathbf{V}\right)\leq\mathbf{u}\right).$

The marginal distribution of $U_{l}$ conditionally on living on an edge
$e_{k}$ is $\mathcal{U}[\alpha_{l,k},\beta_{l,k}]$ with CDF

$F_{U_{l}|K}\left(u_{l};k\right)=\frac{\max\\{\alpha_{l,k},\min\\{\beta_{l,k},u_{l}\\}\\}-\alpha_{l,k}}{\beta_{l,k}-\alpha_{l,k}}$
(37)

where $\alpha_{l,k}$, $\beta_{l,k}$ are defined in Lemma 2.1. For each
$k=1,\ldots,n$ define the sets:

$\displaystyle\mathcal{L}_{k}^{+}$ $\displaystyle=$
$\displaystyle\left\\{l\in\left\\{1,\ldots,d\right\\}:x_{li\left(k\right)}-x_{lj\left(k\right)}\geq
0\right\\}$ $\displaystyle\mathcal{L}_{k}^{-}$ $\displaystyle=$
$\displaystyle\left\\{l\in\left\\{1,\ldots,d\right\\}:x_{li\left(k\right)}-x_{lj\left(k\right)}<0\right\\}.$

######  Theorem 3.1.

$\mathbf{U}=\left(U_{1},\ldots,U_{d}\right)$ in representation (36)
conditioning on $K=k$ has cumulative distribution function
$F_{\mathbf{U}|K}\left(u_{1},\ldots,u_{d};k\right)$ given by

$\displaystyle\mathbb{P}\left(U_{1}\leq u_{1},\ldots,U_{d}\leq
u_{d}\left|K=k\right.\right)=F_{\mathbf{V},\mathcal{L}^{-}_{k}}\left(\mathbf{v}_{k}\right)$

where $\mathbf{v}_{k}=(v_{1,k},\ldots,v_{d,k})$ with
$v_{l,k}=F_{U_{l}\left|K\right.}\left(u_{l};k\right)$.

The following result comes from summing over all possible values of $K$.

######  Corollary 2.

The random vector $\left(U_{1},\ldots,U_{d}\right)$ in representation (36),
has distribution $F_{\mathbf{U}}\left(u_{1},\ldots,u_{d}\right)$ given by

$\displaystyle\mathbb{P}\left(U_{1}\leq u_{1},\ldots,U_{d}\leq
u_{d}\right)=\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\sum^{\left|\mathcal{E}\right|}_{k=1}F_{\mathbf{V},\mathcal{L}^{-}_{k}}\left(\mathbf{v}_{k}\right).$

For random vectors in representation (10) we are able to derive a closed form
expression of the conditional distribution of $\mathbf{U}$.

######  Corollary 3.

$\mathbf{U}=\left(U_{1},\ldots,U_{d}\right)$ in representation (10)
conditioning on $K=k$ has cumulative distribution function

$\displaystyle
F_{\mathbf{U}|K}\left(u_{1},\ldots,u_{d};k\right)=\max\left(v^{+}_{k}+v^{-}_{k}-1,0\right)$

where $v^{+}_{k}=\min\\{v_{l,k},l\in\mathcal{L}_{k}^{+}\\}$,
$v^{-}_{k}=\min\\{v_{l,k},l\in\mathcal{L}_{k}^{-}\\}$ and
$v_{l,k}=F_{U_{l}\left|K\right.}\left(u_{l};k\right)$. The pairs of variables
$U_{l}$ and $U_{l^{\prime}}$ $l\neq l^{\prime}$ have cumulative distribution
function:

(45)

From Corollary 3 and following the definition of Fréchet bound (Fréchet,,
1951), the elements $U_{l}$ and $U_{l^{\prime}}$ of the stochastic
representation (10) are monotonic in the same direction if
$l,l^{\prime}\in\mathcal{L}^{+}_{k}$ or $l,l^{\prime}\in\mathcal{L}^{-}_{k}$
and antithetic if $l\in\mathcal{L}^{-}_{k},l^{\prime}\in\mathcal{L}^{+}_{k}$
or $l\in\mathcal{L}^{+}_{k},l^{\prime}\in\mathcal{L}^{-}_{k}$, conditionally
on living on the $k$-th segment.

###  3.2 Multivariate Kendall’ $\tau$ and Spearman’s $\rho$

Multivariate Kendall’ $\tau$ and Spearman’s $\rho$ are multivariate measures
of concordance introduced in Joe, (1990) as a generalization of the well known
bivariate measures. Fuchs et al., (2018); Ahn and Fuchs, (2020) show that
$d$-CTM vectors have the same minimal multivariate Kendall’s $\tau$ but
different Spearman’s $\rho$. Then, Spearman’s $\rho$ can be used to rank
$d$-CTM vectors in concordance order.

Let $\mathbf{U}$ and $\mathbf{W}$ be independent random vectors with the same
distribution $F_{\mathbf{U}}=F_{\mathbf{W}}=G$. The multivariate Kendall’s
$\tau$ is defined as:

$\displaystyle\tau\left(F_{\mathbf{U}}\right)$ $\displaystyle=$
$\displaystyle\dfrac{2^{d}}{2^{d-1}-1}\left[\int_{\left[0,1\right]^{d}}F_{\mathbf{U}}\left(\mathbf{u}\right)dF_{\mathbf{U}}\left(\mathbf{u}\right)-\dfrac{1}{2^{d}}\right]$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}}{2^{d-1}-1}\left[\mathbb{E}\left[G\left(\mathbf{W}\right)\right]-\dfrac{1}{2^{d}}\right]$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}}{2^{d-1}-1}\left[\mathbb{P}\left(\mathbf{U}\leq\mathbf{W}\right)-\dfrac{1}{2^{d}}\right].$

Fuchs et al., (2018) show that if $\mathbf{U}$ and $\mathbf{W}$ are $d$-CTM
then $\tau\left(F_{\mathbf{U}}\right)$ attains its minimal value

$\displaystyle\tau\left(F_{\mathbf{U}}\right)=\tau_{\min}=-\dfrac{1}{2^{d-1}-1}.$
(46)

We provide an analytical expression of the Kendall’s $\tau$ for vectors with
the stochastic representation of Eq. 36.

######  Proposition 3.1.

Let $\mathcal{S}=(\mathcal{G},\mathbf{X})$ be the segment set with
$\mathcal{G}=\left\\{\mathcal{V},\mathcal{E}\right\\}$. Let $\mathbf{U}$ and
$\mathbf{W}$ two independent copies of a random vectors with the generalized
line segments representation (36) on $\mathcal{S}$ , then:

$\displaystyle\tau\left(F_{\mathbf{U}}\right)$ $\displaystyle=$
$\displaystyle\dfrac{2^{d}}{2^{d-1}-1}\left[\dfrac{1}{\left|\mathcal{E}\right|}\sum^{\left|\mathcal{E}\right|}_{k_{\mathbf{U}}=1}\mathbb{P}\left(R_{\mathbf{V},\mathcal{L}^{-}_{k_{\mathbf{U}}}}\left(\mathbf{V}\right)\leq\mathbf{Y}_{k_{\mathbf{U}}}\right)-\dfrac{1}{2^{d}}\right],$
(47)

with
$\mathbf{Y}_{k_{\mathbf{U}}}=\left(F_{U_{1}\left|K_{\mathbf{U}}\right.}\left(W_{1};k_{\mathbf{U}}\right),\ldots,F_{U_{d}\left|K_{\mathbf{U}}\right.}\left(W_{d};k_{\mathbf{U}}\right)\right)$.

We use the definition of Spearman’s $\rho$ in Joe, (1990)888The Spearman’s
$\rho$ following the definition in Ahn and Fuchs, (2020) is obtained as:
$\frac{\rho\left(F_{U}\right)+\rho\left(F_{W}\right)}{2}$, with
$\mathbf{W}=R_{\mathcal{D},\frac{1}{2}}\left(\mathbf{U}\right)$. for
distributions with standard uniform marginals. Let $F\left(\mathbf{U}\right)$
be the cumulative distribution function of the random vector $\mathbf{U}$;
then the multivariate Spearman’s $\rho$ is

$\displaystyle\rho\left(F_{\mathbf{U}}\right)$ $\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\int_{\left[0,1\right]^{d}}F_{\mathbf{U}}\left(\mathbf{u}\right)d\mathbf{u}-\dfrac{1}{2^{d}}\right)$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\int_{\left[0,1\right]^{d}}\prod^{d}_{l=1}u_{L}dF_{\mathbf{U}}\left(\mathbf{u}\right)-\dfrac{1}{2^{d}}\right)$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\mathbb{E}\left[\prod^{d}_{l=1}U_{l}\right]-\dfrac{1}{2^{d}}\right).$

The attainable lower bound $\rho_{\min}$ can be computed using the lower bound
for $\mathbb{E}\left[\prod^{d}_{l=1}U_{l}\right]$ given in Corollary 4.1 of
Wang and Wang, (2011). We report the values in Table 1.

Table 1: Minimum values of multivariate Spearman’s $\rho$, using the lower bound in Corollary 4.1 of Wang and Wang, (2011) $d$ | 2 | 3 | 4 | 5 | 10 | 20 | 50 | 100
---|---|---|---|---|---|---|---|---
$\rho_{\min}$ | -1 | -0.56 | -0.32 | -0.18 | -0.01 | -1.99$\cdot 10^{-5}$ | -4.53$\cdot 10^{-14}$ | -7.97$\cdot 10^{-29}$

######  Proposition 3.2.

Let $\mathbf{U}$ have the generalized representation (36), $\mathbf{V}$ be
reflection invariant, and $\xi^{*}\leq 1$. Then:

  * (a)

$\displaystyle\rho\left(F_{\mathbf{U}}\right)$ $\displaystyle=$

$\displaystyle\sum^{d}_{m=0}\sum_{\tiny\begin{array}[]{c}\mathcal{L}_{m}\subseteq\mathcal{D}\\\
\left|\mathcal{L}_{m}\right|=m\end{array}}\xi_{\mathcal{L}_{m}}\rho\left(F_{\mathbf{V},\mathcal{D}\setminus\mathcal{L}_{m}}\right)+\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\dfrac{1}{2^{d}}\left(\xi^{*}-1\right),$

(50)

with

$\displaystyle\xi_{\mathcal{L}_{m}}$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\sum^{\left|\mathcal{E}\right|}_{k=1}\left(\prod_{l\in\mathcal{L}_{m}}x_{l,i\left(k\right)}\prod_{l\in\mathcal{D}\setminus\mathcal{L}_{m}}x_{l,j\left(k\right)}\right)$
$\displaystyle\xi^{*}$ $\displaystyle=$
$\displaystyle\sum^{d}_{m=0}\sum_{\tiny\begin{array}[]{c}\mathcal{L}_{m}\subseteq\mathcal{D}\\\
\left|\mathcal{L}_{m}\right|=m\end{array}}\xi_{\mathcal{L}_{m}}=\dfrac{1}{\left|\mathcal{E}\right|}\sum^{\left|\mathcal{E}\right|}_{k=1}\prod^{d}_{l=1}\left(x_{l,i\left(k\right)}+x_{l,j\left(k\right)}\right);$

  * (b)

$\displaystyle\rho\left(F_{\mathbf{U}}\right)$ $\displaystyle=$
$\displaystyle\xi^{*}\rho\left(F_{\mathbf{V}}\right)+\left(1-\xi^{*}\right)\left(-\dfrac{\left(d+1\right)}{2^{d}-\left(d+1\right)}\right);$
(52)

  * (c)

$\rho\left(F_{\mathbf{U}}\right)\leq\rho\left(F_{\mathbf{V}}\right)$.

We provide a simplified formula for
$\rho\left(F_{\mathbf{V},\mathcal{D}\setminus\mathcal{L}_{m}}\right)$ in the
case of two line segment constructions which will be studied later on in this
paper. For the line segment representation (10) we have:

$\displaystyle\rho\left(F_{\mathbf{V},\mathcal{D}\setminus\mathcal{L}_{m}}\right)$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\mathbb{E}\left[V^{m}\left(1-V\right)^{d-m}\right]-\dfrac{1}{2^{d}}\right)$
(53) $\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(B\left(m+1,d-m+1\right)-\dfrac{1}{2^{d}}\right)$
(54)

where $B\left(x,y\right)$ is the Euler’s beta function.

For the representation (36), under the independence assumption $V_{l}$
$l=1,\ldots,d$ iid, the Spearman’s $\rho$ is:

$\displaystyle\rho\left(F_{\mathbf{V},\mathcal{D}\setminus\mathcal{L}_{m}}\right)$
$\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\prod_{l\in\mathcal{L}_{m}}\mathbb{E}\left[V_{l}\right]\prod_{l\in\mathcal{D}\setminus\mathcal{L}_{m}}\mathbb{E}\left[\left(1-V_{l}\right)\right]-\dfrac{1}{2^{d}}\right)$
(55) $\displaystyle=$
$\displaystyle\dfrac{2^{d}\left(d+1\right)}{2^{d}-\left(d+1\right)}\left(\dfrac{1}{2^{d}}-\dfrac{1}{2^{d}}\right)=0.$

Finally note that the multivariate Kendall’s $\tau$ and Spearman’s $\rho$ are
permutation invariant, thus permuted uniform vectors will have the same
measures.

##  4 Special cases

In this section, we discuss several examples starting from the new
constructions proposed in this paper and then reviewing the constructions
proposed in the literature which are special cases of our stochastic
representations (10) or (36). We use multivariate Spearman’s $\rho$ to rank
the proposals in concordance order999To keep the paper to a reasonable length,
we report here the values of the multivariate Spearman’s $\rho$ for dimensions
$d$ from 2 to 5. We have evaluated the ordering of the constructions up to
dimension $d=20$, and no changes were observed..

###  4.1 Circulant Variates

Obtaining the coordinate matrix $\mathbf{X}$ used in Eq. 10 can be costly in
high dimensions especially when numerical procedures are used to solve the
optimization problem stated in Section 2.2. We propose suitable constraints on
the segment set $\mathcal{S}=(\mathbf{X},\mathcal{G})$ to reduce the
computational cost of our procedure. The proposed conditions on the coordinate
matrix $\mathbf{X}$ allows for decoupling the CTM constraint.

First, we assume that the number of vertexes $n$ is equal to the dimension $d$
of the random vector and $\mathbf{X}=d/2\tilde{\mathbf{X}}$ where
$\tilde{\mathbf{X}}$ is doubly stochastic, and obtain:

$\displaystyle\sum^{d}_{i=1}x_{i,k}$ $\displaystyle=$
$\displaystyle\sum^{n}_{k=1}x_{i,k}=\dfrac{d}{2}.$

This assumption allows us to simplify the optimization problem and to search
for independent solutions for each row of $\mathbf{X}$.

We assume further constraints on the matrix $\mathbf{X}$ and on the graph
$\mathcal{G}$ such that the same optimization problem is solved for all rows
of $\mathbf{X}$. We assume the first coordinates of the $d$ vertexes are
arranged in increasing order $x_{11}\leq\ldots\leq x_{1d}$ and compute the
$k$-th coordinates as the $k$-th circular permutation of the first ones:

$\displaystyle\begin{array}[]{ccc}x_{k1}&=&x_{1\,\left(k-1\right)\left(\mathrm{mod}\,d\right)+1}\\\
&\vdots&\\\
x_{ki}&=&x_{1\,\left(i-1+\left(k-1\right)\right)\left(\mathrm{mod}\,d\right)+1}\\\
&\vdots&\\\
x_{kd}&=&x_{1\,\left(d-1+\left(k-1\right)\right)\left(\mathrm{mod}\,d\right)+1}.\end{array}$
(61)

The resulting coordinate matrix $\mathbf{X}$ is a circulant matrix with $i$-th
row sum equal to the $i$-th column sum for all rows. Imposing
$x_{11}\leq\ldots\leq x_{1d}$ implies the same set of $\mathbf{a}_{l}$ is used
for all $l\in\mathcal{D}$, with the same multiplicities
$\left|\mathcal{M}_{l,m}\right|=\left|\mathcal{M}_{1,m}\right|$ but with
different positions $m\in\mathcal{M}_{l,m}$ as effect of the circular
permutation.

Furthermore, we choose the edge set in such way to have the same projected
graph for each set of coordinates, and assume a circulant graph which is
invariant by circular shifts of the vertexes.

######  Definition 4.1 (Circulant Graph).

Given a subset
$\mathcal{L}\subseteq\left\\{1,\ldots,\left\lfloor\frac{d}{2}\right\rfloor\right\\}$
then the $d$-vertex circulant graph $\mathcal{C}_{d}\left(\mathcal{L}\right)$
is a graph with vertexes $1,\ldots,d$ and edge set
$\mathcal{E}_{d,\mathcal{L}}$ is such that
$(i,j)\in\mathcal{E}_{d,\mathcal{L}}$ if either
$\left|i-j\right|\in\mathcal{L}$ or
$\left(d-\left|i-j\right|\right)\in\mathcal{L}$.

Label |  $\mathcal{G}$ |  $\mathbf{x}_{1}$ |  Spearman’s $\rho$
---|---|---|---
$\mathcal{C}_{2}\left(\left\\{1\right\\}\right)$ |  $1$$2$ |  $\left(0,1\right)$ |  $-1$
$\mathcal{C}_{3}\left(\left\\{1\right\\}\right)$ |  $1$$2$$3$ |  $\left(0,\frac{1}{2},1\right)$ |  -0.5
$\mathcal{C}_{4}\left(\left\\{1\right\\}\right)$ |  $1$$2$$3$$4$ |  $\left(0,\frac{1}{3},\frac{2}{3},1\right)$ |  -0.2840
$\mathcal{C}_{4}\left(\left\\{1,2\right\\}\right)$ |  $1$$2$$3$$4$ |  $\left(0,\frac{1}{2}-\frac{1}{2\sqrt{5}},\frac{1}{2}+\frac{1}{2\sqrt{5}},1\right)$ |  -0.2763
$\mathcal{C}_{4}\left(\left\\{2\right\\}\right)$ |  $1$$2$$3$$4$ |  $\left(0,0,1,1\right)$ |  -0.2121
$\mathcal{C}_{5}\left(\left\\{1\right\\}\right)$ |  $1$$2$$3$$4$$5$ |  $\left(0,1/4,2/4,3/4,1\right)$ |  -0.1659
$\mathcal{C}_{5}\left(\left\\{1,2\right\\}\right)$ |  $1$$2$$3$$4$$5$ |  $\left(0,\frac{1}{2}-\frac{\sqrt{3}}{2\sqrt{7}},\frac{1}{2},\frac{1}{2}+\frac{\sqrt{3}}{2\sqrt{7}}),1\right)$ |  -0.1577
$\mathcal{C}_{5}\left(\left\\{2\right\\}\right)$ |  $1$$2$$3$$4$$5$ |  $\left(0,0,\frac{1}{2},1,1\right)$ |  -0.1385
Table 2: Circulant graphs $\mathcal{C}_{d}\left(\mathcal{L}\right)$ of the
Circulant Countermonotonic on Segments up to dimension 5.

The circular symmetry imposed on the vertex coordinates and on the graph
simplifies the optimization problem. Whatever the multiplicities
$\left|\mathcal{M}_{1,m}\right|$, under Assumptions 2 and 1, we obtain the
following results:

$\displaystyle\Psi_{l}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle\Psi_{1}\left(\mathbf{a}_{1}\right)=\Phi_{1}\left(\mathbf{x}_{1}\right)=-\dfrac{1}{2\left|\mathcal{E}_{d,\mathcal{L}}\right|}\displaystyle\sum_{\left(i,j\right)\in\mathcal{E}_{d,\mathcal{L}}}\log\left|x_{1,i}-x_{1,j}\right|,$
$\displaystyle\sum^{n_{1}}_{m=1}\left|\mathcal{M}_{1,k}\right|a_{1,m}=\sum^{d}_{i=1}x_{1,i}.$

This rewriting of constraint and objective function allows one to minimize on
$\mathbf{x}_{1}$ instead of $\mathbf{a}_{l}$. Minimization automatically
excludes vectors that violate Assumption 1 given
$\mathcal{E}_{d,\mathcal{L}}$, because for those cases the objective function
is infinite. The following definitions introduce formally our new proposal
called Circulant Variates (CCV).

######  Definition 4.2 (Circulant Countermonotonic).

We call the circulant matrix $\mathbf{X}\in\mathbb{R}^{d}\times\mathbb{R}^{n}$
whose rows are obtained by the $d$ circular shifts of the first row, a
solution of the Circulant Countermonotonic on Segments problem, on the
circulant graph $\mathcal{C}_{d}\left(\mathcal{L}\right)$, if $x_{1,1}=0$,
$x_{1,d}=1$ and $\left\\{x_{1,2},\ldots,x_{1,d-1}\right\\}$ solve the convex
minimization problem:

$\displaystyle\displaystyle\min_{\left\\{x_{l,m}\right\\}^{d-1}_{m=2}\in\left[0,1\right]^{d-2}}\Phi_{1}\left(\mathbf{x}_{1}\right),$

with

$\displaystyle\Phi_{1}\left(\mathbf{x}_{1}\right)=-\dfrac{1}{2\left|\mathcal{E}_{d,\mathcal{L}}\right|}\displaystyle\sum_{\left(i,j\right)\in\mathcal{E}_{d,\mathcal{L}}}\log\left|x_{1,i}-x_{1,j}\right|,$

subject to:

$\displaystyle\sum^{d}_{i=1}x_{1,i}=\dfrac{d}{2}.$

######  Definition 4.3 (Circulant Variates).

Let
$\mathcal{S}_{d,\mathcal{L}}=\left\\{\mathbf{X},C_{d}\left(\mathcal{L}\right)\right\\}$
a collection of segments such that $\mathbf{X}$ is a solution Circulant
Countermonotonic on Segments problem on
$\mathcal{C}_{d}\left(\mathcal{L}\right)$. Variates obtained from the
components of the $d$-dimensional random vector uniformly distributed on
$\mathcal{S}_{d,\mathcal{L}}$ are Circulant Variates (CCV).

The following corollary is the application of Theorem 2.1 to CCV.

######  Corollary 4.

CCV are marginally standard uniform and constant in sum.

For $C_{d}\left(\left\\{1\right\\}\right)$ a solution can be derived as
follows. The sum in each of the first $d-1$ equations of (10) has only two
terms, and considering the $m$-th equation we have

$\displaystyle\dfrac{1}{\left(x_{1m+1}-x_{1m}\right)}+\dfrac{1}{\left(x_{1d}-x_{11}\right)}$
$\displaystyle=$ $\displaystyle d.$

Substituting the constraints $x_{11}=0$ and $x_{1d}=1$ , we find that the
$x$’s are uniformly spaced on the unit interval and satisfy all the $d-2$
equations in (21). The case of $C_{3}\left(\left\\{1\right\\}\right)$ was
already studied in Nelsen and Úbeda-Flores, (2012) . In their example, the
probability mass has a distribution uniform on the edges of the triangle with
vertexes $\mathbf{x}_{1}=\left(0,1/2,1\right)$,
$\mathbf{x}_{2}=\left(1/2,1,0\right)$, $\mathbf{x}_{3}=\left(1,0,1/2\right)$.
The other two vertexes are different 3-cycles of the first one. This implies
that the row sum of
$\mathbf{x}=(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})$, is equal to its
column sum and both are equal to $3/2$. Lee and Ahn, (2014) show that the
construction in Nelsen and Úbeda-Flores, (2012) is $3$-CTM. Additionally, the
exchangeable version of the $C_{d}\left(\left\\{1\right\\}\right)$
construction is distributionally equivalent to the degenerate random balanced
sampling introduced in Eq. (8) of Gerow and Holbrook, (1996). They propose to
generate $Z_{1}$ as an uniform random variable on $\left[-1,1\right]$ and
obtain the remaining variables according to

$\displaystyle Z_{l}$ $\displaystyle=$ $\displaystyle
c_{l}-\dfrac{Z_{1}}{d-1},\quad c_{l}=-1+\dfrac{2l-3}{d-1}$

$l\in\mathcal{D}$, and then randomly permute the $Z_{l}$. They show that
$\sum^{d}_{l=1}Z_{l}=0$ and that once permuted, the $Z_{l}$’s are uniformly
distributed on $\left[-1,1\right]$. If we set $U_{l}=(Z_{l}+1)/2$ then the
permuted version can be written in terms of permuted
$C_{d}\left(\left\\{1\right\\}\right)$ construction:

######  Proposition 4.1.

The exchangeable version of $U_{l}=(Z_{l}+1)/2$ $l\in\mathcal{D}$ has the same
distribution of the exchangeable version of CCV with dependence graph
$C_{d}\left(\left\\{1\right\\}\right)$.

###  4.2 Rotation Sampling

Fishman and Huang, (1983) rephrase the original Hammersley and Morton, (1956)
proposal for $d>2$ obtaining an equivalent construction with the standard
marginal uniformity that was missing in Hammersley and Morton, (1956). Their
construction was named rotation sampling because the modulo one arithmetic on
which it is based is often associated with circular motion.

######  Proposition 4.2.

The line segment stochastic representation (10) of the rotation sampling in
the Example 1.1 has $2d$ vertexes with coordinate matrix of elements:

$\displaystyle x_{l,m}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ccc}\dfrac{l+m-1}{d}&\hbox{if}&m<d+2-l\\\
\dfrac{l+m-1-d}{d}&\hbox{if}&m\geq d+2-l\end{array}\right.$ $\displaystyle
x_{l,d+m}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ccc}\dfrac{l+m-2}{d}&\hbox{if}&m<d+2-l\\\
\dfrac{l+m-2-d}{d}&\hbox{if}&m\geq d+2-l\end{array}\right.$

and edge set:

$\displaystyle\mathcal{E}^{RS}=\left\\{\left(i,j\right)\in\\{1,2,\ldots,d\\}\times\\{d+1,d+2,\ldots,2d\\}|j=d+i\right\\}.$

######  Corollary 5.

A rotation sampling random vector has standard uniform marginals and is not
$d$-CTM.

We report in Table 3 the analytic values of the Multivariate Spearman’s $\rho$
for rotation sampling vectors up to dimension $d=5$.

Table 3: Values of the multivariate Spearman’s $\rho$ for rotation sampling. $d$ | 2 | 3 | 4 | 5
---|---|---|---|---
$\rho_{\min}$ | -0.5 | -0.33 | -0.2168 | -0.1372

For $d=2$ this construction does not reduce to the usual antithetic variates.
This is suggested by a value of multivariate Spearman’s $\rho$ different from
the value of $-1$ attained by the Fréchet lower bound. Lacking the constant
sum property, the proposal has multivariate Spearman’s $\rho$ larger than the
one of the $d$-CTM proposals considered in this paper.

###  4.3 Arvidsen and Johnson: A Fresh Look

In the pioneering paper of Arvidsen and Johnsson, (1982), the objective of
variance reduction is obtained by designing the first standard uniform $d$-CTM
construction (1.4). Craiu and Meng, (2005) have shown that this construction
is displacing the binary digits of $U_{1}$ and give the name permuted
displacement to its exchangeable version. We show the following relationship
between the two.

######  Proposition 4.3.

For $d=3$, the construction of Arvidsen and Johnsson, (1982) given in Example
1.4 is equivalent to the antithetic proposal of Gaffke and Rüschendorf, (1981)
given in (1.2).

$1$$b^{d-2}+1$$2b^{d-2}$$b^{d-2}$$\vdots$ Figure 5: $\mathcal{G}$ for
$b$-based Arvidsen and Johnson construction

The Ardvisen and Johnson’ construction and the following family of general
constructions admit a line segment representation. The generalization is
useful to show that Ardvisen and Johnson is the only $d$-CTM within this
family and attains consequently the minimal Spearman’s $\rho$ within this
family.

######  Definition 4.4.

Given $U_{1}\sim\mathcal{U}[0,1]$ and $b\in\mathbb{N}$ the base-$b$ Ardvisen
and Johnson construction is:

$\displaystyle\begin{array}[]{ccccc}U_{i}=\left(b^{i-2}U_{1}+1/b\right)\mathrm{mod}\,1&\hbox{for}&i=\left\\{2,\ldots,d-1\right\\}&&U_{d}=1-\left(b^{d-2}U_{1}\right)\mathrm{mod}\,1\end{array}$
(65)

######  Proposition 4.4.

The line segment stochastic representation (10) of the construction in (65)
has $2b^{d-2}$ vertexes with coordinate matrix
$\mathbf{X}=(\mathbf{z}^{T},\mathbf{y}^{T})$ where:

$\displaystyle\begin{array}[]{ccccccccc}\mathbf{y}_{1}&=&\left(0,\dfrac{1}{b^{d-2}},\ldots,\dfrac{b^{d-2}-1}{b^{d-2}}\right)^{T},&\mathbf{z}_{1}&=&\mathbf{y}_{1}+\dfrac{1}{b^{d-2}}\mathbf{1}_{b^{d-1}}\\\
\mathbf{y}_{k}&=&\left(b^{k-2}\mathbf{y}_{1}+\dfrac{1}{b}\right)\mathrm{mod}\,1,&\mathbf{z}_{k}&=&\mathbf{y}_{k}+\dfrac{1}{b^{d-k}}\mathbf{1}_{b^{d-1}}\\\
\mathbf{y}_{d}&=&\mathbf{1}_{b^{d-1}},&\mathbf{z}_{d}&=&\mathbf{y}_{d}-\mathbf{1}_{b^{d-1}}\\\
\mathbf{y}&=&\left(\mathbf{y}_{1},\ldots,\mathbf{y}_{d}\right)&\mathbf{z}&=&\left(\mathbf{z}_{1},\ldots,\mathbf{z}_{d}\right)\\\
\end{array}$

and edge set:

$\displaystyle\mathcal{E}^{AJ}=\left\\{\left(i,j\right)\in\left\\{1,\ldots,b^{d-2}\right\\}\times\left\\{1,\ldots,b^{d-2}\right\\}|j=b^{d-2}+i\right\\}$

######  Corollary 6.

The base-$b$ Ardvisen and Johnson random vector has standard uniform
marginals. Only the case $b=2$ is $d$-CTM.

Table 4 shows that the constructions in Definition 4.4 attain the lowest
Spearman’s $\rho$ for $b=2$.

Spearman’s $\rho$ | $b$ | GR
---|---|---
| 1 | 2 | 3 | 4 | 5 |
$d$ | 2 | -1 | -1 | -1 | -1 | -1 | -1
| 3 | -0.3333 | -0.5000 | -0.3333 | -0.2083 | -0.1200 | -0.5000
| 4 | -0.0909 | -0.2822 | -0.1662 | -0.0869 | -0.0367 | -0.2525
| 5 | 0.0154 | -0.1637 | -0.0933 | -0.0455 | -0.0165 | -0.1538
Table 4: Multivariate Spearman’s $\rho$ for base-$b$ Ardvisen and Johnson
random vectors and the multivariate proposal of Gaffke and Rüschendorf, (1981)
in Example 1.3

Table 4 reports also the multivariate proposal of Gaffke and Rüschendorf,
(1981) (GR) described in Example 1.3. The table shows that GR performs better
than the non $d$-CTM proposal, but worse than the AJ proposal because it has
independence as one of the main ingredients as discussed in the introduction.
This also show the effectiveness of multivariate $\rho$ as a ranking measure
for $d$-CTM vectors.

###  4.4 Latin Hypercube Iterations

In this section, we will reconsider the Iterated Latin Hypercube construction
introduced in Craiu and Meng, (2005) establish a relationship with a
$d$-dimensional generalization of the superstar introduced in Gerow and
Holbrook, (1996), and with our strict countermonotonic on segments
construction.
The latin hypercube sampling, introduced by McKay et al., (1979), and further
developed by Craiu and Meng, (2005) with the goal of obtaining variance
reduction in MCMC sampling, consists in the following steps: for
$t=0,\ldots,T$ take an iid standard uniform $d$-dimensional random vector
$\mathbf{U}_{0}$ and let
$\mathcal{D}^{\pi}_{t}=\left(\pi_{t}\left(0\right),\ldots,\pi_{t}\left(d-1\right)\right)^{T}$
be a permutation of $\\{0,1,\ldots,d-1\\}$ independent of
$\mathbf{U}_{0},\ldots,\mathbf{U}_{t-1}$ and

$\displaystyle\mathbf{U}_{t}=\dfrac{1}{d}\left(\mathcal{D}^{\pi}_{t}+\mathbf{U}_{t-1}\right)$
(67)

If $t=1$ (67) corresponds to the original Latin Hypercube Sampling, and $t>1$
to the Iterated Latin Hypercube procedure introduced in Craiu and Meng,
(2005). It was showed in Craiu and Meng, (2006) that ILH iterations represent
an Iterated Function System with probabilities (IFSP)
$\left(\left[0,1\right]^{d},\left(w_{\pi}\right),p_{\pi}\right)$ with
$w_{\pi}$ similitudes with contraction ratio $\dfrac{1}{d}$ associated to each
permutation $\pi$ of $\left\\{1,\ldots,d\right\\}$

$\displaystyle
w_{\pi}\left(\mathbf{u}\right)=\left(\dfrac{\pi\left(1\right)-1}{d}+\dfrac{u_{1}}{d},\ldots,\dfrac{\pi\left(d\right)-1}{d}+\dfrac{u_{d}}{d}\right).$
(68)

We can show that the $3$-dimensional superstar considered in Gerow and
Holbrook, (1996) can be obtained by using the same IFSP:

######  Proposition 4.5.

The $3$-dimensional superstar proposed in Gerow and Holbrook, (1996):

$\displaystyle\mathbf{X}_{t}=f_{k}\left(\mathbf{X}_{t-1}\right)=\dfrac{1}{3}\mathbf{X}_{t-1}+\dfrac{2}{3}\mathbf{V}_{k}$
(69)

with $\mathbf{V}_{k}$ a random permutation of $\left\\{-1,0,1\right\\}$ and an
initial $\mathbf{X}_{0}$ a $3$-dimensional vector such that
$-1\leq\mathbf{X}_{i0}\leq 1$ and $\sum^{3}_{i=1}\mathbf{X}_{i0}=0$, up to a
change of support, is generated by the same IFSP of the 3-dimensional version
of the Iterated Latin Hypercube construction introduced in Craiu and Meng,
(2005).

######  Proposition 4.6.

The line segment stochastic representation (36) of the construction in (67)
has $2d!$ vertexes with coordinates
$x_{l,k}=\frac{\pi_{k}\left(l\right)+1}{d}$ and
$x_{l,d!+k}=\frac{\pi_{k}\left(l\right)}{d}$, $k=1,\ldots,d!$ and edge set:

$\displaystyle\mathcal{E}^{LH}=\left\\{\left(i,j\right)\in\left\\{1,\ldots,d!\right\\}\times\left\\{d!+1,\ldots,2d!\right\\}|j=d!+i\right\\}$

Craiu and Meng, (2005) obtained standard marginal uniformity for the special
case of an initial vector with iid components. In the superstar case, the
original distribution is concentrated on points but converges to standard
marginal uniforms as showed in Gerow and Holbrook, (1996). In Craiu and Meng,
(2005), it is also shown that in the limit $t\rightarrow\infty$ the ILH is
$d$-CTM . We show that the ILH iterations preserve marginal uniformity and
constant sum in the general case.

######  Corollary 7.

Let $\mathbf{U}_{t-1}$ be a dependent random vector of dimension $d$, whose
coordinates add up to d/2 (a.s.) and each coordinate has a $\mathcal{U}[0,1]$
distribution. Then the random vector $\mathbf{U}_{t}$ in (67) has all its
coordinates marginally $\mathcal{U}[0,1]$ distributed and adding up to $d/2$.

The preservation of strict $d$-CTM property raises the possibility of using
ILH iterations on Arvidsen and Johnson’s and CCV constructions. In addition,
using results in section 3.2, we can give a closed-form expression for the
multivariate Kendall’s $\tau$ and Spearman’s $\rho$ for ILH iteration applied
to $\mathbf{V}$ with i.i.d components or in representation (10) and rank the
obtained random vector in the correlation order. In particular, ILH iterations
on the Arvisen and Johnson construction and CCV have a constant sum and
minimal multivariate Kendall’s $\tau$ ( equation (46)). Spearman’s $\rho$ for
those constructions can be easily obtained using the deterministic composition
and equation (53). Using the segment representation in Proposition 4.6, the
ILH($T$) proposal applied to $\mathbf{V}$ with i.i.d components, can be
expressed as a $T$-fold deterministic composition. Each composition expands
the cardinality of vertex and edge sets by $d!$, resulting in a vertex set of
cardinality $2\left(d!\right)^{T}$ and an edge set of cardinality
$\left(d!\right)^{T}$ . To maintain a feasible notation we substitute the
index $k$ of the different edges with the multi-index
$\left\\{k_{t}\right\\}^{T}_{t=1}$ where each $k_{t}=1\ldots d!$ :

$\displaystyle\alpha_{l,k_{1},\ldots,k_{T}}=\sum^{T}_{r=1}\dfrac{\pi_{k_{r}}\left(l\right)}{d^{T-r+1}},$
$\displaystyle\beta_{l,k_{1},\ldots,k_{T}}=\sum^{T}_{r=1}\dfrac{\pi_{k_{r}}\left(l\right)}{d^{T-r+1}}+\dfrac{1}{d^{T}}$
(70)

The following proposition allows for computing multivariate Kendall’s $\tau$
and Spearman’s $\rho$ for the ILH($T$) case.

######  Proposition 4.7.

Let $U_{0}$ a $d$ dimensional random vector with iid $\mathcal{U}[0,1]$
components and let $\mathbf{U}_{T}$ be the vector obtained by applying the
transformation $T$ given in (67) to it. Then the Kendall’s $\tau$ is

$\displaystyle\tau$ $\displaystyle=$
$\displaystyle\dfrac{1}{2^{d-1}-1}\left(\dfrac{1}{\left(d!\right)^{T}}-1\right)$

and the Spearmans’ $\rho$ has coefficient $\xi^{*}$ of Equation (52) equal to

$\displaystyle\xi^{*}$ $\displaystyle=$

$\displaystyle\sum_{m_{0}+m_{1}+\ldots+m_{T}=d}{d\choose
m_{0},m_{1},\dots,m_{T}}\prod^{T}_{t=1}\left\\{{d\choose
m_{t}}^{-1}\sum_{{\tiny\left\\{i_{1},\ldots,i_{m_{t}}\right\\}\in\mathcal{P}_{\mathcal{D},m_{t}}}}\prod^{m_{t}}_{l=1}\dfrac{2\left(i_{l}-1\right)}{d^{T-t+1}}\right\\}\left(\dfrac{1}{d^{T}}\right)^{m_{0}}$

where $\mathcal{C}_{\mathcal{D},m_{2}}$ denotes the set of combinations of the
elements of $\mathcal{D}$ with $m_{t}$ elements.

In the one iteration case, i.e., $T=1$, we obtain the following expression:

$\displaystyle\xi^{*}$ $\displaystyle=$
$\displaystyle\prod^{d-1}_{l=0}\dfrac{(2l+1)}{d}=\prod^{d}_{l=1}\dfrac{\left(2l-1\right)}{d}=\dfrac{\left(2d-1\right)!}{d^{d}\left(2^{d-1}\left(d-1\right)!\right)}=\dfrac{\prod^{2d-1}_{l=d}l}{d^{d}2^{d-1}}$

and using the arithmetic and geometric means inequality we obtain the bound:

$\displaystyle\xi^{*}_{1}$ $\displaystyle\leq$
$\displaystyle\dfrac{\left(\frac{1}{d}\sum^{2d-1}_{l=d}l\right)^{d}}{d^{d}2^{d-1}}\leq
2\left(\dfrac{3d-1}{4d}\right)^{d}\leq 1.$

For the sake of comparison, in figure 6 we summarize Kendall’s $\tau$ and
Spearman’s $\rho$ for all the constructions discussed in this section. Since
Kendall’s $\tau$ and Spearman’s $\rho$ are permutation invariant, the same
ranking applies to the exchangeable versions of the constructions.

Figure 6: Ratio of minimum value attainable by multivariate association
measures and the values attained by different antithetic random vectors
(vertical axis) as a function of dimension $d$ (horizontal axis). Left: the
Spearman’s $\rho$ measure. Right: the Kendall’s $\tau$ measure. Note: in each
plot, larger values indicate constructions farther from the minimum.

##  5 A Central Limit Theorem

In this section, we study the Central Limit Theorem for our best performing
classes of variates.

######  Definition 5.1 (Generalized Latin Hypercube Sample).

Let $\pi_{i}$, $i=1,\ldots,p$ independent random permutations of
$0,\ldots,d-1$ and $\mathbf{V}^{i}=\left(V_{1}^{i},\ldots,V_{d}^{i}\right)$
$i=1,\ldots,p$ random vectors, independent from the $\pi_{i}$ and from each
other, identically distributed with probability measure $\mu$. A $p\times d$
matrix $\mathbf{U}$ is a Generalized Latin Hypercube Sample if has stochastic
representation

$\displaystyle
U_{l}^{i}=\frac{\pi_{i}\left(l\right)+1}{d}V^{i}_{l}+\frac{\pi_{i}\left(l\right)}{d}\left(1-V^{i}_{l}\right)\,i=1,\ldots,p,l=1,\ldots,d$
(71)

We remark that the most negative dependent constructions ILH(T),
LH$C_{d}\left(\left\\{1\right\\}\right)$, and LHAJ generate Generalized Latin
Hypercube Samples.

The following lemma introcuce the irrelevance of the distribution on the
$\mathbf{V}^{i}$, $i=1,\ldots,p$.

######  Lemma 5.1.

Consider a $a$-LH sample and a $b$-LH sample where $a$ and $b$ are two
different Radon measures. The following relationship holds for every function
$f$ locally integrable with respect to both measures :

$\displaystyle\mathbb{E}_{a-LH}\left[\left(\dfrac{1}{d}\sum^{d}_{l=1}f\left(\mathbf{U}_{l}\right)\right)^{r}\right]-\mathbb{E}_{b-LH}\left[\left(\dfrac{1}{d}\sum^{d}_{l=1}f\left(\mathbf{U}_{l}\right)\right)^{r}\right]=o\left(1\right)$

Given the previous lemma we are able to show that the asymptotic distribution
is the same of the ordinary Latin Hypercube. In particular Stein, (1987)
express the variance of the latin hypercube using the ANOVA decomposition of
the function $f$:

$\displaystyle f\left(\mathbf{u}\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{IID}\left[f\left(\mathbf{U}\right)\right]+\sum^{p}_{i=1}f_{i}\left(u_{i}\right)+r\left(\mathbf{u}\right)$
(72) $\displaystyle f_{i}\left(u_{i}\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{IID}\left[f\left(\mathbf{u}\right)-\mathbb{E}_{IID}\left[f\left(\mathbf{U}\right)\right]\left|U_{i}=u_{i}\right.\right]$

######  Theorem 5.1.

Let $\bar{X}=\dfrac{1}{d}\sum^{d}_{l=1}f\left(\mathbf{U}_{l}\right)$ with
$\mathbf{U}_{l}$ $l\in\mathcal{D}$ from a $\mu$-LH sample with $\mu$ Radon and
$f$ bounded and locally integrable with respect to $\mu$. Then
$\sqrt{d}\left(\bar{X}-E_{IID}\left(X\right)\right)$ converges in distribution
to
$\mathcal{N}\left(0,\int_{\left[0,1\right]^{p}}r\left(\mathbf{u}\right)^{2}d\mathbf{u}\right)$,
where $r\left(\mathbf{u}\right)$ is introduced in (72).

######  Remark 5.1.

The hypothesis of independence $\mathbf{U}^{j}$ independent from
$\mathbf{U}^{k}$ for all $j\neq k$ $j,k=1,\ldots,p$ in the definition of
Generalized Latin Hypercube Sample is not restrictive as it seems. In
practice, one can use the inverse Rosenblatt transform (see for example
Rüschendorf, (2013), Theorem 1.12) to obtain samples from a generic
distribution. Then, if the composition of $f$ with the inverse Rosenblatt
transform is bounded, we are still under the incidence of Theorem 5.1. A
multivariate version of the central limit theorem can be obtained as in
Corollary 1 of Owen, (1992). Concerning the boundedness assumption on $f$, it
is probably too restrictive because in Loh, (1996) a multivariate Berry-Essen
type Bound for the standardized multivariate version of $\bar{X}$, in the case
of Latin Hypercube, is obtained under the assumption that the multivariate
function $\mathbf{f}$ involved is Lebesgue measurable and
$\mathbb{E}\left\|\mathbf{f}\right\|^{3}<\infty$.

##  6 Numerical Illustrations

We illustrate the performance of the methods presented in this paper using
several simulation exercises involving standard Monte Carlo, Markov chain
Monte Carlo and Sequential Monte Carlo algorithms 101010Replication codes for
the examples in the paper can be found at
https://github.com/Frattalol/Livingontheedge.

One of the critical dimension used to rank our countermonotonic vectors is
sampling time, as shown in Figure 7. It is obtained by averaging over 1000
independent replications. In each experiment we sampled 5000 antithetic
vectors of dimension $2\leq d\leq 20$. Samplings schemes have been implemented
in Matlab on a Windows 10 laptop with an Intel i7-6500U CPU and 8 GB of RAM.

Figure 7: Sampling times (vertical axis) as function of the sampling dimension
$d$ (horizontal axis) for different methods (different symbols and colors).
All estimates are averages over 1000 experiments. In each experiments 5000
antithetic vectors of dimensions $d$ are sampled following a given method.

Permuted $C_{d}\left(\left\\{1\right\\}\right)$ outperforms the other
competitors. We report the times for the segmental version that randomly
permutes stochastic representation (10) and the RBS version that uses the
equivalence with random balanced sampling described in Proposition 4.1. The
latter choice is faster and it is also used as the base for the Latin
Hypercube iteration in $LH$ $C_{d}\left(\left\\{1\right\\}\right)$.

Consequently, all the experiment in the section use the exchangeable version
of LH $C_{d}\left(\left\\{1\right\\}\right)$. In fact, LH
$C_{d}\left(\left\\{1\right\\}\right)$ vector reach the lower value of
multivariate Kendall’s $\tau$ and Spearman’s $\rho$ (fig.6 ), but is also a
faster sampling scheme.

###  6.1 Monte Carlo Integration

The first evaluation of our methodology is for the Monte Carlo integration on
the unit hypercube, with integration points chosen according to the three
competing schemes. In standard Monte-Carlo (MC), the sampling points are iid,
and each sampling point corresponds to a random vector of dimension equal to
the number of variables in the function. In our antithetic method (LH
$C_{d}\left(\left\\{1\right\\}\right)$) the variables used to populate the
same coordinate for different sampling points are from an antithetic vector.
Across coordinates, these antithetic vectors are independent of one another.
For the Quasi Monte Carlo’s Sobol scheme (QMC Sobol’), a deterministic
sequence of points uniformly cover the hypercube of dimension equal to the
number of variables.In the evaluation of the integration problem complexity
the effective dimension plays an important role. We refer to Owen, (2003),
Wang and Fang, (2003) and Wang and Sloan, (2005) for the formal definition of
truncation $p_{t}$ and superposition $p_{s}$ dimensions. Owen, (2003) shows
that a low $p_{s}$ is necessary for QMC to surpass the computational
efficiency of MC when the sample sizes are at practical levels. Wang and Fang,
(2003) and Wang and Sloan, (2005) show that the integrands commonly used in
option pricing have $p_{t}\simeq p$ and $p_{s}\leq 2$ and explain, using those
results, QMC’s good performance in this domain. According to our Theorem 5.1
for integrands with $p_{s}=1$, i.e. for functions that are well approximated
by sums of one dimensional functions, LH
$C_{d}\left(\left\\{1\right\\}\right)$ should be efficient in reducing the
variance. Theorem 5.1 also guarantees that our method cannot perform worst
than MC, asymptotically in the number of points. To investigate the role of
effective dimension in the relative performance of the three competing
methods, we use the two parameter function introduced in Wang and Sloan,
(2005):

$\displaystyle
f\left(\mathbf{x}\right)=\prod^{p}_{i=1}\left(1+a\tau^{i}\left(x_{i}-1/2\right)\right).$
(73)

Varying the parameter $a$ has more effect on $p_{t}$ than on $p_{s}$ and
varying the parameter $\tau$ has the opposite effect. We consider an high
dimensional function ($p=100$) and different specifications of effective
dimensions, according to the parameters reported in Table 5.

Table 5: Wang and Fang, (2003) effective dimensions (truncation and superposition dimensions, $p_{t}$ and $p_{s}$, respectively) of the integrand function in Eq. (73) for dimension $p=100$ and different parameter settings (columns). $a$ | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 1 | 1 | 1 | 1 | 1 | 10 | 10 | 10 | 10 | 10
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\tau$ | 0.1 | 0.5 | 0.8 | 0.9 | 1 | 0.1 | 0.5 | 0.8 | 0.9 | 1 | 0.1 | 0.5 | 0.8 | 0.9 | 1
$p_{t}$ | 2 | 4 | 11 | 22 | 100 | 2 | 4 | 11 | 23 | 100 | 2 | 5 | 17 | 39 | 100
$p_{s}$ | 1 | 1 | 1 | 1 | 2 | 1 | 1 | 2 | 3 | 14 | 1 | 3 | 8 | 15 | 96

Effective dimensions are computed using the methods for multiplicative
functions introduced in Wang and Fang, (2003).

Figure 8 shows the mean square error (MSE) (vertical axis) and the computing
time (horizontal axis) of Monte Carlo (red), QMC (black) and circulant
variates LH-$C_{d}(\\{1\\})$ (blue) sampling, for the different effective
dimensions $p_{t}$ and $p_{s}$ (different plots) given in Table 5. Our
LH-$C_{d}(\\{1\\})$ method has the best performance when the superposition
dimension is equal to 1. The performance is decreasing in the truncation
dimension. For $p_{s}=2,3$, QMC performs better than the method proposed here
when the number of points used in the integration is high. QMC advantage
increases in the truncation dimension. At moderate $p_{s}$ QMC dominates. In
those cases, our proposal is slower than MC, but reaches the same MSE. In the
extreme case of $f$ almost full dimensional (right lower corner),
$C_{d}\left(\left\\{1\right\\}\right)$ is performing as bad as MC, but QMC is
doing orders of magnitudes worst. These numerical results are in line with the
result in Theorem 5.1 and indicate that our method should be used when the
superposition dimension is low and when there is no information about
effective dimensions of the integrand since in the worst case, it reproduce
the precision of standard MC estimates.

Figure 8: Mean square error (vertical axis) and computing time (horizontal
axis) for different effective dimensions $p_{t}$ and $p_{s}$ (different
plots). In each plot: Monte Carlo (red), QMC (black) and circulant variates
LH-$C_{d}(\\{1\\})$ (blue) sampling. For each line: different number of
samples from 10 to 1000 (circles). For each setting, all statistics are
averages over 10000 experiments.

###  6.2 Markov Chain Monte Carlo

####  6.2.1 Bayesian inference on Probit van Dyk and Meng, (2001)

The data used are taken from van Dyk and Meng, (2001) and represent the
clinical characteristics summarized by two covariates of 55 patients, of which
19 diagnosed with lupus. The disease indicator is modelled as independent
$Y_{i}\sim\mathcal{B}er\left(\Phi\left(x_{i}^{T}\beta\right)\right)$ where
$\Phi$ is the standard normal CDF and
$\beta=\left(\beta_{0},\beta_{1},\beta_{2}\right)^{T}$ is a the vector of
parameters. The objective is to sample from the posterior distribution
corresponding to the flat prior for $\beta$. We adopt the standard Gibbs
sampler with latent variables
$\Psi_{i}\sim\mathcal{N}\left(x_{i}^{T}\beta,1\right)$ of which we consider
only the sign. To obtain draws from the posterior we repeat the following
alternating two steps. First we sample from
$\left.\beta\right|\psi\sim\mathcal{N}\left(\tilde{\beta},\left(X^{T}X\right)^{-1}\right)$
with $\tilde{\beta}=\left(X^{T}X\right)^{-1}X^{T}\psi$ with $X$ the data
matrix whose $i$-th rows is $x_{i}$. Then from
$\psi_{i}\left|\beta,Y_{i}\right.\sim\mathcal{TN}\left(x_{i}^{T}\beta,1,Y_{i}\right)$
where $\mathcal{TN}\left(\mu,\sigma^{2},Y\right)$ is a the normal distribution
with mean $\mu$ and variance $\sigma^{2}$, truncated to be positive if $Y>0$
or negative otherwise. Further details of the algorithms can be found in
(Craiu and Meng,, 2005).

Figure 9: Asymptotic variance ratio between antithetic and iid sampling
(vertical axis) for different dimension $d$ (horizontal axis) of antithetic
variates, and parameters of the Bayesian probit model (different plots). In
each plot: the average ratio (dots) and its range (vertical segments). Note:
all estimates are based on 1000 independently-replicated experiments. In each
experiment 5000 Gibbs iterations with antithetic vectors of dimensions $d$ are
used.

####  6.2.2 Metropolis within Gibbs Hierarchical Poisson Gelfand and Smith,
(1990)

This second example concerns the counts of failures
$s=\left(s_{1},\ldots,s_{n}\right)$ for $n=10$ pumps in a nuclear power plant.
The time $t=\left(t_{1},\ldots,t_{n}\right)$ of operation are known. The model
assumes $s_{k}\sim\mathcal{P}oi\left(\lambda_{k}t_{k}\right)$ and
$\lambda_{k}\sim\mathcal{G}a\left(\alpha,\beta\right)$. The objective of the
inference are draws from the posterior distribution of the parameters $\alpha$
to which we assign an exponential prior with mean 1 and $\beta$ with a Gamma
prior distribution $\mathcal{G}a\left(0.1,1\right)$. Using the conjugate
priors it is easy to obtain a Gamma distribution for
$\lambda_{1}\left.|\lambda_{2},\ldots\lambda_{n},\beta\right.$ and those
variables can be easily sampled using a Gibbs step. Sampling of $\alpha$ is
slightly more difficult. In fact we have :

$\displaystyle\mathbb{P}\left(\alpha\left.|\lambda_{1},\ldots\lambda_{n},\beta\right.\right)\propto\exp\left[\alpha\left(n\log\beta+\sum^{n}_{k=1}\log\lambda_{k}-1\right)-n\log\Gamma\left(\alpha\right)\right]$

We then sample $\alpha$ with a random walk Metropolis Hastings (MH) step with
a deterministic scan. We show results for the case when we antithetically
couple also the uniform draws for acceptance rejection choice and when we are
not doing it. The former case is the one for which Frigessi et al., (2000)
report the worst performance of the usual two variates antithetic coupling of
chains, in agreement with our results.

Figure 10: Asymptotic variance ratio between antithetic and iid sampling
(vertical axis) for different dimension $d$ (horizontal axis) of antithetic
variates, and parameters of the Bayesian Hierarchical Poisson model (different
plots). In each plot: the average ratio for the MH with (blue dots) and
without (red dots) antithetic acceptance rule and their ranges (vertical
segments). Note: all estimates are based on 1000 independently-replicated
experiments. In each experiment 5000 Gibbs iterations with antithetic vectors
of dimensions $d$ are used.

####  6.2.3 Pseudo Marginal Metropolis Hastings Stochastic Volatility Gerber
and Chopin, (2015)

The last application, target a state of the art methodology, the Pseudo
Marginal MH (PCMH) proposed by Andrieu and Roberts, (2009) and which is able
to estimate models with intractable likelihoods that are approximated using a
particle filter. In particular, following Gerber and Chopin, (2015), we
consider the bivariate stochastic volatility model introduced in Chan et al.,
(2006):

$\displaystyle\mathbf{y}_{t}$ $\displaystyle=$ $\displaystyle
S_{t}^{1/2}\mathbf{\epsilon}_{t}$ $\displaystyle\mathbf{x}_{t}$
$\displaystyle=$
$\displaystyle\mathbf{\mu}+\Phi\left(\mathbf{x}_{t-1}-\mathbf{\mu}\right)+\Psi^{1/2}\nu_{t}$
$\displaystyle S_{t}$ $\displaystyle=$
$\displaystyle\hbox{diag}\left(\exp\left(x_{1t},x_{2t}\right)\right)$
$\displaystyle\left(\mathbf{\epsilon}_{t},\mathbf{\nu}_{t}\right)$
$\displaystyle\sim$ $\displaystyle\mathcal{N}\left(\mathbf{0}_{4},C\right)$

with $\Phi$ and $\Psi$ diagonal matrices and $C$ a correlation matrix.
Following those authors we take the prior:

$\displaystyle\phi_{ii}$ $\displaystyle\sim$ $\displaystyle\mathcal{U}[0,1]$
(74) $\displaystyle\dfrac{1}{\psi_{ii}}$ $\displaystyle\sim$
$\displaystyle\mathcal{G}a\left(10\exp\left(-10\right)/,10\exp\left(-3\right)\right)$
(75)

where $\phi_{ii}$ and $\psi_{ii}$ denote respectively the diagonal elements of
$\Phi$ and $\Psi$, and a flat prior for $\mu$. In addition, we assume that C
is uniformly distributed on the space of correlation matrices. To sample from
the posterior distribution of the parameters, we use a Gaussian random-walk MH
algorithm with covariance matrix calibrated by Gerber and Chopin, (2015) so
that the acceptance probability of the algorithm becomes, as $N$ tends to
infinity, close to 25%. We consider the mean-corrected daily returns on the
Nasdaq and Standard and Poor’s 500 indices for the period ranging from January
3rd, 2012, to October 21st, 2013, so that the data set contains 452
observations. Figure 11

Figure 11: PCMH using Sequential Monte Carlo (SMC), Sequential Quasi-Monte
Carlo (SQMC) and Sequential Antithetic Monte Carlo (SAMC) (different colors).
Acceptance rate of the Metropolis step (vertical axis) versus number of
particles (horizontal axis, Panel A) and computing time (horizontal axis,
Panel B). Maximum and minimum effective sample size (vertical axis) versus
number of particles (horizontal axis, Panel C) and computing time (horizontal
axis, Panel D).

PCMH algorithms using sequential quasi-Monte Carlo and antithetic Monte Carlo
are equivalent in acceptance rate and effective sample size (Panels A and C)
when a low number of particles (up to 20) is used. Nevertheless, antithetic
Monte Carlo achieves larger acceptance rates (AR) and effective sample size
(ESS) with a lower computing time (Panel B and D). When a larger number of
particles is used (above 20), the performances are equivalent in terms of ESS,
whereas SQMC is better in terms of AR.

##  7 Discussion

The development of antithetic constructions has generated a rich class of
methods to accompany the evolution of Monte Carlo sampling algorithms. We
enrich this class with a new antithetic method, the circulant variates (CCV),
that satisfies the countermonotonicity, exchangeability, and marginal
uniformity conditions. In particular, the marginal uniformity condition is
linked to the Kullback-Leibler optimality.

The principle behind the proposal, relying on sampling on segments, leads to a
unification of several classical antithetic constructions: rotation sampling,
Latin hypercube, permuted displacement, and random balanced sampling.

Within this common framework, we evaluate theoretically the resulting
distributions of the antithetic vectors and their concordance measures. The
latter allows us to rank the methods within the class of sampling on segments.
We also demonstrate a central limit theorem in the case of asymptotically
increasing vector size.

Leveraging on iterative Latin Hypercube (LH) properties, we combine the two
methods by using the CCV construction to initialize the LH construction and
reduce the number of iterations. This reduces the simulation cost and improves
the concordance lower bound. The numerical experiments include MCMC,
Sequential MC, Quasi-MC, and classical MC integration. The methods proposed
outperform standard implementations and are competitive with Quasi-Monte Carlo
methods in low effective-dimension scenarios. In addition, Theorem 5.1 is
confirmed by our numerical experiments, which show variance reduction at least
as good as standard MC.

Future work includes possible extensions of the theory for KL optimality
beyond the marginal univariate case. An investigation of the relationship
between the line segment representation and orthogonal array-based latin
hypercubes Tang, (1993) could lead to an improvement of the performance for
superposition dimensions bigger than 1. We notice that different methods that
satisfy countermonotonicity, exchangeability, and marginal uniformity yield
different variance reductions in practice. A more general and ambitious aim is
to propose a mathematical framework directed at identifying the additional
features that produce these differences. Our sampling method has been
successfully applied within the Bayesian estimation framework of the European
Commission’s multi-country model (Albonico et al.,, 2019). We expect that the
proposed simulation technique will find direct application in other fields of
computational mathematics and statistics.

## Acknowledgements

The views expressed by Lorenzo Frattarolo are the author’s alone and do not
necessarily correspond to those of the European Commission. We thank for the
useful comments the participants at: the 1st Italian-French Statistics
Workshop, University Ca’ Foscari of Venice, October 2017, Venice, Italy; the
32nd European Statisticians Meeting, Palermo, Italy, 2019; the Workshop ”Monte
Carlo methods, and approximate dynamic programming”, ESSEC, October 2019,
Cergy, France; the seminar series at the University of Surrey. Radu V. Craiu’s
research was supported by NSERC of Canada grant RGPIN-249547. Roberto Casarin
acknowledges support from the Venice Centre for Risk Analytics (VERA) at the
University Ca’ Foscari of Venice.

## References

  * Ahn and Fuchs, (2020) Ahn, J. Y. and Fuchs, S. (2020).  On minimal copulas under the concordance order.  Journal of Optimization Theory and Applications, 184(3):762–780.
  * Albonico et al., (2019) Albonico, A., Calés, L., Cardani, R., Croitorov, O., Ferroni, F., Giovannini, M., Hohberger, S., Pataracchia, B., Pericoli, F. M., Raciborski, R., Ratto, M., Roeger, W., and Vogel, L. (2019).  Comparing post-crisis dynamics across euro area countries with the global multi-country model.  Economic Modelling, 81:242–273.
  * Aliprantis and Border, (2007) Aliprantis, C. and Border, K. (2007).  Infinite Dimensional Analysis: A Hitchhiker’s Guide.  Springer.
  * Andréasson, (1972) Andréasson, I. J. (1972).  Combinations of antithetic methods in simulation.  Kungliga Tekniska Högskolan. Department of Information Processing. Computer ….
  * Andréasson and Dahlquist, (1972) Andréasson, I. J. and Dahlquist, G. (1972).  Groups of antithetic transformations in simulation.  Kungliga Tekniska Högskolan. Department of Information Processing. Computer ….
  * Andrieu and Roberts, (2009) Andrieu, C. and Roberts, G. O. (2009).  The pseudo-marginal approach for efficient monte carlo computations.  Ann. Statist., 37(2):697–725.
  * Arvidsen and Johnsson, (1982) Arvidsen, N. I. and Johnsson, T. (1982).  Variance reduction through negative correlation, a simulation study.  Journal of Statistical Computation and Simulation, 15(2-3):119–127.
  * Barlow and Proschan, (1975) Barlow, R. and Proschan, F. (1975).  Statistical Theory of Reliability and Life Testing: Probability Models.  International series in decision processes. Holt, Rinehart and Winston.
  * Besicovitch, (1945) Besicovitch, A. S. (1945).  A general form of the covering principle and relative differentiation of additive functions.  Mathematical Proceedings of the Cambridge Philosophical Society, 41(2):103–110.
  * Birkhoff, (1946) Birkhoff, G. (1946).  Three observations on linear algebra.  Univ. Nac. Tacuman, Rev. Ser. A, 5:147–151.
  * Brown, (1966) Brown, J. (1966).  Approximation theorems for markov operators.  Pacific Journal of Mathematics, 16(1):13–23.
  * Bubenik and Holbrook, (2007) Bubenik, O. and Holbrook, J. (2007).  Densities for random balanced sampling.  Journal of Multivariate Analysis, 98(2):350 – 369.
  * Chan et al., (2006) Chan, D., Kohn, R., and Kirby, C. (2006).  Multivariate stochastic volatility models with correlated errors.  Econometric Reviews, 25(2-3):245–274.
  * Craiu and Meng, (2005) Craiu, R. V. and Meng, X.-L. (2005).  Multiprocess parallel antithetic coupling for backward and forward Markov chain Monte Carlo.  The Annals of Statistics, 33(2):661–697.
  * Craiu and Meng, (2006) Craiu, R. V. and Meng, X.-L. (2006).  Meeting hausdorff in monte carlo: A surprising tour with antihype fractals.  Statistica Sinica, 16(1):77–91.
  * Dhaene and Goovaerts, (1996) Dhaene, J. and Goovaerts, M. J. (1996).  Dependency of risks and stop-loss order 1.  ASTIN Bulletin: The Journal of the IAA, 26(2):201–212.
  * Durante et al., (2014) Durante, F., Fernández Sánchez, J., and Trutschnig, W. (2014).  Multivariate copulas with hairpin support.  Journal of Multivariate Analysis, 130:323–334.
  * Durante and Sanchez, (2012) Durante, F. and Sanchez, J. F. (2012).  On the approximation of copulas via shuffles of min.  Statistics & Probability Letters, 82(10):1761 – 1767.
  * Evans and Gariepy, (1991) Evans, L. and Gariepy, R. (1991).  Measure Theory and Fine Properties of Functions.  Studies in Advanced Mathematics. Taylor & Francis.
  * Fishman and Huang, (1983) Fishman, G. S. and Huang, B. D. (1983).  Antithetic variates revisited.  Commun. ACM, 26(11):964–971.
  * Fréchet, (1935) Fréchet, M. (1935).  Généralisation du théoreme des probabilités totales.  Fundamenta Mathematicae, 25(1):379–387.
  * Fréchet, (1951) Fréchet, M. (1951).  Sur les tableaux de corrélation dont les marges sont données.  Ann. Univ. Lyon, 3ê serie, Sciences, Sect. A, 14:53–77.
  * Frigessi et al., (2000) Frigessi, A., Gåsemyr, J., and Rue, H. (2000).  Antithetic coupling of two Gibbs sampler chains.  The Annals of Statistics, 28(4):1128–1149.
  * Fuchs et al., (2018) Fuchs, S., McCord, Y., and Schmidt, K. D. (2018).  Characterizations of copulas attaining the bounds of multivariate kendall’s tau.  Journal of Optimization Theory and Applications, 178(2):424–438.
  * Gaffke and Rüschendorf, (1981) Gaffke, N. and Rüschendorf, L. (1981).  On a class of extremal problems in statistics.  Mathematische Operationsforschung und Statistik. Series Optimization, 12(1):123–135.
  * Gelfand and Smith, (1990) Gelfand, A. E. and Smith, A. F. M. (1990).  Sampling-based approaches to calculating marginal densities.  Journal of the American Statistical Association, 85(410):398–409.
  * Gerber and Chopin, (2015) Gerber, M. and Chopin, N. (2015).  Sequential quasi monte carlo.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(3):509–579.
  * Gerow and Holbrook, (1996) Gerow, K. and Holbrook, J. (1996).  Statistical sampling and fractal distributions.  Math. Intelligencer, 18:12–22.
  * Hammersley and Morton, (1956) Hammersley, D. C. and Morton, K. V. (1956).  A new Monte Carlo technique: antithetic variates.  Mathematical Proceedings of the Cambridge Philosophical Society, 52:449–475.
  * Hammersley and Mauldon, (1956) Hammersley, J. M. and Mauldon, J. G. (1956).  General principles of antithetic variates.  Mathematical Proceedings of the Cambridge Philosophical Society, 52(3):476–481.
  * Handscomb, (1958) Handscomb, D. C. (1958).  Proof of the antithetic variates theorem for n ¿ 2.  Mathematical Proceedings of the Cambridge Philosophical Society, 54(2):300–301.
  * Hardy et al., (1934) Hardy, G., Collection, K. M. R., Littlewood, J., Pólya, G., Pólya, G., and Littlewood, D. (1934).  Inequalities.  Cambridge University Press.
  * Hoeffding, (1940) Hoeffding, W. (1940).  Masstabinvariante korrelationtheorie.  Schriften des Mathematischen Instituts und des Instituts für Angewandte Mathematik der Universitat Berlin, 5:179–233.
  * Isbell, (1955) Isbell, J. R. (1955).  Birkhoff’s problem 111.  Proceedings of the American Mathematical Society, 6(2):217–218.
  * Joe, (1990) Joe, H. (1990).  Multivariate concordance.  Journal of Multivariate Analysis, 35(1):12 – 30.
  * Knott and Smith, (2006) Knott, M. and Smith, C. (2006).  Choosing joint distributions so that the variance of the sum is small.  Journal of Multivariate Analysis, 97(8):1757 – 1765.
  * Lee and Ahn, (2014) Lee, W. and Ahn, J. Y. (2014).  On the multidimensional extension of countermonotonicity and its applications.  Insurance: Mathematics and Economics, 56(Supplement C):68 – 79\.
  * Lee et al., (2017) Lee, W., Cheung, K. C., and Ahn, J. Y. (2017).  Multivariate countermonotonicity and the minimal copulas.  Journal of Computational and Applied Mathematics, 317:589–602.
  * Lemieux, (2009) Lemieux, C. (2009).  Monte Carlo and Quasi-Monte Carlo Sampling.  Springer Series in Statistics. Springer New York.
  * Lévy, (1962) Lévy, P. (1962).  Extensions d’un théorème de d. dugué et m. girault.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 1(2):159–173.
  * Lindenstrauss, (1965) Lindenstrauss, J. (1965).  A remark on extreme doubly stochastic measures.  The American Mathematical Monthly, 72(4):379–382.
  * Loh, (1996) Loh, W. L. (1996).  On Latin hypercube sampling.  The Annals of Statistics, 24:2058–2080.
  * Lu and Yi, (2004) Lu, T.-Y. and Yi, Z. (2004).  Generalized correlation order and stop-loss order.  Insurance: Mathematics and Economics, 35(1):69–76.
  * Matoušek, (1998) Matoušek, J. (1998).  On the $l^{2}$-discrepancy for anchored boxes.  Journal of Complexity, 14(4):527–556.
  * McKay et al., (1979) McKay, M. D., Beckman, R. J., and Conover, W. J. (1979).  A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.  Technometrics, 21:239–245.
  * Mikusiński and Taylor, (2010) Mikusiński, P. and Taylor, M. D. (2010).  Some approximations of n-copulas.  Metrika, 72(3):385–414.
  * Moameni, (2016) Moameni, A. (2016).  Supports of extremal doubly stochastic measures.  Canadian Mathematical Bulletin, 59(2):381–391.
  * Nelsen and Úbeda-Flores, (2012) Nelsen, R. B. and Úbeda-Flores, M. (2012).  Directional dependence in multivariate distributions.  Annals of the Institute of Statistical Mathematics, 64(3):677–685.
  * Owen, (1992) Owen, A. B. (1992).  A central limit theorem for latin hypercube sampling.  Journal of the Royal Statistical Society: Series B (Methodological), 54(2):541–551.
  * Owen, (2003) Owen, A. B. (2003).  The dimension distribution, and quadrature test functions.  Statistica Sinica, 13:1–17.
  * Puccetti and Wang, (2015) Puccetti, G. and Wang, R. (2015).  Extremal Dependence Concepts.  Statistical Science, 30(4):485 – 517.
  * Roach and Wright, (1977) Roach, W. and Wright, R. (1977).  Optimal antithetic sampling plans.  Journal of Statistical Computation and Simulation, 5(2):99–114.
  * Rockafellar, (1970) Rockafellar, R. (1970).  Convex Analysis.  Princeton Landmarks in Mathematics and Physics. Princeton University Press.
  * Rubinstein and Samorodnitsky, (1987) Rubinstein, Y. R. and Samorodnitsky, G. (1987).  A modified version of handscomb’s antithetic variates theorem.  SIAM Journal on Scientific and Statistical Computing, 8(1):82–98.
  * Rüschendorf, (2013) Rüschendorf, L. (2013).  Mathematical Risk Analysis: Dependence, Risk Bounds, Optimal Allocations and Portfolios.  Springer Series in Operations Research and Financial Engineering. Springer Berlin Heidelberg.
  * Rüschendorf and Uckelmann, (2002) Rüschendorf, L. and Uckelmann, L. (2002).  Variance Minimization and Random Variables with Constant Sum, pages 211–222.  Springer Netherlands, Dordrecht.
  * Shepp, (1962) Shepp, L. (1962).  Symmetric random walk.  Transactions of the American Mathematical Society, 104(1):144–153.
  * Stein, (1987) Stein, M. (1987).  Large sample properties of simulations using Latin hypercube sampling.  Technometrics, 29:143–151.
  * Tang, (1993) Tang, B. (1993).  Orthogonal array-based latin hypercubes.  Journal of the American Statistical Association, 88(424):1392–1397.
  * Tukey, (1957) Tukey, J. W. (1957).  Antithesis or regression?  Mathematical Proceedings of the Cambridge Philosophical Society, 53(4):923–924.
  * van Dyk and Meng, (2001) van Dyk, D. A. and Meng, X.-L. (2001).  The art of data augmentation.  Journal of Computational and Graphical Statistics, 10(1):1–50.
  * Vitale, (1990) Vitale, R. A. (1990).  On stochastic dependence and a class of degenerate distributions.  Topics in Statistical Dependence, 16:459.
  * Von Neumann, (1953) Von Neumann, J. (1953).  A certain zero-sum two-person game equivalent to the optimal assignment problem.  Contributions to the Theory of Games, 2:5–12.
  * Wang and Wang, (2011) Wang, B. and Wang, R. (2011).  The complete mixability and convex minimization problems with monotone marginal densities.  Journal of Multivariate Analysis, 102(10):1344 – 1360.
  * Wang and Fang, (2003) Wang, X. and Fang, K.-T. (2003).  The effective dimension and quasi-monte carlo integration.  Journal of Complexity, 19(2):101–124.
  * Wang and Sloan, (2005) Wang, X. and Sloan, I. H. (2005).  Why are high-dimensional finance problems often of low effective dimension?  SIAM Journal on Scientific Computing, 27(1):159–183.
  * Whitt, (1976) Whitt, W. (1976).  Bivariate Distributions with Given Marginals.  The Annals of Statistics, 4(6):1280 – 1289.
  * Wilson, (1979) Wilson, J. R. (1979).  Proof of the antithetic-variates theorem for unbounded functions.  Mathematical Proceedings of the Cambridge Philosophical Society, 86(3):477–479.
  * Wilson, (1983) Wilson, J. R. (1983).  Antithetic sampling with multivariate inputs.  American Journal of Mathematical and Management Sciences, 3(2):121–144.

##  Appendix A Proof of Theorems and Lemmas

###  A.1 Proof of Lemma 2.1

###### Proof.

By (10), conditionally on being on the $k$-th edge
$e_{k}=\left(i\left(k\right),j\left(k\right)\right)$, $U_{l}$ is a standard
uniform on $\left[\alpha_{lk},\beta_{lk}\right]$ if the interval has length
bigger than 0 and a Dirac mass centered in $x_{lj\left(k\right)}$ otherwise
i.e. has a density:

$\displaystyle
f\left(\right.u\left|K=k\right)=\left\\{\begin{array}[]{ccc}\dfrac{1}{\left(\beta_{lk}-\alpha_{lk}\right)}\mathbb{I}_{\left\\{\left[\alpha_{lk},\beta_{lk}\right]\right\\}}\left(u_{l}\right)&\hbox{if}&x_{li\left(k\right)}\neq
x_{lj\left(k\right)}\\\ &&\\\
\mathbb{I}_{\left\\{x_{lj\left(k\right)}\right\\}}\left(u_{l}\right)&\hbox{if}&x_{li\left(k\right)}=x_{lj\left(k\right)}\end{array}\right.$

We get rid of the case $x_{li\left(k\right)}=x_{lj\left(k\right)}$ by
excluding segments contained in hyperplanes parallel to the hyperfaces of the
hypercube, obtaining uniformity on $\left[\alpha_{lk},\beta_{lk}\right]$.

∎

###  A.2 Proof of Theorem 2.1

Before stating the validity of the result in Theorem 2.1 we show the following
lemma.

######  Lemma A.1.

Let $U_{l}$ be the $l$-th component of random vector
$\mathbf{U}=(U_{1},\ldots,U_{d})$ in the stochastic representation (10) and
let

$\mathcal{K}_{u_{l}}=\left\\{k\in\left\\{1,\ldots,\left|\mathcal{E}\right|\right\\}:u_{l}\in\left[a_{l,m_{\alpha}\left(k\right)},a_{l,m_{\beta}\left(k\right)}\right]\right\\}$

for each $u_{l}\in\left[0,1\right]$. Under Assumption 1 and 2 if

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{u_{l}}}\dfrac{1}{a_{l,m_{\beta}\left(k\right)}-a_{l,m_{\alpha}\left(k\right)}},$

then $U_{l}$ has standard uniform marginal density. If $u_{l}$ belongs to the
element $A_{l,m}$ of the partition defined in Equations we get
$\mathcal{K}_{u_{l}}=\mathcal{K}_{l,m}$.

###### Proof.

When $u_{l}\in\left[0,1\right]$ the marginal standard uniform is obtained by
equating the piece-wise uniform PDF of $U_{l}$ to $1$:

$\displaystyle f\left(u_{l}\right)$ $\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum^{\left|\mathcal{E}\right|}_{k=1}f\left(\left.u_{l}\right|K=k\right)$
$\displaystyle=$
$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum^{\left|\mathcal{E}\right|}_{k=1}\dfrac{1}{\left(\beta_{lk}-\alpha_{lk}\right)}\mathbb{I}_{\left\\{\left[\alpha_{lk},\beta_{lk}\right]\right\\}}\left(u_{l}\right)=1,$

which is satisfied by choosing $\mathcal{S}=(\mathcal{G},\mathbf{X})$ and
$\alpha_{lk}$ and $\beta_{lk}$ as in the definition of admissibility of given
in Remark 2.1. ∎

In the previous lemma the elements in the sum depend on $u_{l}$ and the
equation (A.1) has to be checked for infinitely many cases. The proof of
theorem 2.1 apply lemma A.1 and the properties of the finite partition (20) of
the interval $\left[a_{l,1}=0,a_{l,n}=1\right]$, to show that condition (A.1)
need to be verified only for a finite number of points.

###### Proof.

Consider the interval $A_{l,m}=\left[a_{l,m-1},a_{l,m}\right)$, There are two
possibilities
$A_{l,m}\cap\left[a_{l,m_{\alpha}\left(k\right)},a_{l,m_{\beta}\left(k\right)}\right]=\emptyset$,
when $a_{l,m_{\alpha}\left(k\right)}>=a_{l,m}$ or
$a_{l,m_{\beta}\left(k\right)}<a_{l,m-1}$, and
$A_{l,m}\subseteq\left[a_{l,m_{\alpha}\left(k\right)},a_{l,m_{\beta}\left(k\right)}\right]$
when we have jointly $a_{l,m_{\alpha}\left(k\right)}<=a_{l,m-1}$ and
$a_{l,m_{\beta}\left(k\right)}>=a_{l,m}$.

If $u_{l}\in A_{l,m}$ and $m_{\alpha}\left(k\right)+1\leq m\leq
m_{\beta}\left(k\right)$ then
$u_{l}\in\left[a_{l,m_{\alpha}\left(k\right)},a_{l,m_{\beta}\left(k\right)}\right]$.
Thus $u_{l}\in A_{l,m}$ implies $\mathcal{K}_{u_{l}}=\mathcal{K}_{l,m}$ and
$\mathcal{K}_{u_{l}}$ does not depend on $u_{l}$.

$a_{l,m_{\beta}\left(k\right)}$$a_{l,m}$$a_{l,m_{\alpha}\left(k\right)}$$a_{l,m-1}$

Let us show that if (A.1) is valid for the all the $A_{l,m}$, then it is valid
for every other partition of the unit interval. We start by verifying that if
condition (A.1) is verified by each disjoint interval in the partition it is
also verified by intervals from coarser partitions. Consider partitions where
we join two different intervals $B=A_{l,m}\cup A_{l,m+1}$. We should verify
(A.1) for all $u_{l}\in B$ but since the intervals are disjoint it suffice to
verify separately for all $u_{l}\in A_{l,m}$ and for all $u_{l}\in A_{l,m+1}$
and those conditions are satisfied by assumption. Let us show, now, that if
condition (A.1) is verified by each disjoint interval in the partition it is
also verified by intervals from finer partitions. Lets divide $A_{l,m}$ in two
disjoint subsets by choosing an arbitrary point $a_{*}\in A_{l,m}$.

$A^{1}_{l,m}=\left[a_{l,m},a_{*}\right),\quad
A^{2}_{l,m}=\left[a_{*},a_{l,m+1}\right).$

Let us call $\mathcal{K}_{i}=\mathcal{K}_{u_{l}}$ when $u_{l}\in A^{i}_{l,m}$,
$i=1,2$.

Since no segment in $\mathcal{S}$ has a vertex with $l$-th coordinate $a_{*}$,
when $k\in\mathcal{K}_{l,m}$, we have $a_{l,m_{\alpha}\left(k\right)}\leq
a_{l,m}<a_{*}$ and $a_{l,m_{\beta}\left(k\right)}\geq a_{l,m+1}>a_{*}$. Then,
we have jointly $a_{l,m_{\alpha}\left(k\right)}<a_{*}$ and
$a_{l,m_{\beta}\left(k\right)}\leq a_{l,m+1}$ and
$a_{l,m_{\alpha}\left(k\right)}\leq a_{l,m}$ and
$a_{l,m_{\beta}\left(k\right)}>a_{*}$ . This means:

$\displaystyle\mathcal{K}_{1}=\mathcal{K}_{2}=\mathcal{K}_{l,m},$

verifying (A.1) for $A^{1}_{l,m}$ and $A^{2}_{l,m}$ i.e. for the finer
partition. Every other partition can be obtained by iteratively splitting and
joining intervals.

Finally we show that condition (A.1) is satisfied on the last interval
$A_{l,n_{l}}$ if it is verified on every other interval and that the
probability mass is 1. Since $U_{l}$ is a mixture of uniform random variables
on the range $\left[0,1\right]$ the PDF integrates to one and exploiting
$\left[0,1\right]=(\left[0,1\right]/A_{l,n_{l}})\cup A_{l,n_{l}}$ one obtains:

$\displaystyle
1=\int^{1}_{0}f\left(u_{l}\right)du_{l}=\int^{a_{l,n_{l}-1}}_{0}f\left(u_{l}\right)du_{l}+\int^{1}_{a_{l,n_{l}-1}}f\left(u_{l}\right)du_{l}$
(77)

If (A.1) is verified for $\left\\{A_{l,m}\right\\}^{n_{l}-1}_{m=2}$, then:

$\displaystyle\int^{a_{l,n_{l}-1}}_{0}f\left(u_{l}\right)du_{l}=a_{l,n_{l}-1},$
$\displaystyle
1=a_{l,n_{l}-1}+\int^{1}_{a_{l,n_{l}-1}}\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{l,n_{l}}}\dfrac{1}{1-a_{l,m_{\alpha}\left(k\right)}}du_{l}$
$\displaystyle
1=a_{l,n_{l}-1}+\left(1-a_{l,n_{l}-1}\right)\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{l,n_{l}}}\dfrac{1}{1-a_{l,m_{\alpha}\left(k\right)}},$

or equivalently

$\displaystyle\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{l,n_{l}}}\dfrac{1}{1-a_{l,m_{\alpha}\left(k\right)}}=1.$

∎

###  A.3 Proof of Theorem 2.2

In the proof we will use the following lemma.

######  Lemma A.2.

The system in (21) is equivalent to the system :

$\displaystyle G_{l,m}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle\displaystyle-\dfrac{1}{\left|\mathcal{E}\right|}\sum^{n_{l}}_{m^{\prime}=1}\dfrac{n^{l}_{\left(m,m^{\prime}\right)}}{\left|a_{l,m^{\prime}}-a_{l,m}\right|}=0$
$\displaystyle m$ $\displaystyle=$ $\displaystyle 2,\ldots,n_{l}-1$ (78)

###### Proof.

Recall the definition in Eq. (21)

$F_{l,m}\left(\mathbf{a}_{l}\right)=\dfrac{1}{\left|\mathcal{E}\right|}\displaystyle\sum_{k\in\mathcal{K}_{l,m}}\dfrac{1}{a_{l,m_{\beta}\left(k\right)}-a_{l,m_{\alpha}\left(k\right)}}-1,$

then $G_{l,m}(\mathbf{a}_{l})$ writes as:

$\displaystyle G_{l,m}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle
F_{l,m}\left(\mathbf{a}_{l}\right)-F_{l,m+1}\left(\mathbf{a}_{l}\right)$
$\displaystyle m$ $\displaystyle=$ $\displaystyle 2,\ldots,n_{l}-1$ (79)
$\displaystyle G_{l,n_{l}}\left(\mathbf{a}_{l}\right)$ $\displaystyle=$
$\displaystyle F_{l,n_{l}}\left(\mathbf{a}_{l}\right).$

Then: