domain: stat math

text:
Olivier P Faugeras

The Stick-Breaking and Ordering Representation of Compositional Data: Copulas and Regression models

. We propose an alternative view of CoDa as stick breaking processes. The first stick-breaking approach gives rise to a view of CoDa as ordered statistics, from which we can derive "stick-ordered" distributions. The second approach is based on a rescaled stick-breaking transformation, and give rises to a geometric view of CoDa as a free unit cube. The latter allows to introduce copula and regression models, which are useful for studying the internal or external dependence of CoDa. We establish connections with other topics of statistics like i) spacings and order statistics, ii) Bayesian nonparametrics and Dirichlet distributions, iii) neutrality, iv) mixability.

Introduction and summary

Compositional Data (CoDa) analysis deals with statistical analysis of d-variate data which are quantitative descriptions of the parts of some whole, conveying only relative information. Composition of soil in geology, elements in a mixture in chemistry, sources of calories in nutrition, or vote shares in an election are examples of CoDa.

There are several competing ways to describe CoDa. The traditional approach (J. [START_REF] Aitchison | The statistical analysis of compositional data[END_REF], [START_REF] Boogaart | Analyzing compositional data with R. Use R! Springer[END_REF]) considers that one of the key characteristics of CoDa is that the sum of the proportions must always be equal to a constant (w.l.o.g. 1). This means that the different components of a CoDa point are often expressed as percentages or fractions, rather than absolute values. Hence, Aitchison's approach normalizes a raw composition vector by its sum (an operation called closure in the CoDa literature): let y = (y 1 , . . . , y d ) ∈ R d ≥0 be a vector of non-negative absolute values of a 1 composition, its closure is denoted as

C(y) := y ||y|| 1 = y d i=1 y i .
This leads to the consideration of normalised (i.e. after rescaling to unit sum) CoDa element as a vector p = (p 1 , . . . , p d ), taking its values in the d -1 dimensional unit simplex1

∆ d-1 := {p ∈ R d : p i ≥ 0, d i=1 p i = 1}. ( 1
)
CoDa points are then traditionally studied through a variety of log-ratio transforms, an approach pioneered by Aitchison John [START_REF] Aitchison | The statistical analysis of compositional data[END_REF], J. [START_REF] Aitchison | The statistical analysis of compositional data[END_REF]. This gives rise to a special geometry, called Aitchison geometry, which turns the positive simplex into an Euclidean vector space. For recent accounts on the latter, see e.g., [START_REF] Greenacre | Compositional data analysis in practice[END_REF] or [START_REF] Pawlowsky-Glahn | Modeling and analysis of compositional data[END_REF]. For an intrinsic approach to CoDa analysis based on projective geometry, see [START_REF] Faugeras | An invitation to intrinsic compositional data analysis using projective geometry and Hilbert's metric[END_REF].

Inspired by Bayesian nonparametrics, we propose an alternative view of CoDa points, as a stick-breaking process. This give rises to two interrelated geometric views on CoDa points. The first is that of an ordered set on the unit interval. This ordered view allows to define distributions on CoDa points via order statistics and spacings of a latent vector on the unit interval.

The second view is based on a transformation which gives the relative positions of the breaks, yielding a parametrization of the CoDa point as an unconstrained unit cube. These relative positions have an interpretation as conditional probabilities and are related to the concept of neutrality, a natural intra-independence notion for compositions. This second view also allows to define distributions on the simplex, in particular copula distributions, for the study of the intra-dependence of CoDa. In addition, it is useful to study internal (resp. external) regression models, i.e. when one wants to explain/predict a (set of) components by other components acting as predictors (resp. by external covariates).

The outline of the paper is as follows: in Section 2, we introduce our first stick breaking transformation and define stick-breaking distributions for CoDa, based on spacings and order statistics. Several examples are given and numerically illustrated. Section 3 elaborates on the first construction by considering a rescaled version of the stick-breaking process. It turns the constrained CoDa point of the simplex into a free vector of the unit cube [0, 1] d-1 , which can, for positive CoDa, even be transformed to a free Euclidean vector of R d-1 . These approaches yield a triple representation of CoDa. Section 4 and 5 deal with statistical applications of such rescaled stick-breaking transformation for the study of the intra-dependence of CoDa. Section 4 introduces CoDapulas as the analogue of copulas for CoDa, opening the gates of the vast copula literature, tools and methodologies for CoDa. Several examples illustrates how copula models can easily be constructed for CoDa. In particular, completely monotone copulas give interesting complete dependence patterns for CoDa. Section 5 aims at studying intra-dependence of CoDa from the regression viewpoint. A basic example of a parametric regression model on real data set is given. Several extensions and alternatives are discussed. Eventually, we conclude in Section 6, with additional remarks about the choice of ordering of the components, mixability, and CoDa with zeroes.

2 The first stick-breaking view on CoDa: ordered points on the unit interval 2.1 The ordered stick-breaking view

Our approaches are based on the representation of CoDa as a normalised point p in the simplex (1). Instead of considering the (p i ), 1 ≤ i ≤ d, as primary parameters for p, one can consider the accumulated sums s = (s 0 , s 1 , . . . , s d ), defined as

s 0 = 0, s i = p i + s i-1 = i j=1 p j , i = 1, . . . , d -1 (2) s d = d i=1 p i = 1,
as an alternative system of coordinates of (1). Dropping the fixed values s 0 = 0 and s d = 1, this leads to a representation of (1) as

Σ d-1 := {(s 1 , . . . , s d-1 ) ∈ R d-1 : 0 ≤ s 1 ≤ . . . ≤ s d-1 ≤ 1}. (3)
(3) can be interpreted as iteratively breaking the unit stick [0, 1]: first, one picks some s 1 ∈ [0, 1], then, for the second step, one picks some s 2 in the remaining interval [s 1 , 1], and so on for s i to be picked in the interval [s i-1 , 1], i = 1, . . . , d -1. The process stops after d -1 steps. See Figure 1.

0 1 1 1 s 1 s 1 s 2 s 2 s 3 1 s 3 Figure 1: Stick-breaking of the unit interval Remark 1. 1. A normalized CoDa point p ∈ ∆ d-1 can be identified with a discrete probability measure µ p on R, p → µ p (.) := d i=1 p i δ xi (.),
where x 1 < . . . < x d ∈ R denotes (arbitrarily located) distinct components and δ x stands for the Dirac mass at x. The Coda p "forgets" about the locations x 1 , . . . , x d of µ p of the components, to only retains their probabilities p 1 , . . . , p d . As such, parametrization (2) of the simplex by the (s i ) interprets as characterizing the discrete distribution of a r.v. X ∼ µ p by its c.d.f. F (x) = P (X ≤ x) = y≤x P (X = x), while the parametrization by the (p i ) interprets as characterizing the distribution of X by its probability mass function P (X = x).

2.

(3) only uses the order structure of the interval [0, 1], and not the addition operation. This suggests that one can generalise the notion of simplex to arbitrary ordered space, endowed with a top (= 1) and bottom (= 0) element.

3. Instead of breaking the stick from the left to the right, i.e. putting s i ∈ [s i-1 , 1], for increasing i = 1, . . . , d -1, one can also consider a stickbreaking process from the right to the left, i.e. putting s i ∈ [0, s i+1 ] for decreasing i = d -1, . . . , 1. This corresponds to characterizing µ p by its survival function instead of its cumulative distribution function.

CoDa distributions via order statistics

Stick-Ordered distributions

This geometric view of CoDa as a set of ordered points on the unit interval suggests a natural connection with order statistics on the unit interval. This gives an easy way to build distributions on the simplex by taking as s i the order statistics of some u i distributed on the unit interval. More precisely, one can define a "Stick-Ordered" (SO) distribution on the simplex as follows:

Definition 2.1 (Stick-Ordered distribution for CoDa). Let u i ∼ F i , i = 1, . . . d- 1 be independent r.v. with (F 1 , . . . , F d-1 ) a set of univariate c.d.f.s on the unit interval. Set u (1) ≤ . . . ≤ u (d-1)
the corresponding order statistics. Eventually, define

s 0 = 0, s i = u (i) , i = 1, . . . d -1, s d = 1. Then, the CoDa point p ∈ ∆ d-1 corresponding to s in (3) is said to be (F 1 , . . . , F d-1 )- Stick-Ordered distributed, which is denoted by p ∼ SO(F 1 , . . . , F d-1 ).
In case

F i = F , p is said to be F -Stick-Ordered distributed, which is denoted by p ∼ SO(F ).
In other words, the Stick-Ordered distribution of p is the distribution of the spacings corresponding to the order statistics s. The latter has to be computed as the distribution of (possibly non-identically distributed) order statistics.

p ∼ SO(F 1 , . . . , F d-1 ) ⇐⇒      u i ∼ F i , (u 1 , . . . , u d-1 ) independent, s 0 = 0, s i = u (i) , s d = 1, i = 1, . . . , d -1, p i = s i -s i-1 , i = 1, . . . , d

Examples

Example 1. In particular, for F = U [0,1] the uniform distribution, (i.e. F (x) = x, 0 < x < 1), p ∼ SO(U [0,1] ) gives a uniform distribution on the simplex, as shown in Pyke (1965) equation (2.1). (See also S. S. [START_REF] Wilks | Order statistics[END_REF] equation ( 6)). More precisely, (p 1 , . . . , p d-1 ) has a density w.r.t. the d -1 dimensional Lebesgue measure given by

