Gilles R Ducharme
email: gilles.ducharme@umontpellier.frlafaye@unsw.edu.au

Pierre Lafaye De Micheaux

A GOODNESS-OF-FIT TEST FOR ELLIPTICAL DISTRIBUTIONS WITH DIAGNOSTIC CAPABILITIES

Keywords: Primary 62F03, 62H05; secondary 62E10 Diagnostic Information, Elliptical Distribution, Goodness-offit test, Multivariate Laplace distribution, Multivariate logistic distribution, Multivariate normal distribution, Multivariate Pearson type II distribution, Smooth tests

This paper develops a smooth test of goodness-of-fit for elliptical distributions. The test is adaptively omnibus, invariant to affinelinear transformations and has a convenient expression that can be broken into components. These components have diagnostic capabilities and can be used to identify specific departures. This helps in correcting the null model when the test rejects. As an example, the results are applied to the multivariate normal distribution for which the R package ECGofTestDx is available. It is shown that the proposed test strategy encompasses and generalizes a number of existing approaches. Some other cases are studied, such as the bivariate Laplace, logistic and Pearson type II distribution. A simulation experiment shows the usefulness of the diagnostic tools.

1. Introduction. Elliptically contoured (EC), or elliptical for short, distributions have become important tools in the analysis of multivariate data. They retain from the multivariate normal (MVN) distribution the feature of elliptical symmetry about a location µ. They extend the MVN to allow modelling data with short or large tails. Tools of multivariate analysis such as regression, correlation, PCA, discrimination, are easily ported to them. Many inferential procedures from MVN theory remain valid, after slight modifications. They pave the way toward more flexible models, such as skew-elliptic distributions, or more specialized, such as elliptical copulas, elliptical complex and elliptical matrix distributions. Consequently, they are used in many applications, such has portfolio theory, risk management, radioimmunoassay, astronomy, physics, signal processing, etc. (see [START_REF] Chmielewski | Elliptically Symmetric Distributions: A Review and Bibliography[END_REF] for a bibliography).

In this paper we focus on m-dimensional EC distributions with density of the form (1.1) f (x; η) = c m det(V -1 ) φ m ((x -µ) T V -1 (x -µ)),

where c m = Γ(m/2)/ π m/2 ´∞ 0 y m/2-1 φ m (y) dy is a normalizing constant, the density generator φ m (•) is defined over R ≥0 and the parameter η = µ, V -1 ∈ Ξ = R m × M m + , the space of m × m positive definite matrices, is unknown. To give a few examples taken from [START_REF] Bilodeau | Theory of multivariate statistics[END_REF]Chap. 13], the MVN has φ m (y) = exp{-y/2}; the multiStudent with ν degrees of freedom has φ m (y) = (1 -y/ν) -(m+ν)/2 ; the multivariate α-th power exponential distribution has φ m (y) = exp {-y α }. The multivariate slash-elliptical distribution with shape parameter α [START_REF] Wang | The multivariate skew-slash distribution[END_REF] has φ m (y) = y -(m+α)/2 ´y/2 0 t (m+α)/2 e -t dt. A more elaborate version, the multivariate extended slash-elliptical distribution has been introduced recently [START_REF] Rojas | An extension of the slashelliptical distribution[END_REF] to fit heavy-tailed data. Additions to the catalogue of EC distributions are periodically made in response to challenging new data sets.

To exploit these statistical models in applications, a user needs tools to help in selecting a distribution from this growing catalogue. One such tool is a Goodness-of-Fit (GoF) test of the null hypothesis that some data arise from a postulated EC density of the form (1.1). The goal of the present work is to develop such a GoF test having attractive features.

An important feature of a GoF test is its power in detecting departures from the null density. GoF tests can be roughly divided into omnibus and directional. An omnibus test will have power converging to one (i.e. be consistent) for any departure from the null. A directional GoF test will be consistent for only some departures. At first glance, the omnibus property appears essential. But equal power cannot be allocated to each departure and, for finite samples, only some will have significant power (see [START_REF] Shorack | Empirical Processes with Applications to Statistics[END_REF]Chap. 5]). Thus the interest for compromising approaches offering some control over the power allocation, such as the smooth test paradigm introduced by [START_REF] Neyman | Smooth tests for goodness-of-fit[END_REF]. The smooth test is directional along a set of K departures, but the value of K can be selected. This creates a slider between directional and omnibus GoF tests that yields an adaptive form of the omnibus property which has been found, in many simulations, to be very effective. This has brought [44, p. 9] to recommend : "Don't use those other methods-use a smooth test!". In its original form, its main weakness is the lack of criteria to position K on the slider. To correct this, [START_REF] Ledwina | Data-Driven Version of Neyman's Smooth Test of Fit[END_REF] has introduced a version that selects K in a data-driven fashion. The resulting data-driven smooth test constitutes, power-wise, a significant improvement. This has prompted [START_REF] Kallenberg | Data driven smooth tests when the hypothesis is composite[END_REF] to sharpen [START_REF] Rayner | Smooth Tests of Goodness of Fit: An Overview[END_REF]'s recommendation into : "use a data-driven smooth test!".

Another desirable feature of a GoF test is the ability to extract, when the test rejects, some information regarding the aspects of the null density contradicted by the data. As stated in [15, p. 800] regarding the process of model checking: "the purpose of the checking is not simply to "accept" or "reject", but rather to discover what aspects of the data are not being fit well by the model ". A GoF test that can provide such information is said to have diagnostic (Dx) capabilities. The smooth test has Dx capabilities, see [START_REF] Ducharme | Tree based Diagnostic Procedures Following a Smooth Test of goodness-of-Fit[END_REF].

Encouraged by these features, much work ([21, 49, 29, 23, 43, 48] among others) has been done to extend the smooth test paradigm. However, for multivariate densities, few extensions have appeared. One hurdle comes from the fact that EC distributions are closed under the group of affine-linear transformations of the data. The invariance principle states that if a statistical problem is left invariant after a transformation, then its solution should also be invariant under this transformation, otherwise interpretability of the inference can be compromised. [19, p. 469] nicely summarizes the importance of this principle : "... any proposal to use a non-invariant test... must come together with a special justification". For the case of the bivariate normal, [START_REF] Bogdan | Data driven smooth tests for bivariate normality[END_REF] has developed a data-driven smooth test that cleverly combines the smooth test statistics for univariate normality. But it is not invariant to rotations, so that an unscrupulous user could simply rotate the data to reach a pre-specified conclusion. Also her approach breaks down for general EC distributions because a multivariate EC density does not always have the same density generator as its marginals (the so-called inconsistency property, see [4, p. 209]. [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF] develops a smooth GoF test invariant to rotations, but their method is confined to data lying on an hypersphere.

Because of this scarcity, users facing the problem of assessing the fit of a multivariate EC density may feel abandoned by statistical methodology. When the null density is the MVN, numerous GoF tests exist (see [START_REF] Koziol | Assessing multivariate normality: A compendium[END_REF][START_REF] Romeu | A comparative study of goodness-of-fit tests for multivariate normality[END_REF][START_REF] Henze | Invariant tests for multivariate normality : a critical review[END_REF][START_REF] Ozlem | Comparison of some multivariatre normality tests: A simulation study[END_REF] and references therein). However, GoF tests tailored to other EC distributions are almost non-existent (with the noticeable exceptions of [START_REF] Fragiadakis | Goodness-of-fit tests for multivariate Laplace distributions[END_REF] and [START_REF] Ebner | Asymptotic theory for the test for multivariate normality by Cox and Small[END_REF]) and a user must resort to a general-purpose GoF test, such as the multivariate version of the Cramer von Mises approach (see [START_REF] Mcassey | An empirical goodness-of-fit test for multivariate distributions[END_REF] for a short review), or ad hoc methods.

In this paper, we develop a smooth GoF test for EC densities of the form (1.1) with unknown parameter η, that is adaptively omnibus, has Dx capabilities and is affine-linear invariant. As a bonus, the asymptotic reference distribution is a standard χ 2 and its power function can be approximated by a sum of independent non-central χ 2 . The approach is based on the ideas in [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF] but adapted to the case of affine-linear transformations. Section 2 adapts the smooth test paradigm to the case of a general EC density and sets the stage for the several levels of invariance required in our approach: problem invariance, test statistic invariance and Dx invariance. Section 3 makes explicit the orthonormal basis on which our smooth test is based and derives the associated test strategy by exploiting a variant of Rao's score test. Section 4 shows that the smooth test statistic possesses a convenient explicit expression from which emerges a decomposition into a sum of three invariant and asymptotically independent χ 2 distributed components, which we refer to as the Q = UIR decomposition. It is explained that each component can provide interpretable Dx information about the aspects of the null EC density contradicted by the data. It is also explained how the omnibus/directional slider must be set up to preserve invariance in all components of this decomposition. Section 5 computes the Q = UIR decomposition for the multivariate normal (MVN) distribution and discusses how our smooth test extends and relates to a number of proposals for this problem. The R package ECGofTestDx can perform the necessary calculations. The case of the bivariate Laplace is also considered. Section 6 reports on an experiment that was conducted to see how the elements of the Q = UIR decomposition behave under a number of alternative densities somewhat representative of what can be encountered in actual situations. It is seen that the decomposition can indeed help in providing illuminating Dx information. Some points requiring further research are collected in Section 7. Proofs are confined to Appendix A and Appendix B gives the details for the computation of the smooth test for the MVN when m = 2, 3 as well as MATHEMATICA commands to extend in higher dimensions. Appendix C gives the smooth test for two other EC distribution, the bivariate logistic and Pearson type II distributions.

2. The Smooth Test Paradigm for Elliptical Distributions. Let X 1 , . . . , X n be independent and identically distributed m-dimensional observations with density ψ(•). Consider the problem of testing :

(2.1)

H 0 : ψ(•) ∈ F 0 = {f (•; η), η ∈ Ξ} vs H 1 : ψ(•) / ∈ F 0 ,
where f (•; η) has the form (1.1) with a given φ m (•). Here we consider the usual case where η = (µ, V -1 ) is unknown. We set the following :

Assumption A : The support of f (•; η) does not depend on η and the mapping from Ξ to F 0 is one-to-one. Moreover f (•; η) has a moment generating function so that all moments of X 1 exist. Also P[X 1 = 0] = 0.

Assumption B : η is an affine-equivariant estimator (in the sense of Definition 13.1 in [START_REF] Bilodeau | Theory of multivariate statistics[END_REF]) of η such that

√ n(η -η) = 1 √ n n i=1 (X i , η) + o p (1)
, where E 0 ( (X 1 , η)) = 0 and the covariance matrix V 0 ( (X 1 , η)) is finite. Now confine ψ(•) to the embedding family of densities

G = {g(•; η, h(•)) |η ∈ Ξ, h(•) ∈ H and g(•; η, 0(•)) = f (•; η) } (2.2)
where 0(•) is the zero function assumed to belong to a suitable space of functions H. This allows to restate problem (2.1) as that of testing

(2.3) H 0 : h(•) = 0(•) vs H 1 : h(•) = 0(•).
Here we adopt from [START_REF] Ducharme | Goodness-of-fit tests for the inverse Gaussian and related distributions[END_REF] embedding functions of the form

(2.4) g(•; η, h(•)) = f (•; η) × (1 + h(•) -E 0 (h(X))) 2 1 + h(•) -E 0 (h(X)) 2 (f,η)
, where the subscript "0" refers to a statistical operator evaluated under H 0 and h(•) 2 (f,η) = ´h2 (x)f (x; η)dx. These are related to Hellinger's metric. The confinement to G, and thus the choice of H, affects the properties of our test. Regarding invariance, let AL(m) be the group of affine-linear transformations on R m with generic element γ = γ(x) = A(x + b), represented by the pair (A, b) where A is a m × m non-singular real matrix and b ∈ R m . EC distributions are closed under these transformations: the density of γ(X) has the form (1.1) with its parameters transformed from η = (µ, V -1 ) to γ * (η) = (A(µ + b), (A -1 ) T V -1 A -1 ) [4, p. 207]. Our GoF test, and its Dx information, must be invariant to such transformations. To this end, let

L 2 (f, η) = {h(•) | h(•) 2 (f,η) < ∞} and set H = {h ∈ L 2 (f, η) | h • γ(•) (f,η) < ∞, ∀γ ∈ AL(m)}.
With this choice, embedding family (2.2) does not depend on the value of η and is obviously closed under affine-linear transformations.This will be important to ensure invariance. Note that H is a slight restriction of the Banach space L 2 (f, η), so that almost any reasonable ψ(•) can be written as (2.4). This is important for the omnibus property.

We now recall the following terminology. A subspace

H of H is said AL(m)-invariant if h(•) ∈ H implies h • γ(•) ∈ H for all γ ∈ AL(m).
An AL(m)-invariant subspace is further said to be irreducible if it contains no nontrivial AL(m)-invariant subspaces. For any η ∈ Ξ, let •, • (f,η) be the scalar product associated with • (f,η) . H equipped with this scalar product is a Hilbert space denoted H (f,η) . With these definitions at hand, a variant of the argument in [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF] shows that H (f,η) can be decomposed into a sequence of disjoint irreducible AL(m)-invariant subspaces {Π k , k ≥ 0}, each of a finite dimension d m (k) that will be made explicit in Theorem 3.1 below. Hence, (2.5)

H (f,η) = ∞ k=0 Π k ,
where ⊕ is the direct sum operator with respect to •, • (f,η) and Π 0 is the set of constant functions. For each k ≥ 0, let {π k,j, (•; η), (j, ) ∈ B k } be a complete orthonormal basis (CONB) for Π k . Here B k is a set of d m (k) pairs of integers that will be explicited, along with a choice π k,j, (•; η) in Theorem 3.1. Thus Span{π k,j, (•; η), (j, ) ∈ B k } = Π k and for any (j, ), (j , ) ∈ B k , π k,j, (•; η), π k,j , (•; η) (f,η) = δ jj δ where δ is Kronecker's delta. It follows from (2.5) that the whole set of π k,j, (•; η) forms a CONB for H (f,η) and that any h(•) ∈ H (f,η) can be written as

(2.6) h(•) = ∞ k=0 (j, )∈B k θ k,j, (η)π k,j, (•; η),
where θ k,j, (η) = h(•), π k,j, (•; η) (f,η) . With (2.6) the problem of testing (2.3) reduces to that of testing the nullity of all (for an omnibus test within G ) or part (for a directional version) of the θ k,j, (η). Focusing on the latter, fix K and consider the subset of G that contains the functions of the form

(2.7) g K (•; η, θ(η)) = f (•; η) × 1 + K k=1 (j, )∈B k θ k,j, (η)π k,j, (•; η) 2 1 + K k=1 (j, )∈B k θ 2 k,j, (η)
.

Further confining ψ(•) to this subset reduces the problem of testing (2.1) to that of testing (2.8)

H 0 : θ(η) = 0 vs H 1 : θ(η) = 0, where θ(η) = {θ k,j, (η), 1 ≤ k ≤ K, (j, ) ∈ B k }. Under (2.7
), it is easy to see that if θ(η) is small, the dominant term in the maximum likelihood estimator (mle) of θ k,j,l (η) is n -1 n i=1 π k,j, (X i ; η). Hence, it makes statistical sense (for another justification, see [START_REF] Klar | Diagnostic smooth tests of fit[END_REF]) to base a test of (2.8) on

(2.9) π(η) = n -1 n i=1 π k,j, (X i ; η), 1 ≤ k ≤ K, (j, ) ∈ B k T .
If the π k,j, (•; •) are continuously differentiable in η, then Theorem 2.1 a) of [START_REF] Klar | Diagnostic smooth tests of fit[END_REF] ensures the asymptotic normality (with asymptotic expectation 0) of √ n π(η) under H 0 . The asymptotic covariance matrix depends on the choice of η. In particular, if f (•; η) satisfies the conditions in [13, p. 121], then η can be taken as the mle of η, which satisfies Assumption B with V 0 ( (X i , η)) = J -1 η , where J η is Fisher's information for η under H 0 . If J η denotes the matrix with elements Cov 0 (π k,j, (X i ; η), ∂ log f (X i ; η)/∂η a ) where η a is a a-th component of η, then from Theorem 2.1 c) of [START_REF] Klar | Diagnostic smooth tests of fit[END_REF], this covariance matrix is

