Claude Manté
email: claude.mante@mio.osupytheas.fr

Computing the Rao's distance between negative binomial distributions. Application to Exploratory Data Analysis

Keywords: Riemannian manifold, Negative Binomial, geodesics, cut locus

published or not. The documents may come

Computing the Rao's distance between negative binomial distributions.

Its rst derivative will be denoted α. A geodesic curve γ connecting two points p and q of M will be denoted p q, and p s ⊕ s q will denote the broken geodesic [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF] connecting p to q with a stopover at s. We will also consider for any θ ∈ M the local norm V g (θ) associated with the metrics g on the tangent space T θ M :

∀ V ∈ T θ M, V g (θ) := V t .g(θ).V .

(

) 1
The length of a curve α traced on M will be denoted Λ (α). A parametric probability distribution L i will be identied with its coordinates with respect to some chosen parametrization; for instance, we will write L i ≡ φ i , µ i for some negative binomial distribution. In addition, R + * := ]0, +∞[, and M F will denote the Frobenius norm of the matrix M ; logical propositions will be combined by using the classical connectors ∨ (or) and ∧ (and).

Introduction

The statistical analysis of counts of living organisms brings information about the collective behavior of species (schooling, habitat preference, etc), possibly associated with their socio-biological characteristics (aggregation, growth rate, reproductive power, survival rate, etc). Such data consist in general of two-way r × c tables of counts, whose rows are associated with surveys (spatialtemporal positions, generally) and columns are associated with species. Roughly speaking, these tables can be analyzed through two dierent approaches. On the one hand, multivariate methods are widely used to investigate relationships between the community structure and the spatio-temporal variations of the surveys, frequently in connection with explanatory environmental variables (see for instance [START_REF] Manté | A functional data-analytic approach to the classication of species according to their spatial dispersion. Application to a marine macrobenthic community from the Bay of Morlaix (Western English Channel)[END_REF][START_REF] Kidé | Spatio-Temporal Dynamics of Exploited Groundsh Species Assemblages Faced to Environmental and Fishing Forcings: Insights from the Mauritanian Exclusive Economic Zone[END_REF] and the references therein). On the other hand, an alternative way, much earlier used in Ecological Statistics, consists in modeling the rows or the columns count distributions. In the spirit of [START_REF] Manté | A functional data-analytic approach to the classication of species according to their spatial dispersion. Application to a marine macrobenthic community from the Bay of Morlaix (Western English Channel)[END_REF], we propose an intermediary method, consisting in measuring the dissimilarity between species through the probability distribution of some characteristic, and analyzing the obtained dissimilarity table through MDS. In [START_REF] Manté | A functional data-analytic approach to the classication of species according to their spatial dispersion. Application to a marine macrobenthic community from the Bay of Morlaix (Western English Channel)[END_REF], this characteristic was the dispersion of each species while here it will be its abundance. There is a wide range of functional methods to deal with distributional data, fundamentally depending on the chosen metrics on the probabilities set. Recently, multivariate methods with a geometric dominance appeared in the literature, based on Riemannian structures equipping spaces of probability densities: non-parametric Fisher-Rao metrics [START_REF] Srivastava | Riemannian Analysis of Probability Density Functions with Applications in Vision[END_REF], Wasserstein metrics [START_REF] Seguy | Principal Geodesic Analysis for Probability Measures under the Optimal Transport Metric[END_REF]. But all these methods were designed in a non-parametric setting, for absolutely continuous distributions, while our data are discrete. In addition, ecological eld data are typically characterized by a large number of zeros (problematic for most of the methods above); that is why the negative binomial (NB) distribution is widely used to model catches of animals [START_REF] Bliss | Fitting the Negative Binomial distribution to biological data[END_REF][START_REF] O'neill | Use of binary and truncated negative binomial modelling in the analysis of recreational catch data[END_REF][START_REF] Vaudor | Comparing distribution models for small samples of overdispersed counts of freshwater sh[END_REF][START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF]. It is especially relevant for ecologists, because 1. it arises as a Gamma-Poisson mixture, whose parameters depend on the more or less aggregative behavior of the species, and on the eciency [START_REF] Fisher | The relation between the number of species and the number of individuals in a random sample of an animal population[END_REF][START_REF] Anscombe | Sampling theory of the negative binomial and logarithmic series distributions[END_REF] of the trap for catching it 2. it arises as the limit distribution of the Kendall [START_REF] Kendall | On some modes of population growth leading to R. A. Fisher's logarithmic series distribution[END_REF] birth-and-death model; in this setting, its parameters depend on the demography of the species (reproductive power, mortality, immigration rate) 3. one of its limit cases, the log-series distribution, is a natural model for collections (of animals, for instance) [START_REF] Williams | The Logarithmic Series and its application to biological problems[END_REF].

But while the parametric approach is quite sound from the ecological point of view (see [START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF] and the references therein), it is ill-suited for Exploratory Data Analysis (EDA): the visual distance between parameters of several distributions is misleading, because on the one hand it depends on the chosen parametrization and, on the other hand, because these parameters are not commensurable in general (dierent ecological meaning, dierent ranges, ...).

In a seminal paper, Rao [START_REF] Rao | Information and the Accuracy Attainable in the Estimation of Statistical Parameters[END_REF] noticed that, equipped with the Fisher informa-tion metrics denoted g (•), a family of probabilities depending on p parameters can be considered as a p-dimensional Riemannian manifold. The associated Riemannian (Rao's) distance between the distributions of parameters θ 1 and θ

2 is D R θ 1 , θ 2 := ˆ1 0 γt (t) .g (γ (t)) . γ (t)dt ( 2
)
where γ is a segment (minimal length curve) connecting

θ 1 = γ (0) to θ 2 = γ (1)
and γ (t) := dγ dt (t); as any Riemannian distance, D R is intrinsic. Naturally, Rao [START_REF] Rao | Information and the Accuracy Attainable in the Estimation of Statistical Parameters[END_REF][START_REF] Rao | Comment to Kass' paper[END_REF] proposed to use [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF] as a distance between populations or for Goodness-Of-Fit (GOF) testing, followed by a number of authors [START_REF] Carter | FINE: Fisher Information Nonparametric Embedding[END_REF][START_REF] Galanis | Wave height characteristics in the north Atlantic ocean: a new approach based on statistical and geometrical techniques[END_REF][START_REF] Dodson | Some illustrations of information geometry in biology and physics[END_REF][START_REF] Cubedo | A dissimilarity based on relevant population features[END_REF][START_REF] Ilea | Statistical Hypothesis Test for Maritime Pine Forest Sar Images Classication Based on the Geodesic Distance[END_REF][START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF][START_REF] Kass | The geometry of asymptotic inference[END_REF][START_REF] Menendez | Statistical Tests Based on Geodesic Distances[END_REF][START_REF] Cubedo | Hypothesis testing: a model selection approach[END_REF]. The Rao's distance between members of a common family of distributions has been calculated in a number of classical cases [START_REF] Burbea | Informative geometry of probality spaces[END_REF] but it cannot be obtained in a closed form, generally. In such cases, like the NB distributions (when both parameters are unknown), D R must be obtained by numerically solving a second-order nonlinear dierential equation which is frequently hard to integrate.

The outline of this study is as follows. After reminding in Section 2 essential notions of Riemannian geometry, we resume in Section 3.1 a method proposed by Manté and Kidé [START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF] for approximating D R θ 1 , θ Theorem 2. Let γ : I → M be a geodesic with respect to the metric connection ∇. Then γ has constant speed in the local norm (1)

γ g := γ (•) g (γ (•)) = γt (•) .g (γ (•)) . γ (•)
and, for any [a, b] ⊆ I, we have:

ˆb a γt (t) .g (γ (t)) . γ (t)dt = (b -a) γ g .
Geodesics on a p-dimensional Riemannian manifold with respect to ∇ are solutions of the Euler-Lagrange equation [START_REF] Gray | Modern dierential geometry of curves and surfaces with Mathematica[END_REF][START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF][START_REF] Burbea | Informative geometry of probality spaces[END_REF]:

∀ 1 ≤ k ≤ p, γk (t) + p i,j=1 Γ k i,j γi (t) γj (t) = 0 (3)
where each coecient (some Christoel symbol Γ k i,j ) only depends on g and is dened in coordinates by:

Γ k i,j := p m=1 g km 2 ∂g jm ∂θ i + ∂g im ∂θ j - ∂g ij ∂θ m (4)
where g im (resp. g mk ) is some entry of g -1 (resp. g). To determine the shortest path between two points of M, one applies the following result.

Theorem 3. [START_REF] Gray | Modern dierential geometry of curves and surfaces with Mathematica[END_REF][START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF] Let p, q ∈ (M, ∇, g) and suppose α : [a, b] → M is a curve of minimal length connecting p to q. Then, α is a geodesic.

Nevertheless, building the segment connecting L 1 to L 2 is not straightforward, since this theorem only says that a segment is a geodesic. But a geodesic is not necessarily a segment... Theorem 4. [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF][START_REF] Carmo | Riemannian geometry[END_REF] Let p = α (0) be the initial point of a geodesic. Then there is some 0 < t 0 ≤ +∞ such that α is a segment from p to α (t) for every t ≤ t 0 and for t > t 0 thereafter never again a segment from p to any α (t) for t > t 0 . This number t 0 is called the cut value of α and α (t 0 ) is called the cut point of α.

There are only two possible reasons (which can occur simultaneously) for α (t 0 ) to be to be the cut point of α:

• there is a segment from p to α (t 0 ) dierent from α

• α (t 0 ) is the rst conjugate point on α to p (i.e. t 0 α (0) is a critical point of the exponential map (see Denition 5 and Figures S1, S3 and S4).

In addition, the distance function D R (p, •) is not dierentiable at α (t 0 ) [START_REF] Itoh | Cut loci and distance functions[END_REF][START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF].

Remark 1. No matter the cause of the phenomenon, the main point for us is that if t 0 is a cut value of the unit-speed geodesic α,

∀ t ≤ t 0 , D R (p, α (t)) = t while ∀ t > t 0 , D R (α (0) , α (t)) < t. Nevertheless, ∀τ > 0, D R (α (t 0 ) , α (t 0 + τ )) = τ .
This remark is the basis of the method proposed by Manté and Kidé [START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF] for detecting cut points.

Denition 5. [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF] Let M be a Riemann manifold and x ∈ M. The exponential map of M at x is exp x : W x → M, dened on some neighborhood W x of the origin of T x M by: exp

x (V ) := α B(V ) ( V ) (5)
where B (V ) is the projection of V onto the unit ball and α B(V ) is the unique unit-speed geodesic in M such that α B(V ) (0) = x and αB(V

) (0) = B (V ). Remark 2. If α := p q is a segment and V 0 := α (0), because of uniqueness of geodesics, exp p (V 0 ) := α B(V0) (1) = q; reciprocally, if V 1 := -α (1), we have also that exp q (V 1 ) := α B(V1) (1) = p.

The special case of N B(D R )

There is a large number of parametrizations for the NB distribution, and the most classical one is probably

P (X = j; (φ, p)) =   φ + j -1 φ -1   p j (1 -p) φ j ≥ 0 (6)
with (φ, p) ∈ R + ×]0, 1[. Nevertheless, because of its orthogonality, we chose instead the parametrization used by Chua and Ong [START_REF] Chua | Test of misspecication with application to negative binomial distribution[END_REF]:

P (X = j; (φ, µ)) =   φ + j -1 j   µ µ + φ j 1 - µ µ + φ φ , j ≥ 0 (7) (φ, µ) ∈ R + × R + ;
here, µ is the mean of the distribution and φ is the so-called index parameter". In these coordinates, the information matrix is:

g(φ, µ) =   G φφ 0 0 G µµ  
where G µµ = φ µ(µ+φ) , while the expression of G φφ is more complicated:

G φφ = - µ + φ (µ + φ) ( φ /µ+φ) φ -1 ψ 1 (φ) φ (µ + φ) (8)
where ψ 1 is the Trigamma function [START_REF] Abramowicz | Handbook of mathematical functions with formulas, graphs and mathematical tables[END_REF].

One will nd in Burbea and Rao [START_REF] Burbea | Informative geometry of probality spaces[END_REF] the closed-form expression of the Rao's distance for a number of probability families. These authors reported that, when the index parameter of two NB distributions is the same, the Rao's distance is given by:

D N B(p) φ, p 1 , φ, p 2 := 2 φ cosh -1 1 -p 1 p 2 (1 -p 1 ) (1 -p 2 ) (9)
in the parametrization [START_REF] Seguy | Principal Geodesic Analysis for Probability Measures under the Optimal Transport Metric[END_REF]. Of course, if

L 1 = N B φ, p 1 (resp. L 2 = N B φ, p 2 ),
we have necessarily:

D R L 1 , L 2 ≤ D N B(p) L 1 , L 2 . ( 10
)
Due to the complexity of ( 8), D R L 1 , L 2 cannot be obtained in a closed-form.

It must be computed by nding the numerical solution of a the Euler-Lagrange equation ( 3), completed in the parametrization ( 7) by the conditions (boundary value problem)

γ (0) = φ 1 , µ 1 , γ (1) = φ 2 , µ 2 . ( 11
)
Geodesics can be as well be computed by solving (3) under the alternative constraints (initial value problem)

γ (0) = φ 1 , µ 1 , γ (0) = V ∈ R 2 ( 12
)
where V is the initial velocity of the geodesic; this solution is associated with the exponential map at φ 1 , µ 1 .

Numerical approximation of

D R L 1 , L 2 [1]
From now, L i ≡ φ i , µ i will denote some NB distribution parametrized in the ( 7) system, but our purpose could be extended to any parametric family of probabilities. Firstly, all the Christoel symbols (4) were calculated from the expression (8) of G φφ , with the help of Mathematica [START_REF] Inc | Mathematica[END_REF] . Then, the dierential equation ( 3) was numerically solved under the the boundary conditions [START_REF] Fisher | The relation between the number of species and the number of individuals in a random sample of an animal population[END_REF], for a number of distributions of counts of marine species whose parameters had been estimated in [START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF]. In most cases a solution could be found in an acceptable time (four CPU minutes), with a good numerical precision (20 digits), but was each one of the geodesics found a segment? And what about failures met in computation? We indeed had to face various problems detailed in [START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF], where some numerical remedies were proposed. The rst one consisted in inserting a well-placed stopover S between each pair of problematic distributions A and B, in such a way that D R (A, S) and D R (S, B) could be computed in a reasonable time, while D R (A, B) could not. Furthermore, S was placed in order that D R (A, S) + D R (S, B) should by a good approximation of D R (A, B). For sake of brevity, we moved to the supplementary material useful information and illustrations about this previous work. Notice that all references to this supplement will be preceded by an S.

Contribution of geometry: making computations easier thanks to Poisson approximation

From the numerical side, it is noteworthy that the index parameter φ often takes large values, causing diculties in the evaluation of quantities associated with Γ (φ), like formulas [START_REF] Bliss | Fitting the Negative Binomial distribution to biological data[END_REF] and [START_REF] O'neill | Use of binary and truncated negative binomial modelling in the analysis of recreational catch data[END_REF] or the Christoel's symbols [START_REF] Kidé | Spatio-Temporal Dynamics of Exploited Groundsh Species Assemblages Faced to Environmental and Fishing Forcings: Insights from the Mauritanian Exclusive Economic Zone[END_REF].

From the statistical side, the convergence of some NB distribution L ≡ (φ, µ)

towards a Poisson distribution P when φ → ∞ is well-known. More precisely, Majsnerowska [START_REF] Majsnerowska | A note on Poisson approximation by w-functions[END_REF] proved that

d T V (L, P (λ)) ≤ ∆ (φ, µ) := 1 -e -µ µ φ (13)
where λ := φ µ φ+µ and d T V denotes the total variation distance. Thus, we can claim that (φ µ)∨(µ small) ⇒ ∆ (φ, µ) small and conclude that in such cases it may be quite impossible to nd a dierence between L and P (λ), even when the index parameter is small or moderate! This fact suggests to replace the NB model by the Poisson one when both distributions are very close to each other. This is also biologically sound, since the former is well-suited for aggregative species, while the latter is associated to species with a random behavior (see [START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF][START_REF] Manté | A functional data-analytic approach to the classication of species according to their spatial dispersion. Application to a marine macrobenthic community from the Bay of Morlaix (Western English Channel)[END_REF] and the references therein).

Let's focus now on the application

     ω : R + * × R + * → R + * (φ, µ) → λ
associating to each NB distribution L ≡ (φ, µ) the corresponding limit Poisson distribution P(λ).

Lemma 6. The tangent application

T (φ,µ) ω = 1 (φ+µ) 2 µ 2 φ 2
is surjective.

Proof. Let us x some ρ ∈ R + * ; one can easily show that the set of solutions of the equation

T (φ,µ) ω (x, y) = ρ is the line of equation y = -µ φ 2 x+ρ 1 + µ φ 2
As a consequence, ω is a surjective submersion and the ber F λ := ω -1 (λ)

associated with any λ ∈ R + * is a sub-manifold of N B(D R ).
Proposition 7. F λ is dened by either equation:

µ (λ; φ) = λ 1-λ /φ : φ > λ φ (λ; µ) = λ 1-λ /µ : µ > λ . (14)
Proof.

∀ λ ∈ R + * , F λ := {(φ, µ) ∈ R + * × R + * : ω (φ, µ) = λ} = ∅ ∀ (φ, µ) ∈ R + * × R + * , (φ, µ) ∈ F ω(φ,µ) ∀ (λ 1 , λ 2 ) ∈ R + * × R + * , F λ1 ∩ F λ2 = ∅. Proposition 9. Suppose L ≡ (φ, µ) ∈ F λ and ∆ (φ, µ) ≤ δ,
where δ is some threshold chosen for deciding whether L can be identied with P (λ). Then, if L ≡ (φ , µ ) ∈ F λ is another distribution, such that φ > φ, ∆ (φ , µ ) < δ and L cannot be practically distinguished from P (λ) too. Proof. See Appendix 5.1

Let us now x L 0 ≡ φ 0 , µ 0 ; by denition L 0 ∈ F λ 0 , with λ 0 = ω φ 0 , µ 0 . Suppose ∆ φ 0 , µ 0 ≤ δ; then L 0 could be identied with P λ 0 , as well as any distribution of the ber whose index parameter is greater than φ 0 , due to the proposition above. We can now, thanks to [START_REF] Williams | The Logarithmic Series and its application to biological problems[END_REF], nd the distribution

(φ * , µ * ) λ 0 , δ such that φ * λ 0 , δ := arg φ:(φ,µ)∈F λ 0 (∆ (φ, µ) = δ) = arg φ (∆ (φ, µ (λ; φ)) = δ) (15)
and dene :

P λ 0 , δ := (φ, µ) ∈ F λ 0 : φ ≥ φ * λ 0 , δ .
Obviously, L 0 ∈ P λ 0 , δ because φ * ≤ φ 0 , but we have that, more generally, when (φ, µ) ∈ P λ 0 , δ , d T V L, P λ 0 ≤ δ and ω (φ, µ) -λ 0 ≈ 0, simultaneously. Thus, in such cases, N B (φ, µ) and P λ 0 are practically indiscernible from the statistical point of view! Denition 10. We will say that L ≡ (φ, µ) is Poisson-like if (φ, µ) ∈ P (ω (φ, µ) , δ).

We displayed on Figure 1 Poisson-like distributions: ≡ L 2 and L 3 δ ≡ L 2 ; there exists λ 1,2 : L i ∈ P (λ 1,2 , δ) , i = 1, 2 and λ 2,3 : Lemma 8), and these three Poisson-like distributions belong to the same ber. Thus, L 1 δ ≡ L 3 and δ ≡ is transitive Suppose L 1 δ ≡ L 2 belong to a common ber, F λ1,2 . Being indiscernible, these distributions should be necessarily close to each other, and it would be statistically sound to supersede D R L 1 , L 2 by δ.

L 1 δ ≡ L 2 ⇔ ∃λ : L i ∈ P (λ, δ) , i = 1, 2.
L i ∈ P (λ 2,3 , δ) , i = 2, 3. Consequently, L 2 ∈ F λ1,2 ∩ F λ2,3 which is empty if λ 1,2 = λ 2,3 (see
Corollary 12. Consider now two Poisson-like distributions L 1 and L 2 belonging to dierent bers, F λ1 and F λ2 . After computing (thanks to formula 15) φ := max φ * λ 1 , δ , φ * λ 2 , δ , μ1 := µ λ 1 ; φ and μ2 := µ λ 2 ; φ , we can determine the Poisson-like distributions L1 := φ, μ1 ∈ F λ1 and L2 := φ, μ2 ∈ F λ2 . Then, thanks to formula 9, we can easily compute D N B(p) L1 , L2 , which is an upper bound for D R L1 , L2 .

D R Ľ1 , L 2 . Since L1 δ ≡ L 1 δ ≡ Ľ1 and L2 δ ≡ L 2 , it is quite sound to supersede D R L 1 , L 2 by D N B(p) Ľ1 , L 2 or D N B(p)
L1 , L2 . One could also measure the dier- ence between L 1 and L 2 by D P P(λ 1 ), P(λ 2 ) but, since P(D P ) is not a sub-manifold of N B(D R ), there is no clear relationship between the associated Rao's distances (for instance, if L 1 = L 2 belong to the same ber F λ ,

D R L 1 , L 2 > D P ω L 1 , ω L 2 = 0). Nevertheless, one can easily prove that      D N B(p) L1 , L2 = C (μ 1 , μ2 ) D P P(λ 1 ), P(λ 2 ) D N B(p) Ľ1 , L 2 = C (μ 1 , µ 2 ) D P P(λ 1 ), P(λ 2 ) ( 16
)
where C (μ 1 , µ 2 ) and C (μ 1 , μ2 ) ≥ 1 are given by the function dened hereunder (up to a simple change of parametrization).

Lemma 14. C is dened in the parametrization (6) by Proof. In the parametrization ( 6), the mean µ = K p 1-p and thus p = µ φ+µ , while φ = K. Consequently, λ i = K p i and D P P(λ 1 ),

C p 1 , p 2 = cosh -1 1- √ p 1 p 2 √ (p 1 -1)(p 2 -1) -2 p 1 p 2 + p 1 + p 2 ≥ 1 if p 1 = p 2 ,
P(λ 2 ) = 2 √ K p 1 -p 2 .
Then, because of formula ( 9), we have: Suppose rst L 1 and L 2 belong to distinct bers F λ 1 and F λ 2 and each L i ∈ P λ i , δ . Then µ i ≈ λ i and we can use Corollaries 12 or 13 to build a pair of equivalent distributions, whose interdistance is easier to compute.

D N B(p) Ľ1 , Ľ2 D P (P(λ 1 ), P(λ 2 )) = C p 1 , p 2

An exemple

Suppose now L 1 is Poisson-like while L 2 is not; if φ 1 * ≤ φ 2 , we can again build Ľ1 := φ 2 , μ1 , and four distances can be computed: In the next case (Figure S3) L 1 = (0.00996246, 0.121282), while the second distribution is the same as in the previous case. The original distance corresponded to an intricate conguration with a cut point, and to the upper bound 43.1519. We found instead D R L 2 * , L 1 = 3.48809.

D R L 1 , L 2 , D R L 1 * , L 2 , D R Ľ1 , L
In the last case (Figure S4), L 1 = (0.938781, 9.86571), and we found that Ľ2 , L 1 , but we found that D R L 2 * , L 1 = 12.5294 (an intricate conguration with an acceptable rough solution and a stopover), while the original D R L 1 , L 2 , corresponding to an intricate cong- uration with linear interpolation, gave rise to the upper distance 21.351. we found that the proportion of simple congurations was greater than 90%, excepted for the second type of habitat, C2 (85%) (see Table 1). Furthermore, in the intricate cases, the rough solution was always accepted. In our previous study, numerical cut points were rare (less than a pair per class), but now we detect a single numerical cut point (see Table 1)! Poisson-like distributions were quite rare, but they were often so pathological that their replacement by equivalent NB distributions changed a lot the results. This is illustrated by the statistics displayed on Figure 4. Among the 325 distances computed, 305 were simple cases, with a median computation cost of 20"; the remaining 20 cases were intricate, with a median cost of 404". So, the introduction of Poisson approximation simultaneously increased the number of simple conguration, approximately divided by two the number of intricate congurations, and approximately divided by two the computational cost of the corresponding distances.

                     
Figure 4: Statistics of the computational cost for processing the same species as in Figure 3, with Poisson approximation. We display on Figure 5 the estimated parameters of the counts distribution of species sampled in the zone 4 of the MEEZ. This region is of paramount importance: it is a high plankton productivity area, supporting a large variety of sh communities, with many commercial species that sustain shing activities.

                                305
We distinguished two categories of species on Figure 5, according to the index parameter: φ > 1 or φ ≤ 1. Species belonging to the second category being very numerous, we only kept for MDS those which are situated on the convex envelope on the associated cloud. Twenty-six species were selected this way (while 138 species are represented) for computation of the Rao's distances and subsequent MDS of the table of distances. We also plotted on Figure 5 theoretical condence ellipses centered on the spatial medians [START_REF] Sering | Nonparametric multivariate descriptive measures based on spatial quantiles[END_REF] of both species categories.

Representation of the Rao's distance table

Only one of the 26 species retained in C4 was Poisson-like: HISP00 (Hippocampus sp). Consequently, the distance table was dened this way: ∆ i,j = D R L i , L j , excepted when one of these species was HISP00 (index h, say).

In these cases, we set instead: species is represented on Figure 6. Clearly, the counts of the species HISP00 and PHHU76 (Physiculus huloti, a type of cod) are distributed in a very special manner, which was not so obvious on Figure 5. In the inset region we represented the other 24 species separately analyzed by MDS of the restricted table ∆ AC (no point from the rst MDS is hidden). In this case, the median distortion associated with the Schoenberg transformation of the distances table ∆ Exp was 0.401832, greater than the one corresponding to ∆ AC (0.298623),

∆ h,j = min D R L h , L j , D R Lh , L j .
while the relative perturbation of distances was high c * /∆ 2 = 1.2529 .

Results and discussion

Following Rao [START_REF] Rao | Information and the Accuracy Attainable in the Estimation of Statistical Parameters[END_REF], a number of authors used D R (., .) in various statistical settings: either exploratory methods [START_REF] Rao | Comment to Kass' paper[END_REF][START_REF] Carter | FINE: Fisher Information Nonparametric Embedding[END_REF][START_REF] Galanis | Wave height characteristics in the north Atlantic ocean: a new approach based on statistical and geometrical techniques[END_REF][START_REF] Ilea | Statistical Hypothesis Test for Maritime Pine Forest Sar Images Classication Based on the Geodesic Distance[END_REF][START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF][START_REF] Cubedo | A dissimilarity based on relevant population features[END_REF][START_REF] Lebanon | Metric learning for text documents[END_REF] or hypothesis testing problems [START_REF] Menendez | Statistical Tests Based on Geodesic Distances[END_REF][START_REF] Cubedo | Hypothesis testing: a model selection approach[END_REF]. Motivated by the analysis of a large data set of marine species counts collected in the MEEZ, we developed a parameter-free method to compare species counts distributions in the setting of the Riemannian manifold N B(D R ) of negative binomial distributions, equipped with D R .

We focused rst [START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF] on numerical problems met in computing D R L 1 , L 2 : lengthy computations could result from the presence of a cut point on the geodesic L 1 L 2 , requiring to determine a stopover S somewhere between these distributions. D R L 1 , L 2 is then approximated by the upper bound D R L 1 , S + D R S, L 2 . In Section 3.2 we show how Poisson approximation can be used to evaluate more eciently D R L 1 , L 2 when one (at least) of the distributions involved is "Poisson-like". Superseding original NB distributions by equivalent ones (in the sense of Denition 10), we could obtain lower upper bounds of the distances than with the former strategy, with a lower computa-tional cost. In addition, this renement enabled us to get around most numerical issues (numerical cut points, unsuitable rough solutions).

Afterwards, an application to EDA of ecological data is shown. Unfortunately, we had to restricts ourselves to a small number of species, because of the computational cost of D R ; therefore, future research should focus on speeding up more the computation of D R L 1 , L 2 .

Appendices

Proof of Corollary 9

Notice rst that on the ber F λ , because of ( 14), the expression of ∆ dened in ( 13) is ≥ 0. Consequently, ∂∆ ∂φ (φ) is negative and ∆ (φ, µ) is a decreasing function of φ on a ber.

Pre-processing distance tables for MDS

There are several methods for making a distance matrix like ∆ Euclidean (i.e. nd a close distance matrix which can be exactly represented in an Euclidean space) -see for instance [START_REF] Benasseni | On a General Transformation Making a Dissimilarity Matrix Euclidean[END_REF]. The simpler one is the Additive Constant (AC) one [START_REF] Caillez | The analytic solution of the Additive Constant problem[END_REF], consisting in adding an optimal positive perturbation c * to all the extradiagonal terms of ∆ 2 . But other pre-processing methods are worth considering [START_REF] Benasseni | On a General Transformation Making a Dissimilarity Matrix Euclidean[END_REF]: one can search for the smallest positive γ 0 such that the power ∆ γ is Euclidean for γ ≤ γ 0 [START_REF] Joly | Etude des puissances d'une distance[END_REF], or the smallest positive γ * such that 1 -e -γ ∆ is Euclidean (Exp method). It is noteworthy that both these transformations belong to the class of Schoenberg transformations introduced in Data Analysis by Bavaud [START_REF] Bavaud | On the Schoenberg transformations in data analysis: Theory and illustrations[END_REF]. We chose the last one, since it is bounded and rectiable, i.e. a nite length curve is transformed into another nite length curve -fractal curves are not rectiable [START_REF] Tricot | Courbes et dimension fractale[END_REF], for instance. Thus, two element-wise transformations of ∆ were considered:

     ∆ AC i,j := ∆ 2 i,j + c * i = j ∆ Exp i,j := 1-exp(-λ * ∆i,j ) λ * i = j (17)
Naturally, these perturbation should be as small as possible. In the case of ∆ AC the ratio c * /∆ 2 , where ∆ 2 denotes the mean squared distance, is a straightforward and natural criterion. Remember now that since the set of Euclidean matrices of a given size is a convex cone, the solutions proposed in ( 17) can be easily improved by looking for an optimal η minimizing the Frobenius norm ∆ -η ∆ • F . In the applications, the original table ∆ • will be systematically superseded by this optimum. In addition to the specic index c * /∆ 2 for AC, it is interesting to consider statistics of local distortions.

Benasseni et al. [START_REF] Benasseni | On a General Transformation Making a Dissimilarity Matrix Euclidean[END_REF] proposed the criterion max i =j ∆ • i,j ∆ i,j /min i =j ∆ • i,j ∆ i,j which was illsuited for our data, since ∆ i,j was frequently very small. We instead computed for each distance obtained from (17) the index ω i,j :=

     0 if ∆ • i,j = ∆ i,j 1 - ∆i,j ∆ • i,j if ∆ • i,j = ∆ i,j
.

Then, each list Ω of distortions was described by its kernel density estimate, its average and its median.

  Poisson approximation, Multidimensional scaling Notations Consider a Riemannian manifold M, and a parametric curve α : [a, b] → M.

Lemma 8 .

 8 F λ := (φ, µ) : φ µ φ+µ = λ ; thus, the strictly positive parameters λ, φ and µ are linked by the relationship φµ = φλ+λµ, which proves that φ = λ+λ φ µ and µ = λ + λ µ φ . Consequently, λ < min (φ, µ) and lim φ→+∞ µ (λ; φ) = lim µ→+∞ φ (λ; µ) = λ One can easily verify that:

  four examples of such NB distributions (setting δ = 0.01, say). Let us now denote δ ≡ the following relation (for δ xed) between

Corollary 11 .

 11 The relation δ ≡ is an equivalence relation between Poisson-like distributions. Proof. Reexive and symmetric properties are straightforward. Suppose now L 1 δ

Figure 1 :Corollary 13 .

 113 Figure 1: Four instances of Poisson-like distributions; ∆ is given by Formula (13) and vertical bars are associated with NB probabilities while continuous curves are associated with Poisson ones

  and C (p, p) := 1.

Consider L 1 =≡ L 2 (

 12 (172.236, 0.974793) and L 2 = (6, 0.05), and lets us x δ := 0.01. Both these distributions are Poisson-like, with L 1 ∈ F 0.0495868 while L 2 ∈ F 0.969307 . We plotted on Figure 2 interesting portions of these bers. On each one of the panels, the big gray point (of coordinates (λ, λ)) corresponds to the lower bound of φ and µ, while L * := (φ * , µ * ) (λ, δ) is the distribution given by Equation 15. All the distributions situated on the right of L * are Poisson-like. It is the case of L 1 and L 2 , represented on Figure 2 by small gray points. We found that φ 1 * = 63.2096 < φ 1 and φ 2 * = 0.412537; thus, φ = φ 1 * and L1 = L 1 * , while L2 = L 2 * . Next, in accordance with Corollary 12, we com- pute D N B(p) L1 , L2 ≈ 1.53375, which is rather close to D P (P(λ 1 ), P(λ 2 )) ≈ 1.52371. Using Formula 16, one obtains 1.53375 again. Since φ 1 * = 63.2096 < φ 1 , it is also possible to determine the distribution Ľ2 := (φ 1 , μ2 ) δ black point on the lower panel) and, in accordance with Corollary 13, to compute D N B(p) Ľ2 , L 1 ≈ 1.52737 (very close to D P (P(λ 1 ), P(λ 2 )) too).

Figure 2 : 1 .

 21 Figure 2: Portions of bers associated with two Poisson-like distributions (δ = 0.01): F 0.0495868 (lower panel) and F 0.969307 (upper panel).

Ľ2 := φ 1 ,

 1 μ2 and computeD N B(p) Ľ2 , L 1 , but D R L 2 * , L 1 = 3.53253 could be computed straightforwardly (simple conguration, no cut point), while the original D R L 1 , L 2 corresponded to an intricate conguration with a cut point, and to the upper bound 45.1321 (denitions of simple and intricate congurations are remembered in the supplementary material).

L 2 =

 2 (172.236, 0.974793) δ ≡ (63.2096, 0.984403) = Ľ2 . Since φ 2 * = 63.2096 0.938781, we could not compute D N B(p)

Figure 3 :

 3 Figure 3: Statistics of the computational burden for processing (without Poisson approximation) species collected in zone 4 of the MEEZ: number of distances of each type, box-plots of computation length.

Figure 5 : 2 .

 52 Figure 5: Species collected in zone 4 of the MEEZ; the green dotted (resp. orange dashed) closed curve corresponds to the condence region of 0.99 level associated with the spatial median of the rst (resp. second) category species .

Figure 6 :

 6 Figure6: MDS of the Rao's distance between the 26 selected species (big points); the inset graph corresponds to the same analysis, performed after removing both the species HISP00 and PHHU76. .

3 .-λ -ζ e ζ ζ+λ 2 -λ- 1 + ζ + λ 2 (ζ + λ 2 ) 3 whose sign depends on the sign of 1 -e ζ ζ+λ 2 -λ- 1 .

 3212321 (λ 2 +λ-φ)-λ+φ (λ-φ) Since φ > λ, the denominator of this expression is always negative while the numerator is clearly positive, exceptedpotentially if λ 2 + λ -φ < 0. Substituting λ 2 + λ + ζ to φ (with ζ > 0) inthe equation, we get a simpler expression for ∂∆ ∂φ : Since the only solutions of ζ ζ+λ 2 -λ -1 = 0 are λ = 0 ∧ ζ = 0 and ζ = -λ + λ 2 ∧ λ = 0, this expression is negative, and 1 -e ζ ζ+λ 2 -λ-1

  2 and D N B(p) Ľ1 , L 2 . On the contrary, if φ 1 * > φ 2 , we cannot consider the two last ones.

	Three instances.

To illustrate our purpose, look rst at Figures S1, S3 and S4. In all these cases, only one of the pair of distributions is Poisson-like. In the rst case, L 1 = (0.00487399, 0.262591), and we found that L 2 = (592.392, 2.57454) δ ≡ (3.5634, 9.13442) = L 2 * . Since φ 2 * = 3.5634 0.00487399 we could not build

  The median distortion of the Schoenberg exponential transformation of ∆ (see Section 5.2) denoted ∆ Exp , was 0.854805, which is slightly lower than the one corresponding to the traditional Additive Constant method ∆ AC , 0.872421. Nevertheless, since ∆ -η AC ∆ AC Exp minimizes the Frobenius norm), ∆ AC has been chosen. Notice that the ratio c * /∆ 2 was 2.302: the additive constant method greatly altered the distance matrix. MDS of the 26

			F
	was much smaller than ∆ -η Exp ∆ Exp	F	(where η

Acknowledgments

We thank the Mauritanian Institute of Oceanographic Research and Fisheries (especially S. O. Kidé) and the Department of Cooperation and Cultural Action of the Embassy of France in Mauritania for their support for this study.

We also thank all scientists who contributed to eld surveys and data collection.

EDA of eld data: representation of counts distributions of marine species

The Mauritanian coast, situated on the Atlantic side of the northwestern African continent, embeds a wide long continental shelf of about 750 km and 36000km 2 , with an Exclusive Economic Zone (the MEEZ) of 230000km 2 . Manté et al. [START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF] considered the abundance of species of sh and invertebrates collected in the MEEZ during annual scientic trawl surveys since 1997 to now. Because the spatial distribution of groundsh species is strongly inuenced by the physical environment, we split this set into an optimal number (four) of subsets (typical habitats) associated with homogeneous physical conditions determined by available environmental variables (bathymetry, sedimentary type of the substrate, latitude and longitude). The counts associated with each species found in each one of the four habitats were then gathered, and tted by a truncated NB distribution; notice that only a reduced number of species could be satisfactorily tted in each habitat (for further information, see [START_REF] Manté | Fitting the truncated negative binomial distribution to count data. A comparison of estimators, with an application to groundshes from the Mauritanian Exclusive Economic Zone[END_REF]).

Benets of Poisson approximation

Processing these data, we found in [START_REF] Manté | Approximating the Rao's distance between negative binomial distributions. Application to counts of marine organisms[END_REF] an overwhelming proportion of simple congurations (more than 70%), while numerical cut points were quite rare. In the intricate cases, the rough solution was generally accepted (more than 90% of occurrences). We displayed on Figure 3


Mauro De Freitas

Malcolm Egan

Laurent Clavier

Anne Savard

Jean-Marie Gorce

Power Control in Parallel Symmetric α-Stable Noise Channels

Parallel channels form a basic building block for communication systems, including those based on OFDM and CDMA. While parallel Gaussian noise channels have been widely studied, parallel impulsive noise channels have received significantly less attention despite their importance in a range of modern communication systems. In this paper, this problem is addressed and a power allocation strategy is developed for parallel symmetric α-stable noise channels-a key class of impulsive noise channels. We show that our strategy can improve achievable rates by up to a factor of 1.5 over the standard waterfilling algorithm that assumes the noise is Gaussian.

I. INTRODUCTION

Impulsive noise arises in a range of communication systems and is often modeled via the α-stable distribution. In particular, the memoryless additive α-stable noise channel models interference in wireless communications for the Internet of Things (IoT) [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF]. In other wireless network scenarios, α-stable noise models have also been developed [START_REF] Pinto | Narrowband ommunication in a poisson field of ultrawideband interfers[END_REF]. Unlike Gaussian models, the α-stable distribution is characterized by heavy tails, which accounts for a high probability of large amplitude noise. However, due to the lack of closed-form expressions for the noise probability density function, characterizations of achievable rates in the presence of α-stable noise are challenging to obtain. A consequence is that the problem of resource allocation has not been addressed. Despite these challenges, there have recently been several new results characterizing the capacity of real additive α-stable noise channels subject to a range of constraints. In [START_REF] Fahs | A cauchy input achieves the capacity of a cauchy channel under a logarithmic constraint[END_REF], the capacity of the Cauchy noise channel (α = 1) was derived subject to a logarithmic constraint. More generally, capacity bounds for α > 1 subject to absolute moments constraints were obtained in [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF]. The extension to the complex isotropic α-stable channel was also studied in [START_REF] Egan | Achievable rates for additive isotropic α-stable noise channels[END_REF] and vector α-stable channels in [START_REF] Egan | On capacity sensitivity in additive vector symmetric alphastable noise channels[END_REF]. These capacity bounds now provide a basis to optimize systems that experience additive symmetric α-stable noise.

In this paper, we consider the problem of power allocation in parallel symmetric α-stable noise channels (1 < α < 2) subject to a power constraint. In the case of a Gaussian input in the presence of additive Gaussian noise, the rate-optimal solution is the well-known waterfilling algorithm. However, the waterfilling algorithm is tailored to the Gaussian noise channel. As such, it is highly desirable to develop alternative M. de power control strategies that do not rely on the Gaussian noise assumption.

We adopt a two-step approach to the design of power control for the symmetric α-stable noise channel. The first step is to select the input distribution. We note that the optimal input distribution for this channel is discrete and can only be obtained via numerical optimization. Therefore, we evaluate achievable rates via Monte Carlo estimates of the symmetric α-stable noise channel with Gaussian and truncated symmetric α-stable inputs. We show that Gaussian inputs perform comparably or outperform truncated symmetric α-stable inputs, and also nearly achieve a numerical approximation of the capacity obtained via the Blahut-Arimoto algorithm [START_REF] Blahut | Computation of channel capacity and rate-distortion functions[END_REF][START_REF] Arimoto | An algorithm for computing the capacity of arbitrary discrete memoryless channels[END_REF]. This is despite the fact that the truncated symmetric αstable inputs approximately match the input with the noise distribution, and are known to be a good choice with fractional moment constraints [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF].

The second step is to optimize the power control for the selected input distribution and develop a new power allocation scheme for Gaussian inputs. This is a challenging problem due to the fact that the achievable rate is not known in closedform for symmetric α-stable noise channels with Gaussian inputs. To this end, we introduce an approximation via a scale parameter matching method (detailed in Section IV-A) by exploiting the achievable rate with symmetric α-stable inputs obtained in [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF].

The resulting optimization problem is convex and therefore readily solved numerically, but differs from the waterfilling solution arising in Gaussian noise models. In particular, the KKT conditions lead to a system of non-linear equations unlike the linear equations arising from parallel Gaussian noise channels.

In summary, we make two key contributions: 1) We show that Gaussian inputs perform well on scalar additive symmetric α-stable noise channels subject to a power constraint, which is compatible with existing coding schemes. 2) We introduce a new power control scheme for parallel αstable noise channel with Gaussian inputs based on scale parameter matching (detailed in Section IV-A). This leads to a convex optimization problem with numerical results demonstrating that our power control schemes can outperform by up to a factor of 1.5 the rate achieved by waterfilling for Gaussian inputs, where the α-stable noise is assumed to be Gaussian. The remainder of this paper is organized as follows. In Section II, we detail the parallel symmetric α-stable noise model. In Section III, we study the effect of different input distributions. In Section IV, we develop our power control strategy. In Section V, we conclude the paper.