f (p1,...,p d-1 ) (x 1 , . . . , x d-1 ) = (d -1)! if x i ≥ 0, and d-1 i=1 x i ≤ 1 0 otherwise .
On the other hand, p has a singular distribution since

p d = 1 - d-1
i=1 p i , but its restriction to the hyperplane d i=1 p i = 1 admits a density w.r.t. the d -1 dimensional Lebesgue measure given by

f p (x) = (d -1)! if x i ≥ 0, and d i=1 x i = 1 0 otherwise ,
which is symmetric in (x 1 , . . . , x d ) (i.e. the (p i ) are exchangeable). One recognizes the Dirichlet Dir(1, . . . , 1; 1) distribution, see e.g. [START_REF] Rao | Incomplete Dirichlet integrals with applications to ordered uniform spacings[END_REF] or Samuel S. [START_REF] Wilks | Mathematical statistics[END_REF].

These stick-ordered distributions are useful for modeling purposes. They allow to construct CoDa models from classical distributions on [0, 1]. One can consider more examples with other distributions on the unit interval, like the Beta, the Kumaraswamy, (which is similar to the Beta distribution but leads to tractable formulas for the distribution of order statistics, see [START_REF] Jones | Kumaraswamy's distribution: A beta-type distribution with some tractability advantages[END_REF]), or those of [START_REF] Kotz | Beyond beta. Other continuous families of distributions with bounded support and applications[END_REF]. (More generally, any distribution on R can be mapped to a distribution with support included on the unit interval by applying a c.d.f to it). In some cases, analytical formulas can be obtained for the distribution of p, using known results on spacings [START_REF] Pyke | Spacings. (With discussion.)[END_REF]).

Example 2 (Kumaraswamy). The Kumaraswamy distribution [START_REF] Jones | Kumaraswamy's distribution: A beta-type distribution with some tractability advantages[END_REF]) U ∼ Kumaraswamy(α, β) has density

f (u) = αβu α-1 (1 -u α ) β-1 , 0 < u < 1, (4)
and cdf

F (u) = 1 -(1 -x α ) β , 1 < u < 1. (5)
Moreover, for i.i.d. u i ∼ F , i = 1, . . . , d -1, with density f , the marginal distribution of the spacings writes (see e.g. [START_REF] Pyke | Spacings. (With discussion.)[END_REF] p. 399)

f pi (x) = (d -1)! (i -2)!(d -1 -i)! (F (t)) i-2 (1 -F (x + t)) d-1-i f (t)f (x + t)dt, (6)
for 2 ≤ i ≤ d -1, and

f p1 (x) = f u (1) (x) = (d -1)f (x)(1 -F (x)) d-2 . ( 7
)
Applying formulas ( 6) and ( 7) to (4) and (5) leads to computable formulas. For example,

f p1 (x) = (d -1)αβx α-1 (1 -x α ) β(d-1)-1 , 0 < x < 1,
The Markov property of the order statistics can then be used to derive the joint distribution of p.

Generalised Stick-Ordered Distributions

One can also generalise the former definition (2.1) by taking dependent r.v. u i instead of independent ones.

Definition 2.2 (Generalized stick-ordered distribution for CoDa). For u = (u 1 , . . . , u d-1 ) ∈ [0, 1] d-1 with joint distribution function F u , the CoDa point p ∈ ∆ d-1 corresponding to s in (3) is said to be generalized-stick-ordered distributed with generator F u , which is denoted by p ∼ GSO(F u ), viz.

p ∼ GSO(F u ) ⇐⇒      (u 1 , . . . , u d-1 ) ∼ F u s 0 = 0, s i = u (i) , s d = 1, i = 1, . . . , d -1, p i = s i -s i-1 , i = 1, . . . , d
For non-identically distributed or dependent variables, one can compute them e.g. using results of [START_REF] David | Order statistics[END_REF] Chap. 5. [START_REF] Jaworski | On distributions of order statistics for absolutely continuous copulas with applications to reliability[END_REF], [START_REF] Rychlik | Distributions and expectations of order statistics for possibly dependent random variables[END_REF] (See also Balakrishnan N. ( 1998)). However, this often leads to intractable formulas. Nonetheless, it is easy to simulate samples from such distributions.

Example 3 (GSO with Gaussian copula generator). Consider, for d = 3, p ∼ GSO(F u ) with F U a bivariate Gaussian copula (hence with uniform marginals), with correlation ρ. Figure 3 shows ternary diagrams of scatterplots of samples of 1000 realisations of p = (p 1 , p 2 , p 3 ), with varying level of the dependence coefficient ρ. The value of ρ determines the behavior of the distribution of the CoDa element p and generates interesting patterns of dependence between the components. For ρ = 0, one generates a uniform distribution on the simplex (upper right panel). When ρ becomes negative (upper middle panel) and close to -1 (upper left panel), one obtains empirically a CoDa point s.t. p 2 ≈ 1 -2p 1 and p 1 ≈ p 3 . For ρ close to one (ρ = 0.99, lower right panel), the p 2 component is nearly zero and the CoDa point is nearly on the line p 1 = p 3 . 3 The rescaled stick-breaking view: unit cube geometry of CoDa points

Unit cube geometry for CoDa points by rescaling

The second approach we propose is based on a rescaled version of the iterative stick-breaking process of Figure 1: first, one picks some s 1 ∈ [0, 1], as previously. Then, one has to pick s 2 in the remaining interval [s 1 , 1]: in terms of spacings/lengths, the length of s 2 -s 1 = p 2 of the second stick [s 1 , s 2 ] has to be chosen relatively to the length 1 -s 2 of the remaining stick [s 2 , 1], see Figure 2. Similarly, the relative length s i -s i-1 = p i of the interval corresponding to the ith pick s i has to be chosen relatively to the length 1 -s i-1 of the remaining stick [s i-1 , 1]. Following the footsteps of [START_REF] Halmos | Random alms[END_REF], [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF] among others, it is therefore natural to introduce the transformation,

z 1 = s 1 = p 1 , z i = p i /(1 -s i-1 ), i = 1, . . . , d -1, (8)
z d = 1,
with inverse transformation,

p 1 = z 1 , p i = z i i-1 j=1 (1 -z j ), i = 2, . . . , d -1, p d = d-1 i=1 (1 -z i ). ( 9
)
By construction, the z i in ( 8) are in the unit cube, 0 ≤ z i ≤ 1 for i = 1, . . . , d -1, with degenerate z d = 1. Thus the transformation (8) turns the "akward" simplex ∆ d-1 into the unit cube [0, 1] d-1 (we drop z d = 1 as it is equal to 1). This leads to a free unit cube view of CoDa points. Remark 2. Historically, the transformation (8) was introduced by [START_REF] Halmos | Random alms[END_REF] for d = ∞, as a way to distribute gold dust to a countably infinite sequence of beggars, where each beggar receives in turn a fraction z i of the remaining gold. More generally, an infinite sequence (p 1 , p 2 , . . .) defined by ( 8) from a sequence of independent (z i ) ∈ [0, 1] is called a residual allocation model (RAM), see Feng ( 2021) for a review. The case when z i ∼ β(1, θ), θ > 0 is called a GEM distribution and was notably studied in population genetics by [START_REF] Griffiths | Exact sampling distributions from the infinite neutral alleles model[END_REF], [START_REF] Engen | A note on the geometric series as a species frequency model[END_REF][START_REF] Mccloskey | A Model for the Distribution of Individuals by Species in an Environment[END_REF], see e.g. [START_REF] Ewens | Population genetics theory-the past and the future[END_REF]. RAM also appears in Bayesian statistics in connection with the Dirichlet distribution, as the weights of random measures ), where (ξ k ) are i.i.d. and independent of (p k ), see Samuel S. [START_REF] Wilks | Mathematical statistics[END_REF]) p. 178, Ferguson (1973), [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF]. See also Section 3.3 below.

P(.) = ∞ k=1 p k δ ξ k (.

Neutrality and complete neutrality

This interpretation of CoDa points as a relative/proportional iterative stickbreaking process leads to the concept of neutrality, introduced by [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF], which is relevant for the analysis of CoDa. In short, it is a sort of intra-independence concept for a random composition p.

More precisely, neutrality is motivated by the following: if one wants to check whether the first proportion p 1 has an influence on the remaining subcomposition (p 2 , . . . , p d ), the latter has to be rescaled by the remaining mass 1 -p 1 , in order to be a proper normalized CoDa point. One thus has to check for the stochastic influence of p 1 on

p 2 1 -p 1 , . . . , p d 1 -p 1 , (10)
and if p 1 is independent of the latter rescaled subcomposition, one can eliminate p 1 from the analysis of p. Therefore, [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF] defines neutrality as follows: p 1 is said to be neutral if p 1 is independent of (10): p 1 does not influence the manner in which the remaining proportions (p 2 , . . . , p d ) relatively divide the remainder of the unit interval.

A generalisation of neutrality to a vector [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF]). A further generalization of neutrality is complete neutrality: if the z = (z 1 , . . . , z d-1 ) of ( 8) are mutually independent, then the corresponding p is said to be completely neutral, or equivalently (Theorem 2 in [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF]) if p j is neutral for all 1 ≤ j ≤ d -1. These concepts of neutrality are helpful for constructing completely neutral distributions on the simplex: start with mutually independent z i 's each having a specified distribution on [0, 1], and invert (8) to obtain a completely neutral distribution on the simplex. In particular, [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF] construct a generalisation of the Dirichlet distribution from independent z i ∼ β(a i , b i ). [START_REF] Ng | Dirichlet and related distributions[END_REF] Theorem 2.2 use the transformations ( 8) and ( 9) to obtain stochastic representations of the Dirichlet distribution D(a) from independent