I ν Q -J η J -1
η J T η where I v is the identity matrix of order ν and ν Q = K k=1 d m (k). Assume this matrix invertible. Then test statistic (2.10)

Q K = Q K (X, η) = n πT (η)(I v Q -J ηJ -1 η J T η ) -1 π(η)
is under H 0 asymptotically χ 2 ν Q . We refer to Q K as a global test statistic for H 0 . As a by-product of the CONB introduced in the next section, an explicit expression for this test statistic is available (see (4.1) and (4.2)).

For this to make statistical sense requires AL(m)-invariance of version (2.8) of problem (2.1) and AL(m)-invariance of test statistic Q K (X, η). The latter will be tackled in Sections 3 and 4. The former requires that, after transforming into γ(X), the new null θ(γ * (η)) = 0 holds if and only if θ(η) = 0. Because of the properties of the Π k in (2.5), a condition for this is that all π k,j, (•; η), (j, ) ∈ B k must appear in (2.7) with no θ k,j, (η) structurally set to zero (i.e. no term π k,j, (•; η) systematically excluded). As will be seen in Remark 5, some components of π(η) may systematically vanish, i.e. π k,j, (X i ; η) ≡ 0. Then all other terms in {π k,j, (•; η), (j, ) ∈ B k } must be dropped from (2.9) to preserve problem invariance. As will be seen in Section 4, the same condition, coupled with our choice for the π k,j, (•; η) also ensures AL(m)-invariance of Q K (X, η).

Remark 1. Proving the unicity of the mle η can be difficult. Some conditions are given in [START_REF] Kent | Redescending M-estimates of multivariate location and scatter[END_REF]. The verification of regularity conditions ensuring Assumption B (e.g. [START_REF] Ferguson | A course in large sample theory[END_REF]) is also tedious. Some work has been done in [START_REF] Hassan | Estimation dans les modèles elliptiques avec une application à la biostatistique[END_REF] for particular cases of the α-th power exponential distribution. Other estimators could in principle be used but the matrix in (2.10) becomes more complicated. Theorem 2.1.a) in [START_REF] Klar | Diagnostic smooth tests of fit[END_REF] gives its general expression without exploiting the particular structure of elliptical distributions.

Remark 2. In view of (2.10), the power function for a fixed ψ(•) can, under mild assumptions (see [START_REF] Inglot | Power approximations to and power comparison of smooth goodness-of-fit tests[END_REF]), be approximated by

P ψ [Q K > c] = P   K k=1 (j, )∈B k λ k,j, χ 2 1 (nν 2 k,j, ) > c   + O(n -1/2
), (2.11) where ν k,j, are functions of ´πk,j, (x; η)ψ(x)dx and λ k,j, are the eigenvalues of a complicated matrix (see [START_REF] Ducharme | Tree based Diagnostic Procedures Following a Smooth Test of goodness-of-Fit[END_REF]Sec. 3.2] for details in the case where J η = 0, as in Section 5.1). When ψ(•) is of the form (2.10), ν k,j, O(θ k,j, (η)). If one component of θ(η) in (2.8) differs from 0, the test will be consistent. This offers a handle on balancing the directional/omnibus slider. On the one hand, power is lost when, for k > K, one θ k,j, (η) is large. On the other hand, because θ k,j, (η) → 0 with k, a large K could add terms close to a χ 2 1 (0) in (2.11), leading to power dilution. To properly set a balance, some knowledge about ψ(•) is needed; otherwise see Remark 6 for choosing K in a data-driven fashion. But these remarks pertain to the θ k,j, (η) : a compounding difficulty is that some λ k,j, may also be small, with the corresponding λ k,j, χ 2 1 (nν 2 k,j, ) stochastically different from a χ 2 1 (0), which again affects power, in some cases rendering the test biased for small samples; see the α-th power exponential case in Section 6.

Remark 3. The above test can in principle be extended to the case where the density generator φ m (•) in (1.1) involves a shape parameter, such as with the α-th power exponential distribution. We do not pursue this further here as the necessary developments are beyond the scope of the paper.

Remark 4. Convergence toward the asymptotic χ 2 distribution can be slow. A benefit of invariance is that better approximations can be obtained by Monte Carlo resampling from any convenient f (•; η) distribution in F 0 .

3. The smooth test statistic.

Construction of the basis.

Suppose that X has a density of the form (1.1). The representation in [START_REF] Cambanis | On the Theory of Elliptically Contoured Distributions[END_REF] states that the random vector Y = V -1/2 (X-µ) has a spherical distribution (i.e. density (1.1) with η = (0, I m )) with stochastic representation

(3.1) Y = RU , where R = R(X) = V -1/2 (X -µ) is referred to as the radius, independent of U = U (X) = Y / Y ∼ U (Ω m
), the uniform distribution on the unit sphere Ω m . Here • is the Euclidean norm. The density of Y can thus be parametrized in polar coordinates as

(3.2) f R,U (r, u) = c m φ m (r 2 )r m-1 dr dω m (u),
where dω m (u) denotes the area element on

Ω m with ω m (Ω m ) = 2π m/2 /Γ(m/2). This mapping from X ∈ R m to (R, U ) ∈ R + × Ω m depends on η but f R,U (r, u) is invariant to the choice of V -1/2 .
A CONB for H (f,η) will be obtained by combining elements of CONBs associated with the distributions of R and U , so that π k,j, (X; η) can be written as π k,j, (R(X), U (X)).

First consider U . Let O(m) be the group of rotations on R m . From [17, p. 17], the space E m (k) of hyperspherical harmonics of degree k in m dimensions is an irreducible O(m)-invariant subspace of the space of homogeneous polynomials of degree k on Ω m and, from [1, Thm. 5.12, p. 81], the space of square integrable functions on Ω m can be decomposed as :

L 2 (Ω m ) = ∞ k=0 E m (k).
Again from [1, Prop. 5.8, p. 78],

(3.3) dim(E m (k)) = e m (k) =      1 if k = 0 m if k = 1 C m+k-1 m-1 -C m+k-3 m-1 if k ≥ 2 . Let {Ψ k, (u), = 1, . . . , e m (k)} be a CONB for E m (k) with respect to the scalar product h 1 , h 2 dωm = (ω m (Ω m )) -1 ´Ωm h 1 (u)h 2 (u)dω m (u). When m = 2, writing u T = (µ 1 , µ 2 ) = (cos θ, sin θ) we have Ψ k,1 (µ 1 , µ 2 ) = √ 2 cos(kθ), Ψ k,2 (µ 1 , µ 2 ) = √ 2 sin(kθ). When m = 3, Ψ k, (µ 1, µ 2 , µ 3 ) = P j k-j (µ 1 )Ψ k,j (µ 2 , µ 3
) where P j k-j (•) is the associated Legendre function of the first kind of order j and degree k -j, for j = 0, ..., k, j = 1, min{j +1, 2}. This is a particular case of a recurrence formula, explained in [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF], that generates the hyperspherical harmonics in dimension m from those in lower dimensions. Such calculations require symbolic manipulations. In particular, [1, Appendix B] give details about the MATHEMATICA program HFT10.m that computes {Ψ k, (u), = 1, . . . , e m (k)} for any m. To facilitate the application of the methods of the paper, the instructions to generate Ψ k, (u) with this program are detailed in our Appendix B.

Next, we look at the radius R. For each i ≥ 0, consider the set of functions {s j,i (r), j ≥ 0} where s j,i (•) is a polynomial of degree j in r 2 satisfying

(3.4) ω m (Ω m )c m ∞ 0 s j,i (r)s j ,i (r)φ m (r 2 )r m-1+2i dr = δ j,j .
To compute these polynomials, we apply the method described in Section 3 in [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF], but using the scalar product (3.4).

For i ≥ 0, let µ j,i = E 0 (R 2(j+i) ). Write µ j,i = (µ j,i , µ j+1,i , . . . , µ 2j-1,i ) T and M j,i =    µ 0,i • • • µ j-1,i . . . . . . . . . µ j-1,i • • • µ 2j-2,i    .
Set these quantities to 0 when j = 0. Using the argument leading to (4.4) of [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF], we get