II. SYSTEM DESCRIPTION

Consider the real-valued memoryless additive symmetric αstable noise (ASαSN ) channel

Y = hX + N, (1)
where h ∈ R is a constant, X is the channel input, and N is symmetric α-stable noise with 1 < α < 2. We focus on the range 1 < α < 2 as in wireless communications this corresponds to a path loss exponent in the range (2, 4) [START_REF] Egan | Dynamic interference for uplink SCMA in large-scale wireless networks without coordination[END_REF], which captures a range of realistic electromagnetic environments. The α-stable random variables (α < 2) are an important class of random variables with heavy-tailed probability density functions and infinite second moments, which have been widely used to model impulsive signals [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF]. The probability density function of an α-stable random variable is parameterized by four parameters: the exponent 0 < α ≤ 2; the scale parameter σ ∈ R + ; the skew parameter β ∈ [-1, 1]; and the shift parameter δ ∈ R. As such, a common notation for an α-stable distributed random variable is X ∼ S α (σ, β, δ). In the case β = δ = 0, the random variable X is said to be symmetric.

In general, symmetric α-stable random variables do not have closed-form probability density functions. Instead, they are usually represented by their characteristic function, given by

E[e iθX ] = e -σ α |θ| α (2)
As a consequence of the lack of a closed-form probability density function, there are few closed-form characterizations of the capacity for symmetric α-stable noise channels. Nevertheless, it is possible to derive closed-form expressions for achievable rates. In particular, the rate of the ASαSN channel with a symmetric α-stable input was derived in [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF], given by