p k := (p 1 , . . . , p k ), k < d is: (p 1 , . . . , p k ) is a neutral vector if it is independent of ( p k+1 1 -s k , . . . , p d 1 -s k ). Thus, if p j is neutral for j = 1, . . . k, then z k := (z 1 , . . . , z k ) is mutually inde- pendent (Theorem 1 in
z i ∼ β(a i , d k=i+1 a k ), i = 1, . . . , d -1.

Conditional probability interpretation of the rescaled stick-breaking approach and connection with Bayesian priors

The rescaled weights z i interpret as conditional probabilities. The stick-breaking construction appears in the construction of the (finite-dimensional) Dirichlet distribution of [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], see e.g. [START_REF] Ghosal | Fundamentals of nonparametric Bayesian inference[END_REF] p. 30. The latter is used for constructing a prior on a discrete distribution in Bayesian statistics. It is defined as follows: in order to randomly distribute a total mass 1, identified with the unit interval, to the first d integers 1, 2, . . . , d, the stick is first randomly broken by a r.v. 0 ≤ Z 1 ≤ 1, and mass Z 1 is assigned to 1. The remaining mass is 1 -Z 1 and the stick [Z 1 , 1] is broken into two new pieces of relative length Z 2 and 1 -Z 2 , for some 0 ≤ Z 2 ≤ 1. Mass (1 -Z 1 )Z 2 is assigned to the point 2, and the remaining stick has remaining mass (or length) (1-Z 1 )(1-Z 2 ). Iterating, one has defined a random distribution (i.e. a Markov kernel), with values j = 1, . . . , d and (random) probabilities given by (9). Each p i is the probability assigned to i, conditionally on the previous probabilities assigned to the j < i. Indeed, if one denotes by ζ the r.v. with values in 1, . . . , d and (random) probabilities given by ( 9), i.e. s.t.

P (ζ = i) = p i = Z i i-1 j=1 (1 -Z j ), i = 1, . . . , d, Then, Z i = P (ζ = i|ζ ≥ i).
On the other hand, the complete neutrality property expresses the idea that these Z i (or equivalently these conditional probabilities) are chosen independent. In particular, if (Z 1 , . . . , Z d ) are independent with [START_REF] Ghosal | Fundamentals of nonparametric Bayesian inference[END_REF] Corollary G.5.

Z i ∼ β(α i , j>i α j ), then p is Dirichlet Dir(k; α 1 , . . . , α d ) distributed, see

A triple representation of CoDa

One thus has a triple representation of the simplex / of normalized CoDa points: the simplex can be represented as ∆ d-1 with its sum constraint, as the ordered set of points Σ d-1 on the unit interval, or a free cube [0, 1] d-1 via its rescaled representation in z coordinates. Figure 5 shows the different representations as well as the transformations between them (where the arrows between the free unit cube [0, 1] d-1 and Σ d-1 are obtained by composition of the previous transformations). Note that the ordered and rescaled representations are not canonical, as they depend on the order of enumeration of the components of p. . This is also illustrated in Figure 5, (lower-left). The figure shed lights on some results and representations of order statistics and constructions of Dirichlet distribution. For example, it is well-known (see [START_REF] Sukhatme | TESTS OF SIGNIFICANCE FOR SAM-PLES OF THE x2-POPULATION WITH TWO DEGREES OF FREE-DOM[END_REF], [START_REF] Rényi | On the theory of order statistics[END_REF], [START_REF] Pyke | Spacings. (With discussion.)[END_REF]) that the order statistics and spacings of i.i.d (u i ) r.v. uniformly distributed on [0, 1], have a representation as a ratio of (sums of) exponential r.v.: Take y i ∼ Exp(1) i.i.d. in Figure 5, then normalisation by the sum gives the p i which corresponds to spacings, and summing this spacings give the order statistics

Relative/conditional representation z ∈ [0, 1] d-1 Sum representation p ∈ ∆ d-1 ⊂ [0, 1] d Ordered representation s ∈ Σ d-1 ⊂ [0, 1] d-1 y ∈ R d ≥0 u ∈ [0, 1] d-1
s i = u (i) , (u (i) , i = 1, . . . , d -1) d = j≤i y j d j=1 y j , i = 1, . . . , d -1 and (p i , i = 1, . . . , d) d = y i d j=1 y j , i = 1, . . . , d
Also, for y i ∼ γ(α i ) Gamma distributed, Figure 5 allows to explain and visualize the difference between the Dirichlet distribution on ∆ d-1 and its ordered version on Σ d-1 , see Samuel S. Wilks (1962) p. 178, 182 and 238. 3.5 From the unit cube to the free Euclidean space R d-1

If p has no zero components, viz. 0 < p i < 1 for all 1 ≤ i ≤ d, then p is sent to the interior (0, 1) d-1 of the unit cube [0, 1] d-1 by the rescaled stick-breaking transformation (8). In turn, one can then map the open unit-cube representation z ∈ (0, 1) d-1 of the Coda element p ∈ ∆ d-1 to a point ξ = (ξ 1 , . . . , ξ d-1 ) ∈ R d-1 by applying an increasing 2 continuous transformation q : (0, 1) → R to each component z i of z, viz.

ξ i = q(z i ), 1 ≤ i ≤ d -1.
See Figure 6. Examples of q which come to mind include the probit, logit transform, or any quantile function of a distribution on R with positive density (hence the notation q).

p ∈ ∆ d-1 , p > 0 z ∈ (0, 1) d-1 ξ ∈ R d-1 rescaled stick-breaking q q -1 Figure 6: Free Euclidean representation ξ ∈ R d-1 of positive CoDa.
This gives an interesting alternative to the vector space representation provided by Aichison's log-ratio transforms. This variant of the z representation allows to apply standard multivariate analysis techniques designed for Euclidean vectors to CoDa. For example, one can apply classical Principal Component Analysis to the transformed variables ξ for exploratory data analyses of CoDa. Clustering algorithms (i.e. k-means) can be applied on the ξ-representation of CoDa, without further ado. On the modeling side, any classical multivariate distribution for ξ on R d-1 gives, by back-transformation, a corresponding CoDa distribution for p ∈ ∆ d-1 . Potential applications are numerous.

Application of unit cube geometry: Copulas for CoDa

In the complete neutrality view of [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF], each rescaled component z i does not influence the remaining ones. If one moves out of independence, one can generalize in several directions.

CoDapulas: copulas for CoDa

As first generalization, one can construct copula models [START_REF] Nelsen | An introduction to copulas[END_REF]) for CoDa: instead of taking independent z i , one can specify a joint distribution for z by a set of marginals (F zi ), i = 1, . . . , d -1, (each with support the unit interval) and a copula C, viz.

z ∼ C(F z1 , . . . , F z d-1 ).

2 or, more generally, a strictly monotone continuous function.

By back transformation (9), this allows to define general distributions for CoDa points from the specification of a copula and the marginal distributions of the z.

A probabilistic construction of this specification is as follows: let v ∈ [0, 1] d-1 be distributed according to a copula function C, i.e. a multivariate distribution with uniform marginals, and let

Q zi = F -1 zi : [0, 1] → [0, 1], i = 1, . . . , d - 1, be univariate quantile functions with range [0, 1]. Set z i = Q zi (v i ), i = 1, . . . , d -1. Then z = (z 1 , . . . , z d-1 ) ∈ [0, 1] d-1
has copula function C and marginal distributions (F zi ). By back transformation, p ∈ ∆ d-1 is a CoDa point whose distribution is uniquely specified by C an the set (F zi ) of marginal cdfs. Explicitly, (9) yields

p 1 = Q z1 (v 1 ), ( 11
)
p i = Q zi (v i ) i-1 j=1 (1 -Q zj (v j )), i = 2, . . . , d -1, ( 12
)
p d = d-1 i=1 (1 -Q zi (v i )). (13)
Conversely, given a CoDa point p, one can estimate and study its intradependence through the copula of its z-representation: one first transforms p into z by transformation (8), and then standardize the marginals z i to the uniform distribution. The latter operation is obtained, when z is continuous, by the marginal probability integral transforms,

v i := F zi (z i ), i = 1, . . . , d -1,
where F zi is the c.d.f. of z i . Then, v = (v 1 , . . . , v d-1 ) has uniform marginals, i.e. has a copula distribution. See Figure 7.