(3.5) s j,i (r) = r 2j -(1, r 2 , . . . , r 2j-2 )M -1 j,i µ j,i µ 2j,i -µ T j,i M -1 j,i µ j,i
.

The MATHEMATICA commands for these functions in the important case of the MVN distribution are in Appendix B. The following theorem, whose proof is given in Appendix A.1, explains how to construct a CONB (i.e. the π k,j, (•; •) of the previous section) associated with an f (•; η) of the form (1.1).

Theorem 3.1. Fix η ∈ Ξ. Let P k be the space of polynomials of degree k in x ∈ R m . Let Π 0 = P 0 and for k ≥ 1, let Π k be the subspace of

P k such that p(•) ∈ Π k if (3.6) Rm p(x)p k (γ(x))f (x; η)dx = 0 for all p k (•) ∈ Π k , 0 ≤ k < k and γ ∈AL(m). Then, i)
The subspaces Π k are orthogonal with respect to •, • (f,η) . Moreover,

P k = Π 0 ⊕ Π 1 ⊕ • • • ⊕ Π k and (3.7) H (f,η) = ∞ k=0 Π k .
ii) The subspace Π k is an irreducible AL(m)-invariant subspace of the space P k . The polynomials

(3.8) π k,j, (R(x), U (x)) = R(x) k-2j s j,k-2j (R(x))Ψ k-2j, (U (x)), with (j, ) ∈ B k = {j = 0, . . . , [k/2], = 1, . . . , e m (k -2j)} form a CONB for Π k with respect to •, • (f,η) .
Here [z] denotes the integer part of z. Moreover

Π k is of dimension d m (k) = [k/2] j=0 e m (k -2j) = C m+k-1 k
. Finally each Π k can be further decomposed as

(3.9) Π k = [k/2] j=0 Π k,k-2j , where Π k,j is spanned by {π k,j, (R(x), U (x)), = 1, . . . , e m (k -2j)} and is an irreducible O(m)-invariant subspace of Π k whose dimension is e m (k-2j).
Beside giving the expression (3.8) for the π k,j, (•; η), this theorem is crucial for at least two other points. First, (3.7) ensures that an embedding family of the form (2.7) with K properly chosen and with (3.8) as the elements in the basis, can approximate almost any alternatives. Thus the test can be made adaptively omnibus; see Remark 6. Second, as will be seen in Section 4, (3.9) provides conditions ensuring the invariance of Dx tools. Note also that the differentiability of the π k,j, (•; η) required to apply [START_REF] Klar | Diagnostic smooth tests of fit[END_REF]'s results, follows from the structure in (3.8).

3.2.

The test statistic. The smooth test for null hypothesis (2.1) can now be explicited. Estimate η by the mle η = ( μ, V-1 ) with estimated information matrix J η. Set Ŷ i = V-1/2 (X i -μ) and compute Ri = Ŷ i , Û i = Ŷ i / Ŷ i . Inject these into the π k,j, ( Ri , Û i ) to compute πk,j, = n -1 n i=1 π k,j, ( Ri , Û i ). Group these into vector π(η) of (2.9) and compute J η . Test statistic (2.10) can be obtained. Null hypothesis (2.1) is rejected at approximate level α if Q K is greater than the 1 -α-th quantile of the reference distribution, e.g. the Monte Carlo approximation or the

χ 2 ν Q dis- tribution with ν Q = K k=k min C m+k-1 k . Remark 5.
As stated in Section 2, to ensure invariance care must be taken in selecting the elements of π(η). First, any π k,j, (•, •) that is a linear combination of some of the ∂ log f (X; η)/∂η a must be excluded because they will contribute nothing to the test statistic. Next, irreducibility of Π k dictates that to retain invariance, when π k,j, (•, •) has been excluded, then all other {π k,j , (•, •) (j , ) ∈ B k } must also be excluded. Remark 6. One difficulty in applying any smooth test resides in selecting the value of K. For univariate data, [START_REF] Ledwina | Data-Driven Version of Neyman's Smooth Test of Fit[END_REF]'s data-driven smooth test first chooses an integer d(n), performs a selection procedure to get a "good" K ∈ [1, ..., d(n)], and finally computes the associated test statistic. The framework where d(n) → ∞ leads to an adaptively omnibus test and is referred to as the infinite horizon case [48, p. 102]. The available theory about the rate of divergence of d(n) is impressive but hardly translates into a precise value, as it is expressed in terms of o(•). Thus d(n) is in practice fixed by external considerations, the finite horizon framework. Fortunately, simulations show that the power of the data-driven smooth test stabilizes rapidly as d(n) increases.Hence in practice both frameworks lead to the same modus operandi as long as d(n) is not too small.

The methodology that derives from Theorem 3.1 allows to compute the smooth test up to any desired d(n), so it is in principle possible to develop both infinite and finite horizon data-driven tests in our context. Here we only sketch a simple invariant adaptation of [START_REF] Ledwina | Data-Driven Version of Neyman's Smooth Test of Fit[END_REF]'s approach and focus on the finite horizon case. Suppose d(n) is given. Define

K = min{k : 1 ≤ k ≤ d(n), Q k -   k j=1 Card(B j )   log(n) (3.10) ≥ Q -   j=1 Card(B j )   log(n), = 1, ..., d(n)}.
The asymptotic reference distribution of test statistic Q K is a χ 2 with C m 1 degrees of freedom. As d(n) is allowed to increase, Q K spreads its power in a data-driven fashion over an increasing number of directions. The exploration of this and other scenarios is left for future work ; see Section 7.

4. Invariance and the Q = U IR decomposition. We now take a closer look at matrix (

I v Q -J η J -1 η J T η ) -1 in (2.10).
It is explained in Section A.2 that this is a matrix of constants and that after permuting the π k,j, (•, •) in (2.9) according to the values of k -2j, test statistic (2.10) can be conveniently written as (A.17), which leads to

Q K = Q K (X, η) = U K + I K + R K , (4.1)
which we refer to as the Q = UIR decomposition, where

U K = n πU 2 , I K = n πI,1 2 + d 1 tr(c 1 c T 1 Ī1 ĪT 1 ) + πI,2 2 + d 2 tr(c 2 c T 2 Ī2 ĪT 2 ) + πI,3 2 , (4.2) R K = n πR 2 + d 0 tr(c 0 c T 0 πR πT R ) ,
are asymptotically independent χ 2 with degrees of freedom ν U , ν I , ν R given in (A.2) under H 0 . Note that, beyond the moments of R 2 required for the s j,i (•), the computation of this Q = UIR decomposition conveniently requires only the 2 + 2 [K/2] + [(K + 1)/2] quantities described in Section A.2 to get the vectors

√ d 0 c 0 , √ d 1 c 1 and √ d 2 c 2 .
The elements π k,0, (r, u) in U K are proportional to r k Ψ k, (u). The components of πU are thus weighted averages of polynomials in Û i . Note that up to a constant, the π k,0, (r, u) are identical for any EC distribution and thus cannot be related to the distribution of R 2 . Hence U K serves to detect departures from the uniformity of U . When such departures are detected, the true density of the data may not be constant on ellipses centered at µ and a model incorporating this feature should be seeked, e.g. a distribution where U possesses a more complex density on Ω m .

Similarly, because Ψ 0,1 (u) ≡ 1, the elements of R K are π k,k/2,1 (r, u) = s k,0 (r) with k even. Hence this component, which under H 0 is asymptotically χ 2 ν R , serves to detect departures from the distribution of the "radius" R. These departures may then be identified, visually or otherwise, and the null model corrected accordingly. R K can also be used when, as in [START_REF] Kariya | LBI tests for multivariate normality in curved families and Mardia's test[END_REF], it is desired to have a test whose power is directed toward alternatives that are also elliptic about the unknown µ.

Finally, consider the elements in I K . These are products of polynomials in R with hyperspherical harmonics in U . Thus I K , which is approximately χ 2 ν I , detects correlations between U and R. When this occurs, the structure of the true density is complicated and a user could look at more involved densities, e.g. skew-densities of some sort.

Summing up, each component of the Q = UIR decomposition can detect a meaningful, in terms of stochastic representation (3.1), type of departure from the null density.

AL(m)-invariance of Q K and each of its components follows from the following considerations that hold provided all π k,j, (•, •) spanning each Π k , and thus each Π k,k-2j in (3.9), appear in (2.9). For any γ

= (A, b) ∈ AL(m), write X * i = γ(X i ), η * = γ * (η) = (μ * , V-1 * ) defined in Section 2. Also define Ŷ * i = V-1/2 * (X * i -μ * ) and similarly for R * i and Û * i . Observe first that Ri = R * i so that s k,0 ( Ri ) = s k,0 ( R * i ) and component R K , whose el- ements span the Π k,0 of Theorem 3.1, is trivially AL(m)-invariant. For component U K , whose elements span the O(m)-invariant Π k,[k/2] , notice that π k,0, ( R * i , Û * i ) = π k,0, ( Ri , O * Û i ), where O * = (A V-1 A T ) -1/2 A V-1/2 ∈ O(m). It follows that (π k,0, ( Ri , O * Û i ), = 1, ..., e m ([k/2])) T = O * * (π k,0, ( Ri , Û i ), = 1, ..., e m ([k/2]))
Û * i ) = Īa O * * a , it follows that I K is AL(m)-invariant. Then so is global statistic Q K as Q K (γ(X), γ * (η)) = Q K (X, η).
To be diagnostic, the components in (4.1) must be further processed and, in particular, [START_REF] Henze | Do components of smooth tests of fit have diagnostic properties ?[END_REF] shows that they must be scaled. This operation must be done with some care in order for the scaled statistics to retain both AL(m)invariance and their meaningful interpretations. One possibility is the following; let ΣUU be a block diagonal matrix where each block is the empirical covariance matrix of the π k,0, ( Ri

, Û i ) ∈ Π k,[k/2] appearing in π U . Define similarly ΣII , ΣRR for π I , π R . The scaled components U (s) K = nπ T U Σ-1 U U πU , I (s) K = nπ T I Σ-1 II πI and R (s) K = nπ T R Σ-1
RR πR are AL(m)-invariant and diagnostic. See Section 6 for some examples of their usefulness.

Remark 7. Statistic Q 1 is a function of π 1,0, (r, u) ∝ ru and is thus basically a distance between X and the mle μ. As such, it can be useful in discriminating H 0 but provides little Dx information in any senses associated with the Q = UIR decomposition. In the above, this term is bundled into component I K . Similarly Q 2 -Q 1 , which is spread over I K and R K , is related to tr( V -1 S) (S is the empirical covariance matrix of the X i ) and again can be useful in discriminating H 0 but otherwise provides little Dx insights. [START_REF] Csörgö | Consistency of some tests for multivariate normality[END_REF] advises that one should use a powerful test for the null hypothesis, followed by less formal procedures when the null is rejected. Thus it can be a reasonable strategy to consider only the elements Q k -Q 2 at the Dx stage.

5. Applications of the smooth test methodology to some EC distributions. The above smooth GoF test has the desirable features listed in Section 1: with K chosen appropriately, it can be made adaptively omnibus; it possesses Dx capabilities; and it is AL(m)-invariant. In addition, its behavior under both H 0 and H 1 can be conveniently approximated. In this section, we develop the test strategy for two popular EC distributions. Two others are treated in Appendix C. 5.1. The multivariate normal (MVN) distribution. Consider the important problem of testing the null hypothesis H 0 that a sample X 1 , . . . , X n arises from an m-dimensional MVN distribution with density

f (x; η) = (2π) -m/2 det(V -1 ) 1/2 exp -(x -µ) T V -1 (x -µ)/2 . Under H 0 , R 2 has a χ 2
m distribution with j-th moment 2 j Γ(m/2+j)/Γ(m/2). The elements of an orthonormal basis satisfying (3.4) are, via (3.5), (5.1)