R = 1 α log 1 + |h| α σ α X σ α N , (3)
where σ X is the scale parameter of the symmetric α-stable input and σ N is the scale parameter of the symmetric α-stable noise.

The main problem we consider in this paper is power control for K parallel ASαSN channels with a sum power contraint. In this case, the system consists of K channels defined by

Y i = h i X i + N i , i = 1, 2, . . . , K, (4)
where h i ∈ R, X i is the real-valued input to the i-th channel and N i is real symmetric α-stable noise, independent for each i but not necessarily identically distributed. At present, there are no known closed-form expressions of achievable rates for power-constrained inputs in channels with symmetric αstable noise. Nevertheless, if the power constraint is relaxed, it follows from (3) that the sum-rate for parallel channels achieved using a symmetric α-stable input for each channel is given by

R S = K k=1 1 α log 1 + |h k | α σ α X,k σ α N,k . (5)

III. INPUT DISTRIBUTION SELECTION

Although the optimal input distribution for the power constrained additive Gaussian noise channel is well-known to be Gaussian, this is not the case for symmetric α-stable noise channels. As such, it is challenging to optimize the power control in the case of parallel channels and it is highly desirable to obtain input distributions that yield a high achievable rate with a simple parametric form.

We note that the optimal input distribution for an ASαSN channel is discrete [START_REF] Fahs | On properties of the support of capacityachieving distributions for additive noise channel models with input cost constraints[END_REF] and does not have a simple parametric form. While it is possible to obtain numerical approximations of this optimal input, even for a single channel it is time consuming [START_REF] Huang | Characterization and computation of optimal distributions for channel coding[END_REF]. Since the optimal input distribution depends on the power constraint, this means that it is also difficult to implement such an input within the context of a power control algorithm as both the input and power need to be jointly optimized. Using existing numerical optimization methods, this leads to very slow power control, which is not useful when, for example, fading is time-varying.