p ∈ ∆ d-1 z ∈ [0, 1] d-1 v ∼ C Fz i Qz i
Figure 7: Probabilistic construction of a CoDapula (In the non-continuous case, the standardization is obtained by using the marginal distributional transforms instead. The latter is defined as

F zi (x, η) := P (z i < x) + ηP (z i = x), η ∈ [0, 1]. ( 14
)
Then, v is obtained by setting

v i := F zi (z i , η i ), i = 1, . . . , d -1,
where (η i ) is a sequence of i.i.d. uniformly distributed on [0, 1] randomizers, independent of z, see e.g. Ludger [START_REF] Rüschendorf | On the distributional transform, Sklar's theorem, and the empirical copula process[END_REF], [START_REF] Faugeras | Markov morphisms: a combined copula and mass transportation approach to multivariate quantiles[END_REF]).

Let us give a fancy name to the copula of a CoDa point.

Definition 4.1 (CoDapula). Let p ∈ ∆ d-1 a random CoDa point, and z ∈ [0, 1] d-1 be it rescaled stick-breaking representation (8). Then, a CoDapula of (the distribution of ) p is a copula of (the distribution of ) z. In other words, a CoDapula C of p is the distribution of v in the construction of Figure 7.

Thanks to Sklar's Theorem [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF]), a CoDapula always exists. It is unique if z is continuous (see e.g. [START_REF] Nelsen | An introduction to copulas[END_REF]). Definition 4.1 depends on the ordering of the components 1, . . . , d. Hence, a CoDapula of p depends on a permutation π of {1, . . . , d}. Hence, in full rigor, one should have defined a notion of π-CoDapula to stress the dependence on π. We have chosen not to in order to simplify notations. The choice of the ordering, i.e. of π, may depend on the application in view, and will be discussed in Section 6.2.

By (11), the first component p 1 has the same distribution as z 1 , and thus is completely specified by the first marginal distribution function F z1 , (equivalently, quantile function Q z1 ). Note that the marginal distributions of the remaining components p 2 , . . . , p d depend on both the CoDapula and the marginal distributions: this is in contrast with the copula approach for classical Euclidean vectors. Nonetheless, at the z level, one has the the classical copula separation of a multivariate distribution into its marginal distributions F z1 , . . . F z d-1 and the dependence structure embodied in the copula function C.

Examples and numerical illustrations

We illustrate in Figures 8 and9 some CoDa distributions which can be obtained using the specification by a CoDapula of Definition 4.1 and marginal quantile functions for z. In Figure 8, the copula of z is an Ali-Mikhail-Haq copula with parameter α = 0.91, and the marginal distributions are Beta and uniform, F z1 ∼ β(1/4, 2) and F z2 ∼ U [0,1] , while in Figure 9, the CoDapula is a Gumbel-Hougaard copula with parameter θ = 7, and same marginals as in Figure 8. The left panels show scatterplots in the z domains, and the right panels the corresponding ternary plots in the simplex CoDa space for p. Figure 8 give an example of mild dependence with CoDa points spreading above the p 1 = 0 level. As the Gumbel copula approaches the comonotonicity copula as θ → ∞ (while θ = 1 yields the independence copula), the value θ = 7 in Figure 9

Complete dependence of CoDa

The concept of CoDapula opens the gates of the vast copula literature and modeling methodology to CoDa. This is useful to study the intra-dependence of CoDa. The independence copula for z means that p is completely neutral.

At another extreme, complete dependence at the z level induces a specific dependence pattern at the p level, as is shown in the next two examples.

Comonotone CoDapula

Comonotonicity is an extreme form of dependence structure for Euclidean vectors that describes the strongest positive dependence. A comonotone vector is characterised by having as copula the comonotone copula M (x 1 , . . . , x d-1 ) = min(x 1 , . . . , x d-1 ). The comonotone copula corresponds to the distribution of the vector v = (v, . . . , v) ∈ [0, 1] d-1 , with a single v ∼ U [0,1 . In other words,

P (v ≤ x) = P (v ≤ x 1 , . . . , v ≤ x d-1 ) = min(x 1 , . . . , x d-1 ), x ∈ [0, 1] d-1 .
Applied to CoDa, the corresponding z thus writes

z = (Q z1 (v), . . . , Q z d-1 (v)),
where Q zi : [0, 1] → [0, 1] are given quantile functions. This gives, as corresponding p, Coda with components

p 1 = Q z1 (v) p i = Q zi (v) i-1 j=1 (1 -Q zj (v)), i = 2, . . . , d -1.
Example 4 (Comonotone CoDapula, d = 3). For example, for d = 3, one gets

z = z 1 z 2 = Q z1 (v) Q z2 (v)
which translates into

p 1 = Q z1 (v) p 2 = Q z2 (v)(1 -Q z1 (v)) p 3 = 1 -p 1 -p 2 = (1 -Q z1 (v))(1 -Q z2 (v)) (15)
Thus, p 1 is an increasing function of v, p 3 is decreasing, while p 2 switches direction of variation w.r.t v.

Figure 10 shows, for d = 3, the CoDa p corresponding to the comonotone copula, with uniform quantile functions at the z level, viz. The distribution of p is singular, as each component p i is a deterministic function of v ∼ U [0,1] : p lies on the curve shown in the ternary plot. This implies that each pair of components (p i , p j ), 1 ≤ i = j ≤ 3 are totally dependent, i.e. lie on a curve. Figure 11 shows the resulting complete dependence between each pairs of components: (p 1 , p 3 ) are counter-monotone (middle), while (p 2 , p 3 ) (right) is comonotone. (p 1 , p 2 ) (left) switches its sense of variation, being first comonotone, then countermonotone. Note that p 2 also has a limited range of variation p 2 ∈ [0, 1/4].

Q z1 (v) = Q z2 (v) = v in (15), so that p 1 = v, p 2 = v(1 -v), p 3 = (1 -v) 2 .

Counter-monotone CoDapula

Counter-monotonicity is the antithesis of comonotonicity. Note that this notion is well-defined only in two dimensions. We thus restrict our discussion to the case d = 3. The bivariate counter-monotone copula

W (x 1 , x 2 ) = max(x 1 + x 2 -1, 0) is stochastically realized by the vector v = (v, 1 -v), where v ∼ U [0,1]
. This gives a corresponding Coda

p =   p 1 p 2 p 3   =   Q z1 (v) Q z2 (1 -v)(1 -Q z1 (v)) (1 -Q z1 (v))(1 -Q z2 (1 -v))   . ( 16
)
The following example illustrates the case when the quantile functions are the uniform ones.

Example 5 (Counter-monotone CoDapula).

For Q z1 (v) = Q z2 (v) = v, (16) gives p 1 = v, p 2 = (1 -v) 2 , p 3 = v(1 -v).
One thus gets the same parametrization at the CoDa level as in the comonotone case of Example 4, but with the roles of p 2 and p 3 exchanged, see Figure 12. This now translates at the CoDa level into complete dependence between pairs of components, as shown in Figure 13. Notice, however, that the dependence pattern is not the symmetric of the comonotone case of Figure 11: (p 1 , p 2 ) are now counter-monotone, whereas both (p 1 , p 3 and (p 2 , p 3 ) change their direction of variation. Hence, with a counter-monotone CoDapula, only one pair is monotone dependent (viz. (p 1 , p 2 ) counter-monotone), whereas with a comonotone CoDapula, two pairs were monotone dependent (viz. (p 1 , p 3 ) counter-monotone, and (p 2 , p 3 ) comonotone). 5 Application of unit cube geometry: Regression models for CoDa

As a second possible generalization, the rescaled stick-breaking approach (8) can be useful for the intra regression analysis of a CoDa component w.r.t the others. The basic idea is to construct regression models in the z coordinates to iteratively explain one z i component in terms of the other z j . Indeed, the transformation ( 8) is reminiscent of Rosenblatt's generalization of the quantile transform by successive conditioning and the regression representation of a random vector, which we recall now.

Regression representation of an Euclidean random vector

Let X = (X 1 , . . . X k ) ∈ R k be a vector with joint c.d.f. F . If one can transform X into a sequence 1 , . . . , k of independent, identically distributed r.v., with a prescribed distribution λ (say, uniform on [0, 1]), then, one can argue that the distribution of X has been successfully modeled: the transformation

(X 1 , . . . , X k ) φ → ( 1 , . . . , k )
has stripped X of all its stochastic variability and dependence and turned it into white noise. The function φ effectively models the distribution F of X. Rosenblatt (1952)'s transform and its generalizations (see Ludger Rüschendorf ( 2009)) achieves such a reduction: Denote by F i|i-1,...,1 the conditional c.d.f. of X i given (X i-1 , . . . , X 1 ), i = 2, . . . , k, with F 1 the (marginal) cdf of X 1 . Rosenblatt (1952)'s transform is then defined by = ( 1 , . . . k ) with

1 := F 1 (X 1 ) i := F i|i-1,...,1 (X i |X i-1 , . . . , X 1 ), i = 2, . . . , k.
Under an assumption of continuity3 of the successive conditional c.d.f. F i|i-1,...,1 , Rosenblatt's transform turn the vector X into a vector of i.i.d. U [0,1] components (see Ludger [START_REF] Rüschendorf | On the distributional transform, Sklar's theorem, and the empirical copula process[END_REF]).

Conversely, starting from a vector ∼ λ k and applying the successive (conditional) quantile functions F -1 i|i-1,...,1 , viz.

X 1 := F -1 1 ( 1 ) X i := F -1 i|i-1,...,1 ( i |X i-1 , . . . , X 1 ), i = 2, . . . , k, (17)
one obtains a vector X with the desired joint c.d.f. F . Each equation ( 17) interprets as a nonlinear regression equation of X i , given its past covariates X j , j < i, with error/noise/innovation i . This gives a regression representation of X, according to Ludger [START_REF] Rüschendorf | On the distributional transform, Sklar's theorem, and the empirical copula process[END_REF], where ( 17) is a (triangular) stochastic representation of the successive predictive distributions P Xi|Xi-1,...,X1 . In ( 17), the distribution of ( 1 , . . . , d ) is purely conventional, the only constraint is that it be absolutely continuous (so that any distribution of X can be obtained from it by mapping and not by Markov kernels, see [START_REF] Faugeras | Markov morphisms: a combined copula and mass transportation approach to multivariate quantiles[END_REF]). In particular, one can choose the more familiar Gaussian white noise framework by setting

i = φ( i ),
where 1 , . . . , k are i.i.d. standard Gaussian N (0, 1), and φ is the c.d.f. of the N (0, 1) distribution.

The regression representation (17) is the general, exact, nonlinear form of a regression model. In particular, if the conditional quantile functions F -1 i|i-1,...,1 are linear, one obtains the classical linear model (albeit with uniform noise), viz.

X i = a i,1 X 1 + . . . + a i,i-1 X i-1 + i ,
where a i,1 , . . . , a i,i-1 are parameters.

Parametric internal regression models for CoDa

This suggests to make use of this regression representation to construct triangular regression models on the z representation of CoDa, by applying the transformation (17) to the z i of (8) instead of the X i : each z i is explained in terms of the previous z j , j < i, and some extraneous randomness i ∼ λ, for i = 1, . . . , d -1. Then, by back transformation (8), one obtains a (possibly nonlinear) regression model for the original p i , which can be used for internal prediction of a component in term of the others. More precisely, let ∼ λ d-1 be a vector of uniform noise on [0, 1] d-1 . Then, a general nonlinear triangular regression model for the z writes

z 1 = φ 1 ( 1 ) z i = φ i ( i , z i-1 , . . . , z 1 ), i = 2, . . . , d -1, (18)
where

φ i : [0, 1] i → [0, 1] are s.t. i → φ i ( i , z i-1 , . . . , z 1 ) is non-decreasing, left-continuous, with φ i (0, z i-1 , . . . , z 1 ) = 0, φ i (1, z i-1 , . . . , z 1 ) = 1. (i.
e. the φ i satisfy the properties of univariate quantile functions).

For example, a Gausssian (partially )linear triangular model can be obtained by specifying the error distribution as standard multivariate Gaussian ∼ N (0, I d-1 ), and the z i as

z 1 = Φ( 1 ) z i = Φ(a i,1 z 1 + . . . + a i,i-1 z i-1 + i ), i = 2, . . . , d -1 (19)
where Φ, the cdf of the standard univariate Gaussian distribution, is applied to ensure that z i ∈ [0, 1]. More generals models can be constructed via Generalized Linear Models, see e.g. [START_REF] Mccullagh | Generalized linear models[END_REF], and more specifically, [START_REF] Bonat | Regression models with responses on the unity Interval: Specification, estimation and comparison[END_REF] for data on the unit interval [0, 1].

Example: agriculture data.

We provide below a basic example of the construction of a parametric internal regression model in the z space for CoDa, illustrated on a real dataset. The data is taken from the example datasets accompanying Mathematica's (Wolfram Research, Inc. (n.d.)) "TernaryListPlot" command. It gives the raw amount of fertilizers (Nitrogen-Potassium-Phosphate) in a time series from 1960 to 2015. The scatter plots, both at the z level (left panel), and at the compositional level in the ternary plot (right panel) in Figure 14 show a cyclic pattern in the composition of fertilizers. The sinusoidal shape of the transformed data at the z level (left panel) suggests the following model,

z 2 = a 0 + a 1 cos(ωz 1 + φ) + ,
where a 0 , a 1 , ω, φ are parameters to be estimated and the random error. Nonlinear least squares (command "NonLinearModelFit" in Wolfram Research, Inc. (n.d.)) gives as fitted model, p 3 ≈ 1 -p 1 -p 2 ≈ (1 -p 1 )(0.5 -0.047 cos(-21.86 p 1 + 8.76)).

The resulting fitted curve is shown in orange in Figure 14, with extrapolated values on the full range

z 1 = p 1 ∈ [0, 1].

Extensions and alternatives

As shown in the previous example of Section 5.3, the rescaled stick-breaking transformation (8) reduces internal regression analysis of CoDa to classical regression analysis of vector data. Hence, all classical multivariate regression analysis techniques apply to CoDa, in their transformed z representation of the free unit cube. For space constraints, we limited ourselves in the example of Section 5.3 to a very basic illustration with a parametric regression model. Let us thus briefly mention some extensions and alternatives:

• A nonparametric alternative to the above intra-parametric models is to directly start from ( 17) and estimate the conditional distributions of z i given (z i-1 , . . . , z 1 ), or some functional thereof, via some nonparametric estimate. For example, one can look for the mean of these conditional distributions, and estimate the regression function of E[z i |z i-1 , . . . , z 1 ] by a Nadara-Watson, spline, or local polynomial estimator.

• Many applications are interested in explaining/predicting a CoDa point p w.r.t. some covariates X ∈ R k , i.e. in studying the conditional distribution of p|X = x. This can be done via the rescaled z representation (8) by performing a regression of z w.r.t. the covariates X. As said before, one must ensure that the constraint 0 ≤ z ≤ 1 is fulfilled. This can be achieved by a link function which entails the correct normalisation or by mapping z ∈ [0, 1] d-1 into ξ ∈ R d-1 , by using the device explained in Section 3.5. In parametric regression models, one can incorporate these external covariates X by making the parameters of p in its rescaled z representation ( 8), (i.e. the functions φ i in (18) or the coefficients a j,i in ( 19)) as functions of the covariates X. Once a regression model/ or a nonparametric estimate for z given X has been computed, one back transforms the predicted values ẑ into predicted values of p, via the inverse transformation (9).

• Our focus in this paper is on internal dependence analysis of CoDa. However, one can also easily envision external regression analysis such as CoDa to CoDa or CoDa to vector, by conducting a similar analysis in the corresponding z space for the CoDa input/output variables considered.

• In addition to explicit regression models, one can also assess quantitatively neutrality of z 1 by the strength of the regression dependence between z 1 and the remaining z j components of p, for j > 1. This quantification can be achieved through multivariate asymmetric correlation coefficients, like the recent Griessenberger, Junker, and Trutschnig (2022)'s ζ 1 (X, Y ) or Azadkia and Chatterjee (2021)'s T (X, Y )4 , the latter being a multivariate extension of the bivariate measure ξ of Chatterjee (2020). These coefficients quantifies the extent of regression dependence of a univariate random variable Y on a k-dimensional random vector X = (X 1 , . . . , X k ): they are equal to 0 in case of independence, and equal to 1 if Y is measurable function of X. Applied to our context, one can thus quantify the amount of (non-)neutrality of p 1 by computing ζ 1 ((z 2 , . . . , z d ), z 1 ) or T ((z 2 , . . . , z d ), z 1 ).

6 Conclusion and further remarks

Conclusion

We have proposed two related transformations for CoDa based on stick-breaking processes. The first one represents a CoDa point as a set of ordered values on the unit interval, whereas the second one, which originates from [START_REF] Halmos | Random alms[END_REF] and [START_REF] Connor | Concepts of Independence for Proportions with a Generalization of the Dirichlet Distribution[END_REF], removes the unit-sum constraint of the simplex representation and turns a CoDa point into a free vector of the unit cube. Both approaches are useful to construct distributions for CoDa from multivariate distributions of Euclidean vectors. The second approach appears most promising as it allows for a reduction of CoDa points to classical multivariate vectors and thus allows the use of well-established multivariate analysis techniques and models to be directly transferred to CoDa. Such an approach is an interesting alternative to the classical log-ratio coordonatizations techniques of Aitchison and his followers. In particular, we introduced the concept of a CoDapula, a copula for CoDa, and showed how to study the intra-dependence of CoDa with such copulas or via regression techniques.

For length reasons, we have barely scratched the surface of statistical applications based on these transformations. Much more needs to be done to explore the potentialities of the proposed approaches. Being now reduced to vector data, important statistical models and tools like generalised linear models, graphical models, vines/factor copulas, clustering, Principal Component Analysis, nonparametric and semi-parametric techniques, etc. are now at the disposal of the Statistician and beg for their application to CoDa. Let us close the article with some further remarks.

Choice of the ordering of the components

The only possible issue of the transformations (2) and ( 8) is the lack of symmetry w.r.t. the components, as they depend on the ordering of the components 1, . . . , d of the composition. The question thus arises which ordering is most adequate. Several possibilities can be envisioned.

• A first possibility is to let the Statistician decide for himself. This is similar to the "working-in-coordinates" principle in classical log-ratio CoDa analysis: the statistical model is extrinsic and built w.r.t. a given coordinate frame (here, the ordering chosen), and is eventually mapped back to the original simplex. This was the approach chosen in the example of Section 5.3: the sinusoidal regression model in the z space is mapped to the CoDa simplex space and then gives a model which explains/predicts how the remaining components p 2 , p 3 are driven by p 1 .

• The order of the components can be dictated by the type of application in view. For example, in a general regression model Y = r(X, ), there is a natural asymmetry in the vector (X, Y ) between the dependent/predicted variable Y and the independent/predictor variables X: one wants to explain/predict Y from X (with noise ). There is also an asymmetry in the rescaled stick-breaking transformation (9), between the first component p 1 and the remaining ones p 2 , . . . , p d : the first component is identical in the simplex space p representation as in the unit cube z representation, i.e. p 1 = z 1 , whereas the remaining components p 2 , . . . , p d depend on several of the z i . (p i is a function of z j , for 1 ≤ j ≤ i). Thus, z 1 is directly interpretable as one original component and a statistical analysis of z 1 translates into a statistical analysis of the first component p 1 . So, if one is interested in evaluating how a specific component is influenced by the remaining parts, it is sensible to take this component as first one: a regression model z 1 = r(z 2 , . . . , z d , ) at the z level, following the methodology explained in Section 5, with z 1 = p 1 as predicted variable, directly gives a regression model of the first component p 1 in terms of the remaining (rescaled) components, viz. p 1 = r(z 2 , . . . , z d , ).

• One can also envision a data-dependent choice of the ordering of the components: the basic idea of the z transformation is to transform the study of the non-neutrality of the constrained components p i of the composition into a the study of the dependence of the free z i . Thus, it would make sense to order the components by decreasing order of non-neutrality/dependence with the remaining composition. If, w.l.o.g, the first component p 1 is most dependent with the remaining composition ( p2 1-p1 , . . . , p d 1-p1 ), it means that p 1 is the main factor explaining the remaining composition. Having isolated such a component, one can then look within the closed remaining composition ( p2 1-p1 , . . . , p d 1-p1 ) of size d -1, which component is most dependent with the closed subcompositions of size d -2. The process is then iterated, yielding an ordering of the components. In practice, such evaluation of the dependence between a component and a subcomposition can be performed using the estimators of the asymmetric regression dependence coefficients ζ 1 (X, Y ) of [START_REF] Griessenberger | On a multivariate copula-based dependence measure and its estimation[END_REF] or T (X, Y ) of [START_REF] Azadkia | A simple measure of conditional dependence[END_REF], mentioned in Section 5.4. This gives the following algorithm: For a composition p of size d, 1. Select j s.t. T (W j , p j ) (or ζ 1 (W j , p j )) is maximum, where W j := . . . , p k 1-pj , . . . • One can also mix the above approaches, e.g. select as first component the one the Statistician is interested in explaining/predicting, and select the remaining ones in a data-dependent manner.

Connection with mixability: existence of CoDa distributions with given marginals

The proposed approaches were constructive and gave explicit representations of distributions of CoDa points. In case of the ordered approach, the stick-ordered distributions of Definition 2.1, were parametrized by d -1 univariate marginals. Similarly, in the rescaled approach, the distributions are parametrized by either a d -1 dimensional copula and d -1 marginals, or a set of d -1 conditional distributions. These approaches were helpful in constructing distributions for d-dimensional CoDa points. A converse issue is to enquire for the existence of a d-dimensional CoDa distribution with a given set of d marginal distributions. This question is related to the notion of joint mixability, which is a notion mainly investigated in the risk theory literature (See the survey by R. [START_REF] Wang | Current open questions in complete mixability[END_REF]). This connection between mixability and distributions for CoDa does not seem to have been made beforehand by the CoDa community.

The definition of joint mixability (B. Wang and R. Wang (2016)) is as follows:

Definition 6.1. An d-tuple of probability distributions on R, (F 1 , . . . , F d ) is jointly mixable if there exist d random variables X 1 ∼ F 1 , . . . , X d ∼ F d such that X 1 + . . . + X d =: K is almost surely a constant.

Hence, the question of existence of a d-dimensional CoDa distribution with given marginals is a special case of mixability with K = 1. [START_REF] Gaffke | On a class of extremal problems in statistics[END_REF] Theorem 5 give a necessary and sufficient condition. Necessary conditions are given in Theorem 2.1 in B. Wang and R. Wang (2016), and sufficient conditions are given in Theorems 3.1, 3.2, and 3.4 for uniform, monotone, and symmetric-unimodal densities, respectively.

CoDa with some zero components

The well-known additive log-ratio transformation (alr) of J. [START_REF] Aitchison | The statistical analysis of compositional data[END_REF] y i := log(p i /p d ), i = 1, . . . , d -1 and its variants clr and ilr [START_REF] Greenacre | Compositional data analysis in practice[END_REF], Pawlowsky-Glahn, Egozcue, and Tolosana-Delgado ( 2015)) turns a CoDa point p as a vector element of R d-1 . However, it is undefined if p has some zero components and special treatments of the zero components are required, like amalgamation of the finer parts, replacement of the zeroes with small values, or treatment of the zero observations as outliers. However, these treatments have an ad-hoc character and are not completely satisfactory from a modeling perspective. Hence, the statistical literature based on log-ratio analysis usually enforces an assumption of strict non-negativity of the components, which limits its scope of application.

To the contrary, the stick-breaking representations do not require such an assumption: in the ordered representation (3), zeroes translates into ties in s, while in the rescaled z representation, zeroes of p translates into z being sent to the boundary of [0, 1] d-1 . Distributionally, this means that the (u i ) in Definition 2.1 of the ordered representation (3) have common discrete components in their distributions. For the rescaled z representation (8), zeroes of p translates distributionally into z having a singular component on some faces of [0, 1] d-1 in the Lebesgue decomposition of the probability measure of z. One can therefore use mixed/general distributions to model such CoDa points with possibly zero components.

Figure 2 :

 2 Figure 2: Stick-Ordered distribution SO, obtained from the order statistics u (i)

Figure 3 :

 3 Figure 3: Ternary plots of Generalized-Stick-Ordered distribution for d = 3, with Gaussian copula generator, with varying correlation coefficient ρ. From left to right and up to down: ρ = -0.99, ρ = -0.8, ρ = 0, ρ = 0.8, ρ = 0.99.

Figure 4 :

 4 Figure 4: Rescaled stick-breaking

Figure 5 :

 5 Figure 5: Representations of the simplex

Figure 8 :Figure 9 :

 89 Figure 8: CoDa with AMH CoDapula (α = 0.91), and β(1/4, 2), U [0,1] marginal distribution functions: scatter plot at the z level (left) and ternary scatter plot for the resulting p (right), d = 3.

Figure 10 :

 10 Figure 10: Ternary plot of a CoDa with comonotone CoDapula and uniform quantile functions, d = 3, with barycentric axes p 1 (red), p 2 (green), p 3 (blue).

Figure 11 :

 11 Figure 11: Complete dependence between pairs of components of a CoDa with comonotone CoDapula and uniform quantile functions, d = 3: (p 1 , p 2 ) (left), counter-monotone (p 1 , p 3 ) (middle), comonotone (p 2 , p 3 ) (right)

Figure 12 :

 12 Figure 12: Ternary plot of a CoDa with countermonotone CoDapula and uniform quantile functions, d = 3, with barycentric axes p 1 (red), p 2 (green), p 3 (blue).

z 2 ≈

 2 0.5 + 0.047 cos(-21.86 z 1 + 8.76), with sums of squares of 13.43 for the model and 0.0075 for the error. In terms of the original CoDa variables p, this yields a predictive model of the components in term of the first one p 1 , p 1 = p 1 p 2 ≈ (1 -p 1 )(0.5 + 0.047 cos(-21.86 p 1 + 8.76))

Figure 14 :

 14 Figure 14: Internal parametric regression model for agriculture data. Scatter plots (blue points) and fitted sinusoidal model (orange line) in the z space (left) and for the original data (right).

  subcomposition of size d -1 with component j omitted.2. Define a new compositionp = (p k ) of size d = d -1, with, for 1 ≤ k ≤ d, k = j, p k = p k /(1 -p j ), so that p has component j removed.3. If d > 1, return to step 1, with p in lieu of p, and d instead of d.

Note that we allow CoDa points with some null components, contrary to Aitchison's geometry(J.[START_REF] Aitchison | The statistical analysis of compositional data[END_REF]), which is only defined on the interior of the simplex, so that log-coordinates (log p i ), and the subsequent log-ratio transforms (alr,clr, ilr), are well-defined.

For discontinuous conditional cdf, one must use the conditional probability integral transform (14) instead, see Ludger[START_REF] Rüschendorf | On the distributional transform, Sklar's theorem, and the empirical copula process[END_REF]

Note that[START_REF] Azadkia | A simple measure of conditional dependence[END_REF] introduce a more general regression dependence coefficient which allows for covariates and the assessment of conditional independence.

Acknowledgements

[author] acknowledges funding from the French National Research Agency (ANR) under the Investments for the Future (Investissements d'Avenir) program, grant ANR-17-EURE-0010.


Lucas Maisonnave
email: lucas.maisonnave@cea.fr

Cyril Moineau
email: cyril.moineau@cea.fr

Olivier Bichler
email: olivier.bichler@cea.fr

Fabrice Rastello
email: fabrice.rastello@inria.fr

Applying maximum entropy principle on quantized neural networks correlates with high accuracy

Keywords: compression, quantization, deep learning, entropy

Neural network quantization, sparsification or Neural Architecture Search (NAS) have shown great success in reducing model size and computational cost on many different ML tasks and architectures. Appropriately reducing model size without degrading task performance requires a suitable measure for quantifying the amount of important information entailed in the model parameters. Information theory provides such a tool, called entropy, which, from a probabilistic point of view, makes it possible to determine whether a model has effectively extracted important information from the data. Based on information theory and the maximum entropy principle, this paper investigates the influence of the clamping function on the distribution of weights in a quantized neural network. We show that the entropy weights can be increased by using a trainable parameter that evolves during training. We also identify a correlation between high entropy and high performance. Finally, by using a regularizer that enables the model to further increase its entropy we highlighted the importance to split training into a generalization phase and an information optimization phase.

I. INTRODUCTION

The increase in performance of deep neural networks has been driven by the availability of ever-larger databases and the use of hardware adapted to parallelizable calculations like GPUs. This exponential increase in the size of databases and in the number of operations required during training [1] has given rise to a number of problems: high energy consumption and computation times, and an increase in memory footprint.

Today, the trend in deep learning is to constantly increase the size of models and databases. This trend cannot continue indefinitely, and will eventually comes up against the limits of computer memory and computing capacity. In the past years, such challenges have led to the emergence of neural network compression, which aims to reduce the size of AIs without affecting their performance. A great deal of work has been carried out to reduce the size of neural networks, with 4 main categories emerging: sparsification, quantization, automatic architecture search (NAS) and matrix decomposition. The aim of all these methods is to favor a simple model for a given task, in line with the minimal description length (MDL [START_REF] Barron | The minimum description length principle in coding and modeling[END_REF]) principle. It is a very important principle for embedded AI considerations, as it ensures frugality and computational acceleration. The aim is to find the most compact and accurate representation of the data, while maintaining good performance.

Deep Neural Networks can also be studied from an information-theoretic point of view, using Shannon entropy to quantify the amount of important information in a neural network. However, the distribution of weights and activations of many DNNs is not uniform and does not maximize its entropy. There is therefore a difference between the number of bits allocated to encode weights and the amount of information actually present in the distribution. Maximizing entropy to close this gap has been studied and tested in numerous works on DNNs [START_REF] Liu | Self-supervised learning via maximum entropy coding[END_REF] [8] [START_REF] Messerschmitt | Quantizing for maximum output entropy[END_REF] [10] [START_REF] Pereyra | Regularizing neural networks by penalizing confident output distributions[END_REF] and seems particularly well suited to neural networks quantization.

In this paper we study neural network quantization methods using information theory to apply the maximum entropy principles to achieve the most compact data representation possible and make each bit of information as useful as possible. We show that we can increase weights distribution's entropy using a trainable parameter β in the clamping function. We also show strong correlation between high entropy and high performance especially for low-bit quantization.

II. RELATED WORKS AND FONDATIONS A. Quantization

Quantifying parameters of a DNN consists in converting network floating-point 32 bits weights values with a smaller number of bits like 8, 4, or 2 bits (Figure 1) [START_REF] Deng | Model compression and hardware acceleration for neural networks : A comprehensive survey[END_REF]. The main objective is to reduce memory and computing footprint of network while maintaining high performances and favoring faster and less costly integer-only operations.

Quantization methods are numerous [START_REF] Choi | Pact: Parameterized clipping activation for quantized neural networks[END_REF], [START_REF] Jin | Towards Efficient Training for Neural Network Quantization[END_REF], [START_REF] Zhou | Training Low Bitwidth Convolutional Neural Networks with Low Bitwidth Gradient[END_REF], and can be used after (post training quantization) [START_REF] Banner | Post training 4-bit quantization of convolutional networks for rapid-deployment[END_REF] or during training (quantization aware training) [START_REF] Steven | Learned step size quantization[END_REF] , be uniform [START_REF] Choi | Pact: Parameterized clipping activation for quantized neural networks[END_REF] or non-uniform [START_REF] Zhou | Incremental network quantization: Toward lossless CNNs with lowprecision weights[END_REF]. Other quantization schemes introduce randomness to choose the upper or lower bin to round values [START_REF] Chen | A statistical framework for low-bitwidth training of deep neural networks[END_REF]. Some methods allocate a different number of bits depending on the model layer [START_REF] Jin | Fracbits : Mixed precision quantization via fractional bitwidths[END_REF] (Mixed Precision).

In our experiments we will base our method on Scale Adjusted Training (SAT [START_REF] Jin | Towards Efficient Training for Neural Network Quantization[END_REF]) which quantize weights using DoReFa [START_REF] Zhou | Training Low Bitwidth Convolutional Neural Networks with Low Bitwidth Gradient[END_REF] and activations with PArameterized Clipping ac-Tivation (PACT [START_REF] Choi | Pact: Parameterized clipping activation for quantized neural networks[END_REF]). DoReFa introduces a clamping function that bounds weights between [0, 1] (eq. 1). This function uses tanh as a non-linearity to clip outliers before quantization (eq. 2) instead of PACT which uses a trainable parameter as scaling factor to overcome the non symmetric and unbounded activation distributions.

Wl,i = Tanh 0 (W l,i ) = 1 2 tanh(W l,i ) max r | tanh(W l,r )| + 1 (1) W q l,i = 2Q( Wl,i ) -1 = 2 (2 b -1) round((2 b -1) Wl,i ) -1 (2)
W l,i being the i th weight of the l th layer and b the bitwidth. SAT also applies constant variance scaling to quantized weights of linear layers without batch normalisation, n out being the number of weights in the layer :

W * l = W q l n out VAR(W q l ) B.

Information theory

Information theory introduces tools to measure information in a set of elements [START_REF] Cover | Elements of information theory[END_REF] in a probabilistic way by using the log 2 function of the probability distribution of a random variable X ∼ p(x) (eq. 3). This way we can add information of two independent random variables.

I(X = x) = -log 2 (p(x)) (3)
We can now define entropy as the amount of uncertainty in the entire set of discrete values by using the expectation of information across all values:

H(X) = E X∼p(x) (I(x)) = - x∈X p(x) log 2 (p(x)) (4)
Entropy is a fundamental measure of information and has been used in deep learning to quantify information of a dataset filtered by a neural network [START_REF] Nicholas Darlow | What Information Does a ResNet Compress?[END_REF]. Moreover, information theory has been used as an attempt to explain and understand neural network efficiency [START_REF] Schwartz-Ziv | Opening the black box of Deep Neural Networks via Information[END_REF]. In this work entropy and mutual information are used to build a theory linking compression and generalisation of DNNs. Entropy can also be applied to design a regularizer to penalize confident output and maximize uncertainty [START_REF] Pereyra | Regularizing Neural Networks by Penalizing Confident Output Distributions[END_REF].

This leads to a fundamental principle in information theory : maximum entropy principle (MEP) which states that the distribution that best fits the data is the one with the highest entropy. This principle has been used for quantization in order to better optimize the amount of information allocated for weights and activations [START_REF] Liu | Self-supervised learning via maximum entropy coding[END_REF], [START_REF] Li | Q-vit : Accurate and fully quantized low-bit vision transformer[END_REF], [START_REF] Liu | Nonuniform-to-Uniform Quantization: Towards Accurate Quantization via Generalized Straight-Through Estimation[END_REF]. It can be shown that for a Gaussian distribution, maximizing entropy is equivalent to minimizing quantization error [START_REF] Messerschmitt | Quantizing for maximum output entropy[END_REF], which is the L 2 norm between weights and quantized weights :

max Q H(Q(W )) ⇐⇒ min Q ||W -Q(W )|| 2 (5)
This equation links a common information measure to a well studied objective in quantized neural networks like in [START_REF] Rastegari | Xnor-net : Imagenet classification using binary convolutional neural networks[END_REF] where the quantization error is an objective to minimise to get the best representation possible in 1 bit. In the next sections this equivalence will enable us to develop a regularizer to increase entropy.

III. METHOD TANHβ

In the literature DoReFa is a broadly common weights quantization scheme and its clamping function in (eq. 1) is widely used but there is no clear explanation on why tanh clamping would be a more appropriate choice than a naive minmax quantizer (eq. 6).

MinMax(W

l ) = W l -min r (W l,r ) max r (W l,r ) -min r (W l,r ) (6)
This function is simplier but has the same expected effect : to bound weights values between [0, 1]. In fact when we apply MinMax on a MobilenetV1 quantized in 4 bits we can increase Top1 accuracy by 2.14 points and Top5 accuracy by 1.31 points (table II).

One way to explain this gap is to study the amount of information encoded by these different clamping functions. In order to measure the entropy of a neural network we introduce a metric we named H norm which represents the distance .

H norm = 1 N N i=1 H(W i ) b i (7)
N being the number of layers and b i the bitwidth of the i th layer. If H norm = 1 our model reached its maximum entropy which is the uniform distribution for every layer. The normalisation by b i allows to compare layers with different bitwidth especially the first and last layers quantized in 8 bits. H is computed with the empirical distribution of every bin of quantization, then we compute the average normalised entropy of the network. In the table II we can see that MinMax also reaches a higher entropy than Tanh (+0.24). The clamping function seems to have more effect on the distribution and the model performances than just bounding weights values.

Wl = Tanh 0 β l W l σ 2 (W l ) + ϵ (8)
To study the effect of the clamping function on the model we introduce Tanhβ which allows the model to modify the entropy of its weights using a trainable parameter β in the clamping function for each layer (eq. 8), σ 2 being the empirical variance of weights. Indeed, weights values are often very small due to weight decay (figure 1.A) so tanh is almost equivalent to the identity function and doesn't change weight's distribution and the entropy. So the variance rescaling is necessary to escape this linearity of tanh for small values.

Figure 1.B shows the effect of β on a reduced gaussian distribution (µ = 0). We can see that β affects the distribution's shape and tends to be uniform for β = 0.86. We can show that β can make the clamping function transition from Tanh 0 to MinMax when it varies from 0 to 1 (eq. 9, 10).

lim β→0 Tanh 0 (βW ) ∝ MinMax(W ) (9)
lim β→1 Tanh 0 (βW ) = Tanh 0 (W ) (10)
By introducing the parameter β we expect the model to use it during training to flatten its weight's distribution through gradient descent. To help the model to optimize β in this way we can use a regularizer based on (eq. 5) to minimise the quantization error and thus maximise entropy :

R H (W ) = 1 N w∈W (w -Q(w)) 2 (11)
L = L cls + λR H (12)
In (eq. 11) we need to compare FP32 weights with quantized weights which usually have values in [-1, 1] (figure 1.C) and can't be compared to non-quantized weights which values are smaller. So we need to "unclamp" quantized bins to rescale values in the set of FP32 weights (black lines of W l in figure 1.A) and apply the regularizer R H to increase entropy of W q . Then we add this regularizer in the loss function with λ = 10 5 (eq. 12) to make its value significant on Mobilenet models (≃ 10 -1 , tested on one batch of Imagenet). In the following sections, we will first test the method without regularization, then add it to verify its effect on entropy.

Finally we make the hypothesis that :

• Tanhβ clamping can effectively increase entropy of a quantized model. In this section we will validate hypothesis stated previously using SAT quantization scheme based on DoReFa and PACT on MobilenetV1, MobilenetV2 and ResNet34. We will compare all 3 clamping functions introduced previously : Tanh 0 , MinMax, Tanhβ on Imagenet and CIFAR100. Experimental settings In our experiments models are trained from scratch with no pretrained weights to study in detail the influence of clamping during training. We train our models for 150 epochs for Imagenet and 100 for CIFAR100 with a learning rate 0.05, a weight decay 4e -5 and a batch size 256 with SGD optimizer. We initialize weights with a uniform distribution and ∀l, β l = 0.01 so tanh(βx) ∼ 0 id(x) because we don't want to change entropy at the beginning and let the model choose where to converge by itself during training.

A. Datasets and Implementation Details

B. Effect of Tanhβ on entropy

First we study the influence of β on the model's entropy for different bitwidth on CIFAR100. Table I shows Top1 accuracy and entropy for different models and bitwidth and we can see that Tanhβ always has the highest entropy and is very close to the maximum (0.933 for MobilenetV2 3 bits). Whereas MinMax has a lower entropy and a non uniform distribution due to its linearity that can't flatten weights. MinMax is also quite unstable and doesn't always converge in 2bits. With these first results we show that β can increase entropy and the model seems to optimize it in this way without regularization on CIFAR100.

Now we analyze the difference of entropy for each layer between a quantization scheme with Tanh 0 and Tanhβ in order to visualize the effect of β in the network. Figure 2 shows Fig. 2: Entropy gain G of Tanhβ for each layer compared to Tanh 0 for MobilenetV1 in 4bits. On average the entropy gain is equal to G mean = 0.843 bits the difference between the entropy of Tanhβ and Tanh 0 for every layer and G mean represents the average entropy gain. It appears that for every layer except the first and last one Tanhβ has increased entropy with an averaged gain of 0.848 bits for MobilenetV1 and 0.324 bits for MobilenetV2.

C. Correlation between entropy and accuracy

Now we will study the link between high entropy and high accuracy of a model. First we note that most of the time Tanhβ on CIFAR100 (Table I) has the best performances and can outperform non quantized model at very low precision (70.4% on MobileNetV1 3 bits). These results show potential to this new clamping function and seems to confirm our hypothesis.

On ImageNet (Table III) this correlation holds for all models and bitwidths. The best model is the one with the highest entropy which is in line with the maximum entropy principle. This trend seems more important for very low bits (3 or 2 bits) of quantization as we can see a huge drop of performance for MinMax (-10% Top1 accuracy in 3 bits) which is also correlated with a drop of entropy (-0.45 on H norm ).

To quantify this correlation we compute the Pearson's correlation ρ X,Y = cov(X,Y ) σ X σ Y betweeen Top1 Acc and H norm for CIFAR100 and Imagenet and we get 0.7 and 0.58 respectively. Now we want to know if the value of entropy is sufficient to explain this correlation so we analyze the evolution of β and H norm during training

D. Evolution of β

The value of β is very important parameter as it defines the entropy of a layer and, the shape of its distribution. We first notice that β is optimized differently on convolutional and linear layers (figure 3). We see that β converge to a very high value for linear layers which represents a distribution almost binary with 2 long bins at -1 and 1. On the figure 2 we see this drop of entropy on the last layer that reaches -4.1 and -5.1 bits of entropy for MobilenetV1 and MobilenetV2 respectively. We interpret that result as over-parametrization of linear layers which need less information to be efficient. β conv is the average of β for all convolutional layers and β lin is for the last linear layer.

E. Regularization R H

As presented in (eq. 11) we can use a regularizer to increase entropy and help the model to optimize β. As we can see in figure 4, minimizing quantization error seems to be a good objective as it increases entropy for a gaussian distribution quantized in 4 bits. But it is not perfect and entropy can decrease for very low quantization error due to approximations made to compute H in 4 bits of quantization. We now test this regularizer on Imagenet with MobilenetV1 and MobilenetV2 quantized in 4 bits (figure 5). We validate that R H increases entropy compared to Tanhβ without regularizer during training with MBV1 + Tanhβ converging to 0.901 (0.867 without R H ) and MBV2+Tanhβ to 0.932 (0.908 without R H ). Moreover we see that R H increases entropy at the very beginning of training and seems to remain almost constant.

However Top1 accuracy decreased with this regularizer (around 2%) which means that only looking to the value of entropy is not enough to understand its effect. Maximizing entropy at the beginning could be a part of the issue as seen in [START_REF] Pereyra | Regularizing neural networks by penalizing confident output distributions[END_REF] they set a threshold to decide when to start maximizing entropy. It's also in line with the results from [START_REF] Schwartz-Ziv | Opening the black box of Deep Neural Networks via Information[END_REF] where they stated there are 2 phases during training :

• Generalization. A first short phase where the model tries to fit data and increase accuracy • Compression. A long phase where the model compresses its information We state that the model first need to fit labels before optimizing its inner information with an entropic regularizer.

To go further in this study we could try to set a threshold to the regularizer, use a pre-trained model which already converged and could be easier to compress with entropy constraint, use a scheduler to gradually add the regularizer in the loss. In this paper we studied the impact of clamping functions on the entropy of a quantized neural network. We introduced a parameter β in the classic clamping function tanh in order to study its effect on entropy and accuracy. We found that high entropy is correlated to high accuracy in most of the cases on CIFAR100 and Imagenet. On CIFAR100 the model naturally optimizes β to increase its own entropy but it is the not case on Imagenet and it needs to be regularized. We built a regularizer R H minimizing quantization error which is equivalent to maximizing entropy for a gaussian distribution and used it during training. We showed that combining Tanhβ with this regularizer allows us to control and increase the entropy of weights. Entropic regularization also showed that only focusing on the value of entropy is not sufficient to explain the correlation with accuracy and that we need to study entropy over time during training.

Fig. 1 :

 1 Fig. 1: Quantization pipeline with Tanhβ. (1) First we apply our clamping function to W l which flatten the distribution depending on β , (2) then we quantize Wl to get our quantized weights used as our knew weights . (3) In our regularization R H only, we need to "unclamp" our quantized weights to compare them with FP32 W l .

  Datasets. The experiments are carried out on the ILSVRC12 ImageNet classification dataset. The ImageNet dataset is made of 1.2 million training images, and 50k validation images with 1000 classes. CIFAR100 is a simpler dataset than Imagenet and made of 50k training images and 10k validation images for 100 classes. In our experiments, we use the data augmentation method : resize and crop, horizontal flip.

Fig. 3 :

 3 Fig. 3: Evolution of β during training for MobilenetV1 4 bits. β conv is the average of β for all convolutional layers and β lin is for the last linear layer.

Fig. 4 :

 4 Fig. 4: Entropy over Quantization error for a Gaussian distribution using different values of β. 4 bits being the maximum can reach

Fig. 5 :

 5 Fig. 5: Entropy over time for a MobilenetV1 and MobilenetV2 in 4 bits with and without R H

TABLE I :

 I CIFAR100, Top 1 accuracy and H norm for different models, bitwidth and clamping functions. NaN is for a model that didn't converge. Most of the time our method Tanhβ outperform other functions on accuracy and entropy

	Clamping function Top1 Acc (%) Top5 Acc (%)	Hnorm
	Tanh	68.93	88.38	0.66
	MinMax	71.07	89.69	0.9
	Gain	2.14	1.31	0.24

TABLE II :

 II Accuracy and entropy gain for a quantized Mo-bilenetV1 4bits on Imagenet for 150 epoch between our weight distribution's entropy and the maximum it can reach for each layer (eq. 7)

TABLE III :

 III ImageNet, Top 1 accuracy and H norm for different models, bitwidth and clamping functions. NaN is for a model that didn't converge. We notice a correlation between high entropy and high accuracy IV. EXPERIMENTS