s j,i (r) = (-1) j j! Γ(m/2) 2 i Γ(m/2 + j + i) L m/2+i-1 j (r 2 /2),
where L α j (•) is the j -th generalized Laguerre polynomial of order α; see Appendix B for MATHEMATICA commands to generate these quantities along with those for {Ψ k,j (•), j = 1, . . . , e m (k)}. Tables 3 and4 list all π k,j, (r, u) for k = 3, 4, 5 and m = 2, 3.

Next, it is easy to see that for the MVN, ζ(•) = 1 (defined in Section A.2) and c 0 , c 1 , c 2 vanish so that the Q = UIR decomposition (4.2) takes a particularly simple form. To get its explicit expression, first compute Ŷi , Ri and Û i with μ = X and V-1 = S -1 , where ( X, S) is the mle of (µ, V). Then all πk,j, = 0 when k = 1, 2 and in view of Remark 5, this leads to

Q 1 = Q 2 = 0. Upon writing C 2 k,j = n em(k-2j) =1
(π k,j, ) 2 , (2.10) becomes:

(5.2)

Q K = K k=3 [k/2] j=0 C 2 k,j , which is asymptotically χ 2 v Q under H 0 , where ν Q = K k=3 C m+k-1 k . Also (5.3) R K = K k=3 k:even C 2 k,k/2 is, under H 0 , asymptotically χ 2 ν R where ν R = [(K -2)/2]. Moreover, (5.4) U K = K k=3 C 2 k,0
has for reference distribution a χ 2 ν U , where ν U = K k=3 e m (k). Finally, (5.5)

I K = K k=3 [(k-1)/2] j=1 C 2 k,j . is approximately χ 2 ν I with ν I = ν Q -ν U -ν R .
To scale the components in the Q = UIR decomposition, let Σk,j denote the empirical covariance matrix of the n random vectors (π k,j, ( Ri , Û i ), = 1, . . . e m (k -2j)) with empirical mean πk,j = (π k,j, , = 1, . . . e m (k -2j)).

The scaled components are (C

(s) k,j ) 2 = n πT k,j Σ-1 k,j πk,j . These (C (s)
k,j ) 2 , bundled into U (s) , I (s) or R (s) are diagnostic to identify the aspects of the MVN not supported by the data. The R package ECGoFTestDx computing the above smooth test for the MVN has been deposited on CRAN. An application to a data set is shown in Appendix D.

When m = 1, Π k is spanned by π k,[k/2],1 (r, u), with u = ±1. These are the Hermite polynomials of [START_REF] Rayner | Smooth tests of goodness-of-fit using R, Second edition[END_REF]'s smooth test of univariate normality. Thus our approach generalizes their test.

For univariate distributions, the components of smooth test statistics are often related to GoF tests that have been introduced from other principles. The same occurs here and some components of Q K turn out to be wellknown GoF test statistics for the MVN. In particular,

Q 4 = C 2 3,0 + C 2 1,1 (= Q 3 ) + C 2 4,0 + C 2 4,1 + n(π 4,2,1 ) 2 . (5.6)
Inspection of these terms shows that R 4 = n(π 4,2,1 ) 2 = is nb 1,m /6 where b 1,m is [START_REF] Mardia | Measures of multivariate skewness and kurtosis with applications[END_REF]'s multivariate measure of skewness. Moreover, n(π 4,2,1 ) 2 = n(b 2,mm(m+2)) 2 /(8m(m+2)), where b 2,m is [START_REF] Mardia | Measures of multivariate skewness and kurtosis with applications[END_REF]'s multivariate measure of kurtosis. These are popular measures in multivariate analysis, whose usefulness in providing Dx information arises by extending the meaning of their univariate counterparts: in particular, b 1,m vanishes under elliptical symmetry and b 2,m quantifies the tails of the distribution. Consequently, they have been used as GoF tests for the MVN and, in particular, [START_REF] Romeu | A comparative study of goodness-of-fit tests for multivariate normality[END_REF] have concluded that tests based on them rank among the bests in terms of power. [START_REF] Koizumi | On Jarque-Bera tests for assessing multivariate normality[END_REF] have combined them (i.e. Q 3 +n(π 4,2,1 ) 2 ) to get a multivariate version of the popular Jarque-Bera test of univariate normality. [START_REF] Mardia | Rao score tests for goodness-of-fit and independence[END_REF] have developed a test of multivariate normality as an example of some simplifications that occur with Rao's score method when a suitable group structure exists. By brute force calculations, they were able to obtain Q 4 and state to have been unable to generalize to K > 4. Here, by making use of recent advances in computational tools for harmonic analysis, we can easily go beyond this limitation. Moreover, as a by-product, we get for any K the informative Q=UIR decomposition.

The idea of using the elements of stochastic representation (3.1) for testing MVN is not new. Many authors have proposed testing the distribution of R and/or U individually (see [START_REF] Koziol | Assessing multivariate normality: A compendium[END_REF][START_REF] Henze | Invariant tests for multivariate normality : a critical review[END_REF] for reviews), thus in effect crafting in an ad hoc fashion fragments of the Q=UIR decomposition. [START_REF] Koziol | An alternative formulation of Neyman's smooth goodness-of-fit tests under composite alternatives[END_REF], followed by [START_REF] Rayner | Smooth tests of goodness-of-fit for regular distributions[END_REF] and [START_REF] Best | A test for bivariate normality[END_REF], have taken a different route and proposed a smooth test based on the fact that the multivariate normal density is the product of m univariate normal densities when V = I m . For this particular case, a CONB can be obtained by the tensor product of elements of CONB associated with each univariate distribution [30, p. 51, Theorem 4.3]. However, this approach breaks down for general elliptical distributions, as in the following application, because V = I is not associated with independence [START_REF] Muirhead | Aspects of multivariate statistical theory[END_REF] and in view of the inconsistency property of many EC distributions.

The bivariate Laplace distribution.

There are several bivariate extensions of the univariate Laplace ; we take here the variant discussed in [START_REF] Naik | A Kotz-type distribution for multivariate statistical inference[END_REF] obtained from the power exponential distribution in Section 1 by setting α = 1/2. This distribution has larger tails than the MVN. The j-th moment of R 2 is Γ(2(j + 1)). We take K= 5, a reasonable value, and from (3.5), the s j,k-2j (•) required for Q 5 are :

s 0,1 (r) = ( √ 6) -1 ; s 0,2 (r) = (2 √ 30) -1 ; s 0,3 (r) = (12 √ 35
√ d0c0 = (7) -1 {- √ 614, 4 10 3 }, √ d1c1 = ( √ 1203) -1 { √ 4102, -4 2051 55 , 48 6 55 }, √ d2c2 = ( √ 19) -1 {-2 √ 35, 4}.
After computing the mle η = ( μ, V-1 ) and in turn, Ŷ i , Ri and Û i , plug these into the π k,j, ( Ri , Û i ) to get the πk,j, . Setting πU = (π 3,0,1 , π3,0,2 , π4,0,1 , π4,0,2 , π5,0,1 , π5,0,2 ), πR = (π 2,1,1 , π4,2,1 ), πI,1 = (π 1,0,1 , π3,1,1 , π5,1,1 , π5,2,1 , π 1,0,2 , π3,1,2 ), πI,2 = (π 2,0,1 , π4,1,1 , π2,0,2 , π4,1,2 ) and πI,3 = (π 5,1,2 , π5,2,2 ) yields all the elements required to compute the various test components. Their degrees of freedom are those in Section 4. For a GoF test for another bivariate Laplace distribution, see [START_REF] Fragiadakis | Goodness-of-fit tests for multivariate Laplace distributions[END_REF].

6. Performance of the Q = U IR decomposition in deriving Dx information. There can be many directions in which an alternative may depart from a null model. The present Dx tool, based on the Q = UIR decomposition, pertains to departures from the easily interpretable representation (3.1). To evaluate its performance, we need to determine how well it detects those departures, ideally in a context where other departures coexist to complicate matters. Performance here refers to the confidence allocated to the Dx information extracted from a sample, i.e. sensibility and specificity. As our tests are performed at a given level (here 5%), specificity is fixed and we concentrate on sensitivity, the probability that a departure from representation (3.1) will be detected, which is the power of the tests based on the components of the Q = UIR decomposition.

To start building this confidence, a small experiment was performed. We adopt the context of Section 5.1 where the interest lies in assessing the null hypothesis of bivariate normality. Because representation (3.1) is common to all EC distributions, the results below could be representative of what could be obtained with other EC null models and higher dimensions. Samples were generated from thirty distributions taken from [START_REF] Mecklin | A Monte Carlo comparison of the Type I and Type II error rates of tests of multivariate normality[END_REF][START_REF] Bogdan | Data driven smooth tests for bivariate normality[END_REF][START_REF] Johnson | Multivariate Statistical Simulation[END_REF]

. Statistics Q K , U (s) K , I (s)
K and R (s) K , (K = 3, ...,12) were computed and compared to their null Monte Carlo (based on 20,000 replications) approximation, as explained in Remark 4. To yield interesting power values,this was replicated 5,000 times and the sample size n was adapted to each alternative.

We stress that our goal is to appreciate the usefulness of the Dx procedure based of the Q = UIR decomposition and how these components relate to the global Q K . It is not to thoroughly compare the power of Q K to its competitors, as we have made no attempt to optimize the choice of K. Nevertheless, to offer some perspectives we did compute the power of the BHEP test [START_REF] Baringhaus | A consistent test for multivariate normality based on the empirical characteristic function[END_REF] with the tuning parameter set at 1.41; this test has been recommended in some simulation studies [START_REF] Mecklin | A Monte Carlo comparison of the Type I and Type II error rates of tests of multivariate normality[END_REF] but offers no Dx information. We also computed the power of the multivariate Jarque & Bera test [START_REF] Koizumi | On Jarque-Bera tests for assessing multivariate normality[END_REF] that does offer some Dx information via its skewness-kurtosis components.

The alternatives were chosen to somewhat resemble the bivariate normal and to depart along one or two directions identifiable by the Q = UIR decomposition, in addition to other departures. To screen the alternatives along these criteria, we generated 100 000 observations (X 1 , X 2 ) from each, transformed them into (Y 1 , Y 2 ) and then into (R, U ), and then again U into θ = arctan(U 2 /U 1 ). Plots of the various univariate and bivariate densities were made, and various statistics and characteristics were computed. We then excluded the densities where all UIR components are large, because the resulting knowledge (i.e. all is wrong) is not easy to use in the iterative process of correcting a null model after Dx analysis.

From the results from these thirty alternatives, we report here on five that cover most of the behavior we have observed. Also, we report here on the case K = 5.

The first alternative is a member of the Khintchine family of distributions [START_REF] Johnson | Multivariate Statistical Simulation[END_REF]Chap. 8]. Synthetic data from this distribution are obtained by generating Z ∼ Γ(1.5, 1) and then setting

(X 1 , X 2 ) = 3Γ(1.5) Γ(1.5 + 2 × 0.3998935) Z 0.398935 × 2(V 1 -0.5, V 2 -0.5),
where V 1 , V 2 are independent U (0, 1). In the above expression, the various constants are such that the marginals of (Y 1 , Y 2 ) are nearly independent N (0, 1) and the marginal density of R is approximately a χ 2 2 . However U is non-uniform, as the density of θ oscillates as sin(4θ) with modes at ±π/4, ±3π/4. The density of (Y 1 , Y 2 ) has squarish contours but otherwise resembles the bivariate normal. These oscillations induce a complicated relationship between R and U which, at the first order, can be assimilated to near independence, the regression of R on θ involving sin(kθ), cos(kθ), k = 1, ..., 3 being nearly constant with a Spearman's correlation coefficient (ρ S = -0.002) 1 . Because of these, we expect : power(R ). Note that this, and further, expectation derives from the 100 000 samples generated when screening the alternatives. This was blinded in the following simulation study. The sample size is n = 400.