To overcome this difficulty, we investigate the choice of the input distribution for the ASαSN channel in (1) subject to a power constraint. We search for the distribution that allows an achievable rate close to the capacity obtained through the mutual information optimization problem:

maximize µ∈P I(X; Y ) subject to E µ [X 2 ] ≤ P, ( 6
)
where P is the set of probability measures on (R, B(R)).

In order to investigate the choice of the input distribution, we consider the following four choices:

(i) Zero-mean Gaussian input X G with probability density function

p X G (x) = 1 √ 2πσ 2 exp - x 2 σ 2 . (7)
(ii) Truncated symmetric α-stable inputs: let X S be a symmetric α-stable random variable, then the truncated symmetric α-stable input X T with truncation level T is constructed via

X T = X S , |X S | ≤ T sign(X S )T |X S | > T. ( 8
)
The power of the truncated symmetric α-stable input is given by

E[X 2 T ] = T -T x 2 p X S (x)dx + 2 ∞ T T 2 p X S (x)dx, (9)
where p X S is the probability density function of the symmetric α-stable random variable X S .

(iii) Truncated Gaussian inputs: let X G be a Gaussian random variable. The truncated Gaussian input X G,T can be constructed as the truncated α-stable input based on (8), replacing X S by X G . The power of this input is given by (9), replacing p X S (x) by p G (x), the probability density function of the Gaussian random variable X G . (iv) Numerical approximation of the capacity via the Blahut-Arimoto algorithm [START_REF] Blahut | Computation of channel capacity and rate-distortion functions[END_REF][START_REF] Arimoto | An algorithm for computing the capacity of arbitrary discrete memoryless channels[END_REF]. We selected the Gaussian input due to the fact that existing codes are typically constructed to approximate Gaussian inputs. The truncated symmetric α-stable input is chosen because it approximately matches the noise distribution (nearly optimal for fractional moment constraints [START_REF] Freitas | Capacity bounds for additive symmetric α-stable noise channels[END_REF]) and also satisfies the finite power constraint. It is also an appropriate choice of input for the case where the channel is both power and amplitude constrained. The truncated Gaussian input is selected as it forms a natural choice in the case of both power and amplitude constraints. The numerical approximation of the capacity forms a baseline to assess the optimality of the other three schemes.