The second distribution is the generalized Burr-Pareto-Logistic (with α = 1, β = 0) in [START_REF] Johnson | Multivariate Statistical Simulation[END_REF]Chap. 9]; see his Figure 9.8 for a sketch of the density and p. 167 for an algorithm to generate synthetic data. The distribution of (X 1 , X 2 ) has dependent N (0, 1) marginals with contours of triangular shape. Here the distribution of R is almost a χ 2 2 , U is clearly non-uniform with a trimodal distribution, while R and U are slightly correlated (with ρ S ≈ 0.04). Hence we should find : power(R [START_REF] Johnson | Multivariate Statistical Simulation[END_REF]Chap. 4]. R is again very close to 1 Pearson's correlation coefficient is not a good measure of the dependency between these quantities A normalized version of mutual information seems better adapted but we have found its scale difficult to appreciate from one problem to another. We have settled on Spearman's correlation coefficient mainly because it remains the same between θ and R or R 2 and because a sophisticated user can in general make some intuitive sense of differences between its values a χ 2 2 , U has a distribution that slightly differs from uniformity while R and U are slightly correlated (ρ S = -0.03). Again n = 400 and we expect to find : power(R The fourth density is a Laplace-type (because of its marginals) distribution generated by the following scheme : take

(s) 5 ) ≈ 5% < power(I (s) 5 ) ≤ power(U (s) 5 ). The sample size is n = 250. The third distribution is the contaminated binormal : 0.8×M V N 2 (0, I 2 )+ 0.2 × M V N 2 ((1, 1), Diag{1, 2}), see
W 0 , W 1 , W 2 ∼ Exp(1) and form (X 1 , X 2 ) = (W 1 -W 0 , W 2 -W 0 ).
Here, U is non-uniform, R departs from the χ 2 2 but R and θ are approximately independent (ρ S = -0.004), so one can expect: power(I

(s) 5 ) ≈ 5% ≤ power(U (s) 5 ), power(R (s) 5
). We have taken n = 75 to mitigate the high power of Q 5 with that of its components.

The last distribution is a 2-th power exponential distribution generated by taking R ∼ Γ(2, 2), θ ∼ U (0, 2π) and forming (X 1 , X 2 ) = R × (cos θ, sin θ). This yield an EC distribution where R is markedly different from a χ 2 2 . Thus here we should find : 5% ≈ power(U

(s) 5 ) ≈ power(I (s) 5 ) ≤ power(R (s) 5 ).
We have taken n = 100.

Table 1 shows the power (in %) of the various tests. The shaded cells are located where the power of the components of the Q = UIR decomposition are expected to be greater than nominal (5%). A first general conclusion is that no global test dominates. An interesting comparison is between J-B and Q 5 where, as stated in Section 5.1, the first is a fragment of the second. For some alternatives (e.g. Burr-Pareto-Logistic), the extra components have little power and lead to power dilution of Q 5 with respect to J-B. For other alternatives (e.g. Khintchine), power of J-B test is much lower because the extra components in Q 5 detect departures that hardly translate into skewness or kurtosis. A data-driven approach to select the value of K could be one way to further increase the power of Q K , but see Section 7. One atypical case is the 2-power exponential. This will be discussed below.

Regarding the power (i.e. sensitivity) of the components, the Dx information derived from the Q = UIR decomposition is consonant with the expectations related to each alternative. Furthermore, the magnitudes of the power in the shaded cells somewhat reflect the severeness of the departures (when grossly quantifiable). For the Khintchine, both R have power approximately equal to level and the severe non-uniformity of U is detected with good power. For the Burr-Pareto-Logistic, the expectation : power(R

(s) 5 ) ≈ 5% < power(I (s) 5 ) ≤ power(U (s) 5
) is realized while the small power of I 1 Power in % (based on 5,000 replications at level α = 5%) of the test of bivariate normality based on Q5 and its UIR components. The reference distribution has been approximated by 20,000 Monte Carlo samples of size n from the M V N (0, I2). Alternatives are K=Khintchine, BPL=Burr-Pareto-Logistic, CMVN=Contaminated MVN, Lt=Laplace -type and 2PE=α = 2-power exponential. Also shown is the BHEP (Baringhauss-Henze-Epps-Pulley) test and the J-B (multivariate Jarque-Bera) test along with its components β1, β2 (both using Monte Carlo quantiles). The shaded cells correspond to those where we expect power > 5%.

of U is detected with moderate power. For the Laplace-type case, the departures on the marginal distributions of R and U are correctly detected. For the 2-power exponential, power

) ≈ 91.6, showing that the departures along the distribution of R are correctly captured. The two other components have power close to nominal, as they should.

To summarize, one can have some confidence that our Dx procedure can be a useful tool in the iterative process of modeling a data set.

In some cases, at least one component yield better power than the global test statistic from which it is extracted. Thus the natural two-stage strategy where one first test the global null hypothesis using Q K and, if significant, proceed at the Dx level, would lose valuable information. This can be explained by a number of factors. First, when only one component ought to be significant, as with the Khintchine, the degrees of freedom associated with this component (6 in this case) is much smaller than that of the global test (here 15), thus diluting power. This effect is further confounded by the scaling of the components, which also affects power, and suggests looking at a scaled global statistic. Here power(Q

(s) 5 = U (s) 5 + I (s) 5 + R (s)
5 ) = 83.5, greater than the other global tests, which suggest replacing Q 5 by its scaled version. However the benefit of this replacement is not universal : for the contaminated normal, power(Q 5 ) = 54.5 while power(Q (s) 5 ) = 45.6. Thus it is not clear which of the above factors (degrees of freedom vs scaling) has more impact on power. An extreme case of this is the 2-th power exponential where the power of Q 5 is almost zero and the J-B test is not much better. Here scaling has a large effect on power as power(Q (s) 5 ) = 47.5. This is an instance where the values of the λ k,j, in Remark 2 renders the global test biased. Scaling partly corrects the problem, as explained in [START_REF] Ducharme | Tree based Diagnostic Procedures Following a Smooth Test of goodness-of-Fit[END_REF]. A pro-cedure that would select in a data-driven fashion between Q K and Q (s) K , in addition to a proper value of K globally or at the component level, would perhaps reduce such variations. See Section 7.

The J-B test also offers some Dx information via its b 1,m , b 2,m components. But they are not scaled and, in particular, R (s) 5

is the scaled version of b 2,m . Comparing their powers shows the pitfalls in attempting to interpret unscaled Dx statistics. However, being fragments of our components, they can be rescaled and exploited to further refine the Dx information, a problem we now look into from another angle.

In the above, we have taken K = 5 and shown that our procedure can provide useful Dx information supplementing the global test. Now we go further and vary K to see if more precise Dx information can be extracted from a set of components. Table 2 presents further results from our experiment, namely the values of (U 4 pick up some important departure from uniformity. To see what as been detected, we refer to the definition of projection-based Dx information in [START_REF] Ducharme | Tree based Diagnostic Procedures Following a Smooth Test of goodness-of-Fit[END_REF]. Recall that the elements in U (s) 3 have the form sin(kθ), cos(kθ), k = 1, 2, 3 and thus have at most 3 maximums in [0, 2π). Those added in U (s) 4 have the form sin(4θ), cos(4θ) and the increased power indicates that they capture higher oscillations, suggesting that the density of U (or θ) could have at least four maximums. This is consonant with the fact that θ has a distribution that oscillates as sin(4θ). Also, recall that the scatter plot of (R, θ) reveals that these quantities have a complicated relationship. A similar increase of power from I (s) 5 to I (s) [START_REF] Boulerice | Smooth tests of goodness-of-fit for directional and axial data[END_REF] and afterwards indicates again that some correlation has been detected at the order of functions of the form cos(6θ), sin(6θ). Indeed a regression of R on θ with basis functions sin(kθ), cos(kθ), k = 1, ..., 7 retains the functions cos(4θ) and sin(6θ) as significant. Hence the results in Table 2 indicate that beyond the basic analysis of the Q = UIR components, a finer analysis of the subcomponents could be exploited to pinpoint with better precision the nature of the departures.

7.

Conclusions. This paper develops a smooth test of goodness-of-fit for EC distributions that, through the invariant Q = UIR decomposition, allows to extract illuminating Dx information. Some applications are worked out, in particular the important case of the MVN, where the usefulness of the Q = UIR decomposition is demonstrated. For most EC distributions, there exists no competitor to our approach. For the MVN, a small experiment Proof of i ). A change of variable in 3.6 shows that, without loss of generality, we can set η = (0, I m ). By (3.6), when k = k , the spaces Π k and Π k are orthogonal with respect to •, • (f,η) . By the lemma, all polynomials p(•) in P k have an expansion of the form (A.1) where from (A.8) we can take γ = e to simplify the derivations, while all p k (γ(•)) ∈ Π k have an expansion of the form A.3. Substituting these expansions into (3.6) shows, after careful identification and using the orthonormality of the s j,i (•), Ψ i, (•), that p(•) ∈ Π k if and only if its coefficient in (A.1) satisfies a j,i, (e) = 0 for every i, j such that i = k -2j. Thus, after rearranging the indices and extending to any γ ∈ AL(m), each p(•) ∈ Π k has an expansion of the form (A.9)

p(x) = [k/2] j=0 em(j) =1 a j,k-2j, (γ)s j,k-2j ( γ -1 (x) ) Ψk-2j, (x; γ -1 ).
Because the set of polynomials s j,k-2j ( x ) Ψk-2j, (x) are orthonormal with respect to •, • (f,η) and span Π k , it follows from (A.1) and (A.9) that P k = Π k ⊕ P k-1 . This proves the first part of i ). The second part comes from the fact that the polynomials are dense in L 2 AL(m) (f, η). Proof of ii ) By construction, Π k is an AL(m)-invariant subspace. From (A.9), for any γ ∈ AL(m), there is only one element of Π k , up to a multiplicative constant, that is invariant over the subgroup of rotations. It follows from Schur's lemma that Π k is irreducible. Indeed, Schur's lemma [17, p. 390] asserts that if two representations of a group are irreducible then every commuting non-zero linear map between them is an isomorphism and that a representation is irreducible if and only if every such linear map is a scalar multiple of the identity element. In the present context, showing that the space Π k is irreducible amounts to show that the representation of AL(m) induced by the coefficients a k,j, (γ) of its basis in (A.9) is irreducible. The identity element is then the element that is invariant over rotations which are identified with the commuting maps.

Setting

g -1 = (V -1/2 , -V -1/2 µ), the π k,j, (r, u) = s j,k-2j ( g -1 (x) ) Ψk-2j, (x, g -1 )
are orthonormal with respect to •, • (f,η) and, from (A.9), form a complete basis for Π k . Let Π k,j (g -1 ) be the space generated by the span of {s j,k-2j ( g -1 (x) Ψk-2j, (x; g -1 ), = 1, . . . , e m (k -2j)}. Because, for fixed (k, j), Π k,j (g -1 ) is closed under rotations, the functions that span Π k,j (g -1 ) are, up to multiplicative constants, the same as those that span E m (k -2j). Thus Π k,j (g -1 ) is an irreducible rotation-invariant subspace of Π k of dimension e m (k-2j). Combining these functions shows that for g -1 = (V -1/2 , -µ),

Π k = Π k,0 (g -1 )⊕• • •⊕Π k,[k/2] (g -1 )
, with respect to •, • (f,η) . The dimension of Π k is computed from this direct sum as a combinatorial exercise.

A.2. Derivation of the Q = U IR decomposition. We derive (4.1) and (4.2). First reparametrize η = (µ, V -1 ) into the m + m(m + 1)/2 vector η = (µ, lvec(V -1 )) T , where lvec(•) is the "lower vec" operator. Introduce the m 2 × m(m + 1)/2 matrix B with element :

b ij,k = 1 if (i = k and j = ) or (i = and j = k) 0 otherwise ,
where i first goes from 1 to m and then j goes from 1 to m while k varies from to m followed by varying from 1 to m. This matrix, which has rank m(m+1)/2, is such that vec(V -1 ) = B•lvec(V -1 ). Now some algebra yields

∂ log f (x; η) ∂η =N η Z,
where (⊗ denotes Kronecker's product)

N η = V -1/2 0 0 1 2 B T (V -1/2 ⊗ V -1/2 )B , Z = Z 1 Z 2 = R × g(R 2 ) × U lvec(I m -ζ(R 2 )U U T ) , (A.10)
with g(r 2 ) = -2φ m (s)/φ m (s)| s=r 2 and ζ(r 2 ) = r 2 × g(r 2 ). Note that N η is invertible and the distribution of Z does not depend on η; it has expectation 0 and, in view of results about the U (Ω m ) distribution (see [START_REF] Bilodeau | Theory of multivariate statistics[END_REF]p. 239,Ex. 4], its variance is block-diagonal with the inverse of these blocks being :

V -1 (Z 1 ) = mI m /σ 1 , V -1 (Z 2 ) = σ 2 [I m(m+1)/2 - 1 2 Diag(lvec(I m ))] - σ 2 (1 -σ 2 ) 4 + 2m(1 -σ 2 ) lvec(I m )lvec T (I m ), (A.11)
where

σ 1 = E 0 (R 2 × g 2 (R 2 )) and σ 2 = m(m + 2)/E 0 (ζ 2 (R 2 )). Then J η = N η V(Z)N T η is also block-diagonal. Moreover, the lines of J η are Cov 0 π k,j, (R, U ), ∂ log f (X; η)/∂η T = Nη E 0 R k-2j s j,k-2j (R 2 )g(R 2 )Ψ k-2j, (U )Z T .
Hence the elements of J η J -1 η J T η do not involve η and (

I v Q -J η J -1 η J T η ) de- pends only on the distribution of R. Defining Ψ k = (Ψ k, (U ), = 1, . . . , e m (k)) T , it follows that U = Ψ 1 / √ m while lvec(U U T ) = A m (Ψ T 0 , Ψ T 2 ) T ,
E 0 R k-2j s j,k-2j (R 2 )Ψ k-2j Z T 1 = E 0 s j,1 (R 2 )ζ(R 2 ) Im/ √ m if k -2j = 1 0 otherwise , (A.12)
while adding the orthogonality of the s j,i (•),

E 0 R k-2j s j,k-2j (R 2 )Ψ k-2j Z T 2 =      -E 0 s j,0 (R 2 )ζ(R 2 ) A T m(1) if k -2j = 0 -E 0 R 2 s j,2 (R 2 )ζ(R 2 ) A T m(2) if k -2j = 2 0 otherwise . (A.13)
Collect the constants in these expressions into the following vectors : for

k -2j = 0, c T 0 = -E 0 s j,0 (R 2 )ζ(R 2 ) , j = 1, ..., [K/2] , for k -2j = 1, c T 1 = E 0 s j,1 (R 2 )ζ(R 2 ) / √ m, j = 1, ..., [(K + 1)/2] and for k -2j = 2, c T 2 = -E 0 R 2 s j,2 (R 2 )ζ(R 2 ) , j = 1, ..., [K/2] . Permute the components of (π k,j, (•, •), k = 1, ..., K, (j, ) ∈ B k ) T into (π T U , π T I , π T R ) T so that all π k,0, (•, •) with k ≥ 3 are in π U in increasing order of and then k, all π k,[k/2],1 (•, •)
with k ≥ 2 are in π R while all other terms are regrouped in π I in the following way : let I 1 be the matrix with rows (π k,j, (r, u) | k -2j = 1, = 1, ..., e m (k -2j)) ordered from row to row according to k. Similarly define I 2 with k -2j = 2 and regroup all other terms (those with k -2j > 2) into π I,3 ordered lexicographically. Then π T I = (vec(I 1 ) T , vec(I 2 ) T , π T I,3 ) T = (π T I,1 , π T I,2 , π T I,3 ) T . Form (π T U , πT I , πT R ) T , the permuted version of (2.9) ordered in the same way with πT I = (π T I,1 = vec( Ī1 ) T , πT I,2 = vec( Ī2 ) T , πT I,3 ) T . Easy calculations show that J η with rows permuted accordingly has the form (A.14) Table 3 lists the polynomials π k,j, (r, u) for the case m = 2, k = 3, ..., 5 that are used in the experiment of Section 6.

J η = N η         0 ν U × m 0 ν U × m(m + 1)/2 τ 11 ν I × m τ 12 ν I × m(m + 1)/2 0 ν R × m τ 22 ν R × m(m + 1)/2         , where ν R = [K/2] , ν U = K k=1 e m (k) , ν I = ν Q -ν R -ν U where v Q = K k=1 C m+k-1 k . Recall that the τ ij depend only on the distribution of R π3,0,1 = r 3 u 1 (3-4u 2 1 ) 2 √ 6 π3,0,2 = r 3 u 2 (1-4u 2 1 ) 2 √ 6 π3,1,1 = r(r 2 -4)u 1 2 √ 2 π3,1,2 = r(r 2 -4)u 2 2 √ 2 π4,0,1 = r 4 (1-8u 2 1 +8u 4 1 ) 8 √ 3 π4,0,2 = r 4 u 2 (u 1 -2u 3 1 ) 2 √ 3 π4,1,1 = r 2 (r 2 -6)(1-2u 2
Table 4 lists the π k,j, (r, u) required for the Open/Closed book example of Section 6 (m = 3 and k = 3, 4 and 5 )

APPENDIX C: SMOOTH TEST FOR THE BIVARIATE LOGISTIC AND PEARSON TYPE II DISTRIBUTIONS

C.1. The bivariate logistic distribution. We consider the bivariate logistic distribution in [START_REF] Lemonte | Multivariate elliptical models with general parametrization[END_REF] with density generator φ m (y) = e -y /(1 + e -y ) 2 . This is again a competitor to the MVN but with shorter tails. The moments of R 2 have explicit but complicated expressions and it is shorter to use numerical approximations. Consequently, the s j,k-2j (•) required for Q 5 are : s 0,1 (r) = 0.849322; s 0,2 (r) = 0.551329; s 0,3 (r) = 0.30403; s 0,4 (r) = 0.148319; s 0,5 (r) = 0.0654688; s 1,0 (r) = -1.18523 + 0.854964 r 2 ; s 1,1 (r) = -1.36758 + 0.576276 r 2 ; s 1,2 (r) = -1.0461+0.318116 r 2 ; s 1,3 (r) = -0.646026+0.153749 r 2 ; s 2,0 (r) = 1.24468-1.86588 r 2 +0.407913 r Continue as in the bivariate Laplace case.

C.2. The bivariate Pearson type II distribution. We consider the bivariate Pearson type II distribution described in [START_REF] Johnson | Multivariate Statistical Simulation[END_REF]Sec 6.2] with density generator φ m (y) = (1 -y) α with y ∈ [0, 1]. For α > 0, this is another EC distribution that somewhat resembles the MVN. The j-th moment of R 2 is (2j + α)B(2 + α, 1 + j) and from (3.5), the s j,k-2j (•) required for Q 5 are : s 0,j (r) = ( (j + 1)B(2 + α, j + 1) -1 for j = 1, ..., 5; s 1,j (r) =

√ α + 6 √ 3 }, √ d2c2 = - √ α -1 2 √ 2 { (3 + α) √ α + 1 √ 3 , √ 5 + α}.
Continue as in the previous cases. Note that here α is supposed known, but the above can serve to prefigure the difficulties to be encountered with an unknown shape parameter.

APPENDIX D: AN APPLICATION

As an application of the methods of the paper, consider the "Open-book closed-book examination" data set [34, p. 3-4] corresponding to examination marks in Mechanics, Vectors, Algebra, Analysis and Statistic for a sample of 88 students. Here we consider only the marks in Vectors, Algebra and Statistics, so that m = 3 and we wish to test a trivariate MVN. We apply 5 Examination marks (n = 88) in Vectors, Algebra and Statistics from the "Open book-Closed book examination" data set [34, p. 3-4]. The null hypothesis is a trivariate normal distribution; d.f. refers to degrees of freedom of the asymptotic χ 2 approximation, p-value χ 2 are computed from these reference distributions, while p-value MC refers to the p-value computed from a Monte Carlo approximation (20 000 replications)

our test strategy with K = 5. All polynomials required for these computations appear in Table 4. The results are shown in Table 5 which lists test statistic Q 5 and its scaled UIR components, as well as the p-values obtained from the reference asymptotic χ 2 distribution and from Monte Carlo (20 000 replications) approximations. The null hypothesis of trivariate MVN is rejected at the 5% level. Inspection of the scaled components shows that both the distributions of R and U appear consonant with the MVN and that rejection comes from correlations between R and U . If the confidence build in this section for m = 2 can be transferred to the case m = 3, a refined model should try to take into account the dependencies between these random quantities.

) - 1 ;

 1 s 0,4 (r) = (72 √ 70) -1 ; s 0,5 (r) = (720 √ 77) -1 ; s 1,0 (r) = (r 2 -6)/(2 √ 21); s 1,1 (r) = (r 2 -20)/(4 √ 165); s 1,2 (r) = (r 2 -42)/(60 √ 42); s 1,3 (r) = (r 2 -72)/(144 √ 665); s 2,0 (r) = (7r 4 -360r 2 + 1320)/(24 √ 10745) and s 2,1 (r) = (11r 4 -10926r 2 +12600)/(72 √ 225610). Along with the expression for Ψ k, (•) in Section 3.1 when m = 2, it is an easy exercise to compute the π k,j, (•, •). Next, because ζ(r 2 ) = r 2 × 1 r 2 , one finds σ 1 = 1, σ 2 = 6 and

5 ( 21 .

 521 8%) reflects the slight correlation between R and U . Again, the severe non-uniformity of U is almost always detected. For the contaminated MVN, R (s) 5 has, as expected, trivial power, I (s) 5 picks up rather well the dependency between R and U while the slight non-uniformity

4 .

 4 (s) k , k = 3, ..., 7) and (I (s) k , k = 3, ..., 7) in the case of the Khintchine distribution (n = 400). Notice the sharp increase in power between U This indicates that the hypersphericals in the added components of U (s)

1 m 2 - 1 ,

 121 where A m is a square matrix of order m(m + 1)/2) partitioned as (A m(1) A m(2) ) with A m(1) = lvec(I p )/m being its first column. The expression of A m(2) is not required beyond the relationships lvec(I p ) T A m(2) = 0 and A T m(2) A m(2) = (m+2) I m(m+1) from the expression for V(lvec(U U T )) in [4, p. 239, Ex. 4]. It follows from the orthogonality of the Ψ k-2j, (•) that

1 ) 4 √ 3 π4, 1 , 2 = r 2 (r 2 - 6 )u 1 u 2 ) 2 √ 3 π4, 2 , 1 = r 4 -8r 2 +8 8 π5 6 π5, 2 , 1 = r(r 4 -12r 2 +24)u 1 8 √ 3 π5, 2 , 2 = r(r 4 -12r 2 +24)u 2 8 √ 3 3

 1431222623214862142132242233 TableTableof the polynomials π k,j, (r, u) of Theorem 3.1 in the case of a MVN null hypothesis, for m = 2 and k = 3, 4 and 5.generates the basis of dimension e m (k) for E m (K). For example, typingsetDimension[u,5]; Do[Print[TableForm[Ψ(k, u)]], {k,1,4}];prints the {Ψ k,j (u) | j = 1, . . . , e 5 (k)}, for k = 1, ...,4. To generate the s j,i (r) for the case of the null MVN distribution, the command iss[j , i , r ] := (-1)^j*Sqrt[j!*Gamma[m/2]/(2^i*Gamma[m/2 + j + i])]* LaguerreL[j, m/2 + i -1, r^2/2];Finally, to generate the π k,j, (r, u), the commands aree[m , k ] := Which[k == 0, 1, k == 1, m, k >= 2, Binomial[m + k -1, m -1] -Binomial[m + k -3, m -1]]; πfunct[k , j , , r ] := r^(k -2*j)*s[j, k -2*j, r]*Ψ[k -2j, u][[ ]]; CONBΠ[k ] := Flatten[Table[πfunct[k, j, , r], {j, 0, Floor[k/2]}, {l, 1, e[m, k -2*j]}]];

4

 4 and s 2,1 (r) = 1.79127-1.51119 r 2 + 0.230011 r 4 . Next, ζ(r 2 ) = 2r 2 tanh(r 2 /2), σ 1 = 3.18173, σ 2 = 0.82306 and √ d0c0 = {-9.45511, -0.77618}, √ d1c1 = {12.55, 3.40145, -1.89893}, √ d2c2 = {-8.04922, -1.2599}.