To understand how the choice of input distribution affects the achievable rate, Fig. 1 plots the achievable rates using a Gaussian input, a truncated symmetric α-stable input, and also a truncated Gaussian input. In each case, the power is constrained to be P = 3 and the figure shows the impact on the truncation level for each input distribution. As achievable rates of additive symmetric α-stable noise channels are not known, in the experiments they are estimated via Monte Carlo simulation. In particular, we use 5 • 10 6 input samples, the entropy of the output and the noise are obtained by estimating the corresponding probability density functions via the kernel method [START_REF] Bowman | Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations[END_REF], which was performed by using a grid of 10 6 points and support [-200, 200 We observe in Fig. 1 that the Gaussian input outperforms both the truncated Gaussian and truncated symmetric α-stable inputs. Similarly, for most choices of the truncation level, the truncated Gaussian input also outperforms the truncated symmetric α-stable input. Moreover, the truncation level rapidly has no effect on the achievable rate for the truncated Gaussian input. Importantly, the Gaussian input also yields a rate that is very close to a numerical approximation of the capacity obtained via the Blahut-Arimoto algorithm [START_REF] Blahut | Computation of channel capacity and rate-distortion functions[END_REF][START_REF] Arimoto | An algorithm for computing the capacity of arbitrary discrete memoryless channels[END_REF].

We remark that based on extensive numerical experiments, we have observed that these trends hold for a wide range of channel parameters. This suggests that as in the Gaussian noise channel, a Gaussian input is a reasonable choice for the symmetric α-stable noise channel.

IV. POWER CONTROL STRATEGY

In this section, we develop a power control strategy for Gaussian inputs in parallel ASαSN channels. The lack of analytical expression for capacity or even achievable rates in the case of power constrained input makes the problem complex. Motivated by the results obtained in Section III, we propose to view the Gaussian inputs as approximations of symmetric α-stable inputs. This is possible since both of these inputs lie in the α-stable family and a scale parameter matching is used to approximate the Gaussian distribution via an α-stable distributions. This justifies the use of ( 5) to approximate the sum-rate with Gaussian inputs and a simple analytical form for the power optimisation problem. We verify the performance of our strategy via numerical simulation.

A. Scale Parameter Matching

Due to the fact that the input is Gaussian, no tractable closed-form expressions are known for the rate over α-stable noise channels. As a consequence it is highly desirable to approximate the Gaussian statistics by another α-stable random variable so that the sum-rate in (5) can be exploited. In order to do this, we propose a scale parameter matching method which selects the power of the Gaussian input by optimizing the scale parameter of another α-stable input.

Any zero-mean Gaussian random variable X ∼ N (0, Var(X)) is also a stable random variable S α (σ X , 0, 0), where the link between the scale parameter of the stable distribution σ X and the variance of the Gaussian notation is given via

σ 2 X = 1 2 Var(X). ( 10
)
The idea of scale parameter matching is to approximate the Gaussian input distribution (X ∼ N (0, 2σ 2 Z )) by a symmetric α-stable random variable (Z ∼ S α (σ Z , 0, 0)).

For our power control problem, this scale parameter matching method allows for the approximation of the sum-rate with Gaussian inputs by the sum-rate with symmetric αstable inputs. As the rate of symmetric α-stable inputs for the additive symmetric α-stable noise channel is known, a tractable power control problem can be formulated as detailed in Section IV-B.

To justify the scale parameter matching approximation, it is necessary to establish how close the rate with Gaussian inputs is to the rate with symmetric α-stable inputs. To this end, we have the following bound.

Theorem 1. Let X G ∼ N (0, 2σ 2 X ), X S ∼ S α (σ X , 0, 0), Y G = X G + N and Y S = X S + N . Then, |I(X G ; Y G ) -I(X S ; Y S )| ≤ max h min - 1 α log 1 + √ P √ 2σ N α , log(2eb) log 2 -h(N ) - 1 α log 1 + √ P √ 2σ N α =: E B (P ), (11)
where b = 2P π + 2 π Γ 1 -1 α σ N , h min = max{h(X G ) - h(N ), 0} and h(X G ) = 1
2 log(2πeP ). Proof. Omitted due to space constraints.

An observation from [START_REF] Bowman | Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations[END_REF] is that as P → ∞, E B (P ) tends to a constant. This provides evidence that the scale matching approach gives consistent results.

In order to gain further insights into the behavior of the approximation, Fig. 2 plots both the actual relative error and the relative errors using the bound in [START_REF] Bowman | Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations[END_REF] for varying power levels P . A key observation is that the relative error in both cases is approximately constant. This implies that the bound and actual relative error have the same qualitative behavior. Moreover, despite the fact that the relative error is greater than 10%, we note that for a wide range of values of P

|I(X

G ; Y G ) -I(X S ; Y S )| 1 α log 1 + √ P √ 2σ N α ≈ κ, ( 12
)
where κ is a non-negative constant. It then follows that

I(X G ; Y G ) ≈ (κ + 1) 1 α log 1 + √ P √ 2σ N α . (13)
For the purpose of optimizing P , the factor of (1 + κ) does not affect the solution. This provides a justification of the scale parameter matching method, which is further validated by numerical results in Section IV.

B. Optimization Problem Formulation

The scale parameter matching method provides a means of approximating the rate of a symmetric α-stable noise channel with a Gaussian input by a symmetric α-stable noise channel with an α-stable input. This method yields a power control optimization problem given by maximize σ X,k , k=1,2,...,n n k=1

1 α log 1 + |h k | α σ α X,k σ α N,k (14)
subject to n k=1 2σ 2 X,k ≤ P max σ X,k ≥ 0, k = 1, 2, . . . , n. In particular, the parameter σ X,k is the parameter for a symmetric α-stable input. Using the scale parameter matching method, the Gaussian inputs are assumed to have the same parameters σ X,k and as such, the power levels of the inputs are obtained via [START_REF] Huang | Characterization and computation of optimal distributions for channel coding[END_REF]. Note that this relationship is consistent with the fact that the constraint in ( 14) is in fact a sum power constraint for Gaussian inputs.

To solve (14), we apply the transformation ρ k = σ 2 X,k , which yields the problem maximize ρ k , k=1,2,...,n n k=1 We observe that the problem in (15) is convex, which follows from the fact that the function ρ α/2 k is concave for 0 < α < 2 and the linearity of the constraints or from the computation of the Hessian for the objective in (14).

1 α log 1 + |h k | α ρ α/2 k σ α N,k (15)
Observe that unlike the Gaussian noise model, the system of equations arising from the KKT conditions cannot be reduced to simple function parameterized by λ, able to be found using the bisection method. That is, either the non-linear system of equations must be solved numerically or a general purpose convex optimization solver must be applied to (14).

C. Numerical Results

We compare the proposed approach with the waterfilling algorithm designed for Gaussian noise in the case of two parallel channels. In applying the waterfilling algorithm, we assume that the system does not know the noise is non-Gaussian. As such, the variance of the noise is estimated by observing N S = 5 • 10 6 samples and applying the estimator

σ2 G,k = 1 N S -1 N S i=1 n 2 i,k , k = 1, 2, (16)
where n i,k is the i-th noise sample on the k-th channel. Note that since the variance of α-stable noise is infinite, it follows that the variance estimate in (16) does not converge. Nevertheless, (16) provides a means of systematically choosing the noise variance parameter required for the waterfilling algorithm, corresponding to the behavior of a system that does not know the noise is non-Gaussian.

In order to provide a fair comparison with power allocation based on our proposed method, the exponent α and the dispersion are also estimated based on N S = 5 • 10 6 samples, using the characteristic function method [START_REF] Mcculloch | Simple consistent estimators of stable distribution parameters[END_REF]. This is to ensure that noise parameters are estimated rather than assumed known. In the experiments, the scale parameter of the symmetric stable noise is σ N,k = 0.1, k = 1, 2 and 5•10 6 Gaussian input samples are generated. Fig. 3 shows the estimated achievable rate for each choice of α, channel h and power allocation method. The rates are estimated for several choices of the channel vector h using the same procedure as for Fig. 1 with 50 Monte Carlo iterations. Observe that our proposed strategy implemented in CVX [START_REF] Grant | CVX: Matlab software for disciplined convex programming, version 2.0 beta[END_REF] outperforms the waterfilling algorithm for each choice of parameters. In particular, for α = 1.4 and h 1 ∈ {0.1, 0.5} and h 2 varying between 0.1 and 0.9 an increase by up to a factor of 1.5 is achieved.

Fig. 4 plots the achievable rates for varying α and a fixed channel h = [0.9 0.7]. Observe that the proposed strategy outperforms waterfilling in all cases. However for α ≈ 2, there is negligible difference between the two approaches and our proposal closely matches the optimal waterfilling.

V. CONCLUSION

We have considered the problem of power control for parallel symmetric α-stable noise channels. We have shown that Gaussian inputs are a good choice, consistent with the Gaussian noise case. We then developed a new power control strategy for Gaussian inputs tailored to symmetric α-stable noise. This strategy significantly outperforms the rate achieved when the impulsive nature of the noise is ignored. A natural extension of our model is to MIMO systems, which is a target for future work.

  Freitas, L. Clavier and A. Savard are with Université de Lille 1, IEMN UMR CNRS 8520, IRCICA USR CNRS 3380, Villeneuve d'Ascq, France. M. Egan and J.-M. Gorce are with Univ. Lyon, INSA Lyon, INRIA, CITI, France. L. Clavier and A. Savard are also with the Institut Mines-Telecom, TELECOM Lille 1, France.

Fig. 1 .

 1 Fig. 1. Comparison of achievable rates using a truncated symmetric α-stable input (α = 1.4, E[X 2T ] = 3), a Gaussian input and a truncated GaussianInput (E[X 2 G ] = E[X 2 G,T ] = 3) in the presence of symmetric α-stable noise (α = 1.4, σ N = 0.1).

Fig. 2 .

 2 Fig. 2. Relative error of bound for α = 1.5 and σ N = 0.1.

  P max ρ k ≥ 0, k = 1, 2, . . . , n.

Fig. 4 .

 4 Fig. 4. Achievable rates varying α and h = [0.9 0.7]

  ].

		2.2
		2
		1.8
		1.6
	I(X,Y) (nats)	1.4					Gaussian Input Truncated -Stable Input Truncated Gaussian Input Blahut-Arimoto Capacity Approximation
		1.2
		1
		0.8
		0.6	0	2	4	6	8	10	12	14
						Truncation level T

  Fig. 3. Throughput gain, where h= [h 1 h 2 ] for h 1 ∈ {0.1, 0.5}

		1.6
		1.5
	(nats)	1.4							h 1 h 1	= 0.5 = 0.1
	waterfilling	1.3
	/I
	proposed	1.2
	I
		1.1
		1
		0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
						h 2
		3.8
		3.6
		3.4
	Achievable rate (nats)	2.6 2.8 3 3.2						proposed waterfilling
		2.4
		2.2
		2
		1.4		1.5	1.6		1.7	1.8		1.9

ACKNOWLEDGEMENTS

This work has been (partly) funded by the French National Agency for Research (ANR) under grant ANR-16-CE25-0001 -ARBURST.