(. 1 ,

 1 2+j+α)r 2 -(1+j) √ B(2+α,j+2)(1+α)(2+j+α) for j = 0, ..., 3; s 2,0 (r) = (2+r 2 (3+α)(r 2 (4+α)-4) Next, because ζ(r 2 ) = 2αr 2 1-r 2 , one finds σ 1 = 4α(α+1)α-σ 2 = 1 -α -1 (hence α > 1 in the sequel) and

  T , for some O * * ∈ O(e m ([k/2]) by standard properties of Wigner d-matrices. Hence each (π k,0, , = 1, ..., e m ([k/2])), and thus U K , is invariant to AL(m)transformations. A similar argument applies to the parts πI,j 2 , j = 1, 2, 3 of component I K which span the intermediate O(m)-invariant Π k,k-2j . Finally, consider the term tr(c a c T a Īa ĪT

a ), a = 1, 2. Because Īa averaged over all ( R * i ,

  Test Statistics Value d.f. p-value χ 2 p-value MC

	Q5	98.62	46	0.004	0.027
	U 5 (s)	37.38	27	0.088	0.163
	I	(s) 5	33.14	18	0.016	0.028
	R	(s) 5	0.77	1	0.381	0.466
					Table

G. DUCHARME ET AL.

 2Power in % (based on 5000 replications) of some components in the UIR decomposition for the Khintchine alternative. The level is α = 5%. The reference distribution has been approximated by 20 000 Monte Carlo samples of size n = 400 from the M V N (0, I2).

suggests that our smooth test can be competitive powerwise.

In the present work, no attempt is made to optimize the selection of the hyperparameter K. his topic deserves a separate analysis because the present context opens up many new possibilities beyond the obvious extension to select among Q 1 , Q 2 , ..., Q K in a data-driven fashion via (3.10). One of these is to increase sensitivity of the Dx by selecting among the components, say U 1 , U 2 , ..., U K . In addition, it was shown in Table 2 that this could also prove useful in better understanding the nature of the detected departure. Another is to select among {U K , U K +I K , U K +I K +R K = Q K } or some other (datadriven) permutation of the components. Yet another is to select between the raw components and their rescaled versions.

In multivariate contexts, the computation of the mle may be difficult or the information J η may not exist. This is the case for example in the Laplace discusses by [START_REF] Fragiadakis | Goodness-of-fit tests for multivariate Laplace distributions[END_REF], thus forcing the use of moment estimators. Deriving the Q = UIR decomposition in such cases remains to be investigated.

Finally, EC distributions often involve (e.g. the α-th power exponential), in addition to η, a shape parameter that offers more flexibility in adjusting the data. The shape parameter must in general be estimated and this does not, in principle, affect the use of a test statistic similar to (2.9). However, an open question is how does this estimation impact the Q = UIR decomposition.

APPENDIX A: PROOFS AND AUXILIARY RESULTS

A.1. Proof of Theorem 3.1. Let k ≥ 0 and x∈ R m . For any fixed γ = (A, b) ∈ AL(m), define the functions Ψi, (x; γ) = Ψ i, (γ(x)/ γ(x) ) γ(x) i , where the Ψ i, (•) are the spherical harmonics evoked in Section 3.1. To avoid trivialities, set Ψi, (-b; γ) = 0 for i > 0 and Ψ0,1 (-b; γ) = 1. It follows that Ψi, (x; γ) is an homogeneous polynomial of degree i in the components of γ(x) = A(x + b) and thus a polynomial of degree i in x ∈ R m . We need the following lemma which states that for any given γ, any polynomial in R m can be expressed as a linear combination of the π k,j, (r, u) of (3.8) and explains how the coefficients of this linear combination are affected by the choice of γ.

Proof. Let x = ru where r = x and u = x/ x . To avoid trivialities again, when x = 0, set r = 0, define u arbitrarily and agree to take r i = 1. Consider the polynomial p(γ(ru)) on R + × Ω m . For fixed r, this is a polynomial of degree at most k in the components of u. Now, because the spherical harmonics {Ψ i, (•), = 1, . . . , e m (i), i ≥ 0} form a CONB with respect to the scalar product •, • dωm for L 2 (Ω m ) , and thus for the restriction to Ω m of the space of polynomials of degree k on R m , it follows that

Setting e = (I m , 0) as the identity element of AL(m), this shows that

where w i, ( x ; γ) = x -i q i, ( x ; γ). Thus

Now, because Ψi, (x; γ -1 ) is a polynomial of degree i in x and p(γ(x)) is a polynomial of degree at most k, it follows that each w i, ( x ; γ) is a polynomial of degree at most [(k -i)/2] in x 2 . Hence in turn, it can be written as a linear combination of the polynomials s j,i ( x

Combining (A.4) and (A.5) yields (A.1). Now using (3.2) and (3.4) with (A.5), we have

upon transforming from polar to cartesian coordinates. Finally, substitute (A.4) into the right-hand side of (A.2). Simple algebra using the orthogonality of the Ψ i, (•) yields the right-hand side of (A.8). This concludes the proof of the lemma.

under H 0 . Thus

We focus on Σ IR = -τ 12 V -1 (Z 2 )τ T 22 . The lines of matrix τ 22 have the form -E 0 s j,0 (R 2 )ζ(R 2 ) (lvec(I p )) T . Combining this with the expression of V -1 (Z 2 ) and the fact that the lines of τ 12 either vanish or, from (A.13) have the form

. Using standard properties of the vec and operators, we get

) and similarly for the other term. Collecting these, we finally get : We give the MATHEMATICA commands to generate the {Ψ k,j (u) | j = 1, . . . , e m (k)} of Section 3.1. The package HFT10.m must first be downloaded from the site given in Appendix B of [START_REF] Axler | Harmonic function theory[END_REF] and loaded into MATHEMATICA via the command <<''.../.../.../HFT10.m'', where .../.../.../ is the path leading to where the package HFT10.m is stored on the computer; after typing <<, one can use the "File Path..." command in the "Insert" menu to automatically generate this path.

Once the package is loaded (the text "* You can now use the functions in this package." will appear on the MATHEMATICA notebook), the user needs only to set the dimension m via the command setDimension[u, m]. Here u = (u 1 , ..., u m ) is the vector of variables in which the spherical harmonics will be expressed. Then the command

π3,0,3 =

π5,0,7 =

π5,0,8 =

2 ) 6 π5,0,9 =


Alexandr Klimchik
email: alexandr.klimchik@emn.fr

David Daney
email: david.daney@inria.fr

Stephane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Anatol Pashkevich
email: ana-tol.pashkevich@emn.fr

Geometrical patterns for measurement pose selection in calibration of serial manipulators

Keywords: calibration, design of experiments, geometrical patterns, serial robot

The paper is devoted to the accuracy improvement in geometric calibration of serial manipulators. Particular attention is paid to the optimal selection of measurement poses, which reduce measurement noise impact on the parameters identification precision. In contrast to previous works, the proposed approach yields simple geometrical patterns that allow user to take into account the joint and workspace constraints and to find measurement configurations without tedious computations. The advantages and practical significance of the proposed approach are illustrated by an example that deals with 6-dof serial manipulator.

Introduction

Since most of industrial serial robots operate without position feedback, the precision of a geometric model used in the control algorithm should be high enough. In order to meet this requirement in practice, the calibration technique is usually applied. For this reason, the problem of robot calibration has always been in the focus of the robotic community. There exist open-loop and closed-loop techniques that are suitable for serial and parallel manipulators [START_REF] Elatta | An overview of robot calibration[END_REF][START_REF] Hollerbach | Model identification[END_REF]. However, to apply them in a real industrial environment, it is necessary to reduce the impact of the measurement noise that corrupts the input data and is evaluated via the terms of the covariance matrix [START_REF] Mooring | Fundamentals of manipulator calibration[END_REF].

In order to reduce the impact of measurement noise, two approaches can be generally applied: (i) increasing the number of experiments and (ii) applying design of experiment theory for selection of measurement configurations. The first approach is quite a simple, it does not require any computation, but increases a lot the measurement time and related cost. In contrast, the second approach allows us to increase identification accuracy without increasing the number of measurements [START_REF] Borm | Determination of optimal measurement configurations for robot calibration based on observability measure[END_REF]. However this approach requires additional knowledge and some computational efforts. Nevertheless, because of its advantages the second approach is more attractive and is in the focus of this paper.

Although the design of experiment theory has been used for more than a century in different areas, the obtained results cannot be applied directly in robotics since they are mostly suitable for linear regression models. Therefore, to obtain similar results for robot calibration that rely on highly non-linear models, the existing techniques should be essentially revisited. In addition, to make the technique attractive for practicing engineers, it is reasonable to obtain some geometrical patterns allowing the user to generate optimal measurement configurations without any computational efforts.

Problem of measurement pose selection in robot calibration

To show the importance of the problem of interest, let us present first a simple motivation example. For the purpose of simplicity, let us limit our study by a 2-dof planar manipulator with two actuated revolute joints and link lengths l 1 = 1.0 m, l 2 = 1.0 m . This manipulator has four geometrical parameters to be identified: (i) two link length deviations and (ii) two joint encoder offsets.

For comparison purposes, let us examine two sets of measurement configurations assuming that the measurement noise is Gaussian with zero mean and the standard deviation σ = 0.1 mm. Simulation results are summarized in Fig. 1, where the rootmean-square errors ρ of the end-effector position after calibration throughout the robot workspace are presented. As follows from Fig. 1, in the worst manipulator configuration, the Set #1 provides a positioning error equal to 2.29 mm while the Set #2 reduces the worst positioning error down to 0.14 mm, i.e. by a factor of 16. Hence, this simple example clearly shows that the selection of measurement configurations is a very important issue in robot calibration. In fact, poorly chosen measurement configurations may have a negative effect and even reduce the robot accuracy after calibration. This motivates careful planning of the calibration experiments in order to increase calibration efficiency, which is in the focus of the paper.

In the frame of this paper, a general planar manipulator with revolute joints is considered, whose the end-effector position is computed using the following expressions where l 0 i are the nominal link lengths, ∆ l i are their deviations, q 0 j are nominal joint coordinates, the variables θ 0 i are defined as ∑ i j=1 q 0 j , ∆ θ i are the joint offsets and n is the number of links. Collecting the unknown parameters (∆ l i and ∆ θ i ) into the vector ∆ π and the measurements (x k , y k ) into the vector ∆ P k , expression (1) can be rewritten as ∆ P = J ∆ π, where J is a Jacobian matrix. Then, one can get unknown parameters using the least-square technique that leads to

x = ∑ n i=1 l 0 i + ∆ l i • cos(θ 0 i + ∆ θ i ); y = ∑ n i=1 l 0 i + ∆ l i • sin(θ 0 i + ∆ θ i ) (1)
∆ π = ∑ m k=1 J k T J k -1 • ∑ m k=1 J k T ∆ P k (2)
where the superscript k indicates the experiment number, and m is the number of measurements.

Taking into account that each measurement is corrupted by the unbiased random Gaussian noise with standard deviation (std) σ , the identification accuracy of the parameters ∆ π can be evaluated via the covariance matrix, which is computed as follows:

cov(∆ π) = σ 2 ∑ m k=1 J k T J k -1 (3)
The latter allows us to evaluate the impact of measurement noise on the parameters identification accuracy and can be used to estimate the quality of measurement configurations. Based on the covariance matrix, it is possible to choose the measurement configurations that yield parameters less sensitive to measurement noise. In engineering practice, this procedure is referred to as the design of calibration experiments. In order to compare different plans of experiments (i.e. their efficiency), different performance measures have been proposed that deal with the norm, trace, etc. of the covariance matrix or information matrix (its inverse) [START_REF] Atkinson | Optimum experimental designs, with SAS[END_REF]. In [START_REF] Klimchik | Design of experiments for calibration of planar anthropomorphic manipulators[END_REF] the authors proposed an optimality condition of the calibration plan where the main idea is to ensure that the information matrix is diagonal. This condition corresponds to the D-optimal plan of experiments and is satisfied when

∑ m k=1 cos ∑ i s= j q k s = 0, ∑ m k=1 sin ∑ i s= j q k s = 0, ∀ i, j = 2, n, i ≥ j (4)
i.e. the sums of sines and cosines of links orientation q k s with respect to any joint s should be equal to zero.

It has been proved that the above presented equations define the desired set of the optimal measurement configurations, but in general case solution of this system requires essential efforts. In fact, for n = 2, m ≥ 2 expressions (4) provide two simple constraints only ∑ m k=1 cos q k 2 = 0 , ∑ m k=1 sin q k 2 = 0, which are imposed on the 2m design variables q k 1 , q k 2 | k = 1, m . So, relevant configurations can be easily found geometrically [START_REF] Klimchik | Design of experiments for calibration of planar anthropomorphic manipulators[END_REF]. However, even in the case of n = 3, m ≥ 3 there are 6 rather complicated trigonometric constraints imposed on 3m design variables q k 1 , q k 2 , q k 3 | k = 1, m . Therefore, the solution of the corresponding underconstrained system of algebraic equations becomes non-trivial. It is apparent that the complexity of the problem essentially increases with n (the number of manipulator links). Additional difficulties arise when the joint limits should be taken into account (i.e.

q min i ≤ q k i ≤ q max i | i = 1, n ; k = 1, m ).
Hence, it looks reasonable to develop relevantly simple and efficient technique allowing practical engineers to select desired measurement configurations in accordance with equations (4) without any tedious computations. This problem is in the focus of this paper.

Geometrical patterns for measurement pose selection

Before defining patterns, let us obtain some important properties of the optimality conditions (4) that allow us to reduce the problem complexity.

Property 1: Superposition of optimal plans gives also an optimal plan for this.

Proof of this property is obvious and is based on the additivity of the operations included in [START_REF] Hollerbach | Model identification[END_REF]. Using this property it is possible to generate optimal plans with a large number of measurement configurations using simple sub-sets. In a trivial case, it is possible to repeat experiments using the same set of optimal configurations. The latter is very important in practice, since it is much less time consuming to repeat experiments without changing manipulator configurations.

Property 2: The angles q 2 , ...q n can be rearranged in the optimal plans in an arbitrary manner without losing the optimality conditions (4).

Proof. Let us assume that for an n-link manipulator optimal plan consists of m measurement configurations. Let us now change in the optimal plan column a and b (a, b > 1) and analyze the sets of angles for which sums of sines and cosines should be satisfied in the optimal plan ∑ i s= j q k s . In this case conditions (4) will be satisfied for i, j = 2, a and i, j = b, n since all arguments of sines and cosines remain the same. For case a ≤ j ≤ i ≤ b conditions also do not change. However, in the cases a ≤ j ≤ b ≤ i and j ≤ a ≤ i ≤ b optimality conditions ∑ i s= j q k s will be replaced ∑ i s= j q k s = ∑ i s= j q k s + q k aq k b by and ∑ i s= j q k s = ∑ i s= j q k s + q k bq k a respectively, where prime indicates the angle corresponding to the rearranged new set of angles. Assuming that plan ∑ i s= j q k s is optimal, plans ∑ i s= j q k sq k a + q k b and ∑ i s= j q k s + q k aq k b will be also optimal. Further, taking into account that the superposition of optimal plans gives optimal plan and applying it to the last two plans gives plan 2 ∑ i s= j q k s that should be also optimal. This plan can be split into two sub-plans ∑ i s= j q k s for which the sums of sines and cosines should be equal. Taking into account that double sum is equal to zero, the plan ∑ i s= j q k s will be also optimal. This property can be useful when the joint limits are narrower than the minimum range required for the optimal plan. In this case, it is reasonable to change joint coordinates with the ones where joint limits are not so critical (more than 240 • ) and where the required range of joint variations is lower.

Property 3: Optimal plan for n-link manipulator can be obtained using two lower-order optimal plans for n 1 -and n 2 -link manipulators, where n 1 + n 2 = n + 1.

This property gives us an elegant technique to generate optimal plans of calibration experiments without straightforward solution of system (4). In this case the number of measurement configurations is defined by n 1 × n 2 . In fact, by a sequential splitting of the original n-link kinematic chain, it is possible to reduce the problem complexity and to replace the problem of generating optimal plan for n-link manipulator by several sub-problems for 2-, 3-and 4-link manipulators. The latter motivates us to develop typical geometrical patterns that can be used directly to generate optimal plans for complex manipulators. Now let us introduce some geometrical patterns for typical serial manipulators that can be used to generate optimal plans. In the frame of these patterns, all variables α i , β i and γ, δ are treated as arbitrary angles.

For n = 2, m = 2, the complete set of solutions can be expressed as:

q 1 1 = α 1 ; q 1 2 = β q 2 1 = α 2 ; q 2 2 = β + π (5)
For n = 2, m = 3, the desired set of solutions is

q 1 1 = α 1 ; q 1 2 = β q 2 1 = α 2 ; q 2 2 = β + 2π/3 q 3 1 = α 3 ; q 3 2 = β -2π/3 (6)
For n = 3, m = 3, the geometrical pattern can be presented as

q 1 1 = α 1 ; q 1 2 = β ; q 1 3 = γ q 2 1 = α 2 ; q 2 2 = β + 2π/3; q 2 3 = γ + 2π/3 q 3 1 = α 3 ; q 3 2 = β -2π/3; q 3 3 = γ -2π/3 (7)
For n = 3, m = 4, an optimal solution can be expressed as According to the geometrical patterns presented above, the optimal plans for the first and second sub-chains can be generated using Eqs. ( 7) and (9), respectively. For the entire manipulator, the optimal plan can be obtained using Property 3 (it consists of 12 measurement configurations). However, such a plan of experiments requires a large range of variations for q 2 (equal to 2π/3), which cannot be included in the joint limits. This difficulty can be overcome using Property 2 allowing us permutation of q 2 and q 4 . This provides us with the following plan of calibration experiments, which satisfies Eq. ( 4) and, consequently, insures the covariance matrix diagonality:

α 1 ; δ 1 ; γ β ; χ; ϕ α 1 ; δ 1 + π; γ β ; χ; ϕ + π α 1 ; δ 2 ; γ β ; χ + π 2 ; ϕ + δ 1 -δ 2 α 1 ; δ 2 + π; γ β ; χ + π 2 ; ϕ + δ 1 -δ 2 + π α 2 ; δ 1 ; γ + 2π/3 β + 2π/3; χ; ϕ α 2 ; δ 1 + π; γ + 2π/3 β + 2π/3; χ; ϕ + π α 2 ; δ 2 ; γ + 2π/3 β + 2π/3; χ + π 2 ; ϕ + δ 1 -δ 2 α 2 ; δ 2 + π; γ + 2π/3 β + 2π/3; χ + π 2 ; ϕ + δ 1 -δ 2 + π α 3 ; δ 1 ; γ -2π/3 β -2π/3; χ; ϕ α 3 ; δ 1 + π; γ -2π/3 β -2π/3; χ; ϕ + π α 3 ; δ 2 ; γ -2π/3 β -2π/3; χ + π 2 ; ϕ + δ 1 -δ 2 α 3 ; δ 2 + π; γ -2π/3 β -2π/3; χ + π 2 ; ϕ + δ 1 -δ 2 + π (10)
Taking into account given joint limits, the arbitrary angles in the above expression have been set as

α 1 = 0; α 2 = -π/3; α 3 = -2π/3; β = 0; δ 1 = -2π/3; δ 2 = -π/2; γ = 0; χ = π/6; ϕ = π/2 (11)
Using this plan of experiments, the desired geometrical parameters can be identified with a precision of 0.029 mm for the link lengths and 4.13 mdeg for the joint offsets (assuming that the measurement system precision is equal to 0.1 mm). In more details, relevant results are presented in Table 1. For comparison purposes, the identification accuracy has been also evaluated for regular and random plans of experiments that contain the same number of measurement configurations. It should be noted that for a 6-dof manipulator it is not possible to obtain a regular plan with 12 measurement configurations. For this reason the plan has been generated using 5 measurement configurations within the joint limits from which 12 configurations have been selected randomly. In order to reduce the factor of particular set of measurement configurations, simulations for regular and random plans have been repeated 1000 times and the results have been averaged.

Hence, simulation results confirm advantages of the proposed approach. In this study, the identification accuracy for the optimal plan is better by 34 -72% comparing with the random plan and by 38 -48% compared with the regular plan of experiments. In addition, maximum positioning errors have been reduced by 44% and 57%, respectively. It should be stressed that the proposed approach does not require any computation for optimal measurement pose selection and is able to improve the identification accuracy using a small number of measurement configurations.

Conclusions

The paper presented a new approach for the design of calibration experiments for robotic manipulators that essentially simplifies the optimal pose selection procedure. The main theoretical results are expressed as a set of several properties and rules, which allow user to obtain optimal measurement configurations without any computation, simply using superpositions and permutations of the proposed geometrical patterns describing optimal measurement configurations for 2-, 3-and 4link manipulators. The efficiency of the developed approach has been confirmed by an illustrative example that deals with the calibration of a 6-dof manipulator.

Fig. 1

 1 Fig. 1 Manipulator accuracy after calibration for two different sets of measurement poses.

Table 1

 1 Identification accuracy for different plans of calibration experiments. , [mdeg]δ L i , [mm] δ q i , [mdeg] δ L i , [mm] δ q i , [mdeg] δ L i , [mm]

	Joint	Proposed plan	Regular plan	Random plan
	δ q i 1	1.18	0.029	1.70	0.041	1.90	0.046
	2	1.50	0.029	2.21	0.042	2.58	0.050
	3	1.65	0.029	2.32	0.041	2.46	0.043
	4	2.07	0.029	3.04	0.042	3.32	0.046
	5	2.76	0.029	4.04	0.042	4.39	0.046
	6	4.13	0.029	5.73	0.040	5.55	0.039
	ρ max	0.61		0.88		0.96

Acknowledgements

The work presented in this paper was partially funded by the project ANR (Project ANR-2010-SEGI-003-02-COROUSSO), France and FEDER ROBOTEX, France.

q 1 1 = α 1 ; q 1 2 = β 1 ; q 1 3 = γ q 2 1 = α 2 ; q 2 2 = β 1 + π; q 2 3 = γ q 3 1 = α 3 ; q 3 2 = β 2 ; q 3 3 = γ + π q 4 1 = α 4 ; q 4 2 = β 2 + π;

For n = 4, m = 4, the desired geometrical pattern can be written as

When the optimal plan for the entire n-link manipulator is defined, it is required to fix all arbitrary variables. In case there is no joint limit, they can be set to any value. In contrast, in case the joint angles are limited, these variables should be properly selected using the following rules:

Rule 1. If the joint limits for the angle are narrower than 2π/3 while plan of experiments requires higher width for joint variations, the joint coordinates of this joint should be changed with the one for which the joint coordinates can be set within the interval of the length π.

Rule 2. If the joint limits for the angle are narrower than π, some optimality conditions in (4) cannot be satisfied and the plan of experiments should be modified. In this case, the best results can be achieved when the joint coordinates are set to the joint limits.

Rule 3. For the joint with a range of variation equal to π, an arbitrary parameter can be fixed in the interval [q min ; q maxπ], where q min and q max are the corresponding joint limits.

Rule 4. For the joint with a range of variation equal to 2π/3, an arbitrary parameter can be fixed in the interval [q min + π/3; q maxπ/3].

Hence, the above defined properties and rules allow us to generate optimal plan of experiments for the complex manipulator using simple patterns without any tedious computation. In the next section the obtained results will be applied for the generation of a plan of experiments for a 6-dof manipulator.

Case study: optimal measurement poses for a 6-dof robot

To show the efficiency of the developed optimal measurement pose selection technique, let us consider a 6-dof serial manipulator with six revolute joints and six links (l 1 = 1.4 m, l 2 = 1.1 m, l 3 = 1.0 m, l 4 = 0.8 m, l 5 = 0.6 m, l 6 = 0.4 m), whose joint limits are -140

In order to simplify the determination of measurement configurations, the manipulator can be split into two sub-chains: (i) 3-link chain corresponding to links #1 -#3, and (ii) 4-link chain composed of a virtual link/joint and links #4 -#6 of entire manipulator.