Fadoua Balabdaoui

Cécile Durot
email: cecile.durot@gmail.com

Babagnidé François Koladjo

Testing convexity of a discrete distribution

published or not. The documents may come

The testing problem

Modeling count data is an important issue in statistical research, see e.g. [START_REF] Gómez-Déniz | The discrete lindley distribution: properties and applications[END_REF]. A popular parametric model for such data is the Poisson model.

While non-parametric extensions are conceivable, those that only assume a shape constraint of the underlying probability mass function (pmf) may offer more flexibility. Recent papers on estimating a pmf under a shape constraint are [START_REF] Jankowski | Estimation of a discrete monotone distribution[END_REF], [START_REF] Durot | Least-squares estimation of a convex discrete distribution[END_REF][START_REF] Durot | Nonparametric species richness estimation under convexity constraint[END_REF], [START_REF] Balabdaoui | Asymptotic distribution of the discrete log-concave mle and some applications[END_REF], [START_REF] Giguelay | Estimation of a discrete probability under constraint of k-monotony[END_REF], [START_REF] Chee | Nonparametric estimation of species richness using discrete k-monotone distributions[END_REF]. In any case, it is sensible to validate the chosen model using a goodness-of-fit test. Goodnessof-fit tests to validate a model connected to the Poisson distribution are given in [START_REF] Karlis | A simulation comparison of several procedures for testing the poisson assumption[END_REF], [START_REF] Meintanis | A class of count models and a new consistent test for the poisson distribution[END_REF], [START_REF] Ledwina | On charlier polynomials in testing poissonity[END_REF] and the references therein. However, to the best of our knowledge, no goodness-of-fit test has been proposed to validate an assumed shape constraint for discrete data.

Motivated by the biological application in [START_REF] Durot | Nonparametric species richness estimation under convexity constraint[END_REF], where the number of species is estimated assuming a convex abundance distribution, we develop here a goodnessof-fit test for convexity of the underlying pmf on N. To the best of our knowledge, this is the first attempt to build a convexity test for count data. Precisely, based on i.i.d. observations X 1 , . . . , X n from some pmf p 0 on N, we test the null hypothesis H 0 : "p 0 is convex on N" (i.e. p 0 (k + 1) -p 0 (k) ≥ p 0 (k) -p 0 (k -1) for all integers k ≥ 1) versus H 1 : "p 0 is not convex".

The test is based on the convex least-squares estimator p n := argmin p∈C 1 p n -p , where C 1 is the set of all convex pmfs on N, q 2 = j∈N (q(j)) 2 for any sequence q = {q(j), j ∈ N}, and p n (j) = n -1 n i=1 1 {X i =j} , j ∈ N, is the empirical pmf. It is proved in Durot et al. (2013, Sections 2.1 to 2.3) that p n exists, is unique, and can be implemented with an appropriate algorithm. We reject H 0 if {T n > t α,n } where T n = √ n p n -p n and t α,n is an appropriate quantile, chosen in such a way that the test has asymptotic level α.

In the sequel, we assume that p 0 has a finite support in {0, . . . , S} with an unknown integer S > 0 and we consider two different constructions of t α,n . First, we define t α,n as the (1-α)-quantile of a random variable whose limiting distribution coincides with the limiting distribution of T n under H 0 . Next, we calibrate the test under a least favorable hypothesis (when the true pmf is triangular). Theoretical justification requires knowledge of the limiting distribution of T n under H 0 . This needs some notation. For all p = {p(j), j ∈ N}

and k ∈ N\{0} we set ∆p(k) = p(k + 1) -2p(k) + p(k -1) (hence p is convex on N iff ∆p(k) ≥ 0 for all k) and a given k ∈ N\{0} is called a knot of p if ∆p(k) > 0.
For all s > 0 and u = (u(0), . . . , u(s + 1)) ∈ R s+2 , we set u 2 s = s+1 k=0 (u(k)) 2 . Also, let g 0 be a (S + 2) centered Gaussian vector whose dispersion matrix Γ 0 has component (i + 1, j + 1) equal to 1 {i=j} p 0 (i) -p 0 (i)p 0 (j) for all i, j = 0, . . . , S + 1, and g 0 the minimizer of g -g 0 S over g ∈ K 0 := g = (g(0), . . . , g(S + 1)) ∈ R S+2 : ∆g(k) ≥ 0 for all k ∈ {1, . . . , S} such that ∆p 0 (k) = 0 .

Existence, uniqueness and characterization of g 0 are given in Balabdaoui et al. (2017, Theorem 3.1). The asymptotic distribution of T n under H 0 is given below.

Theorem 1.1. (i)The distribution function of T 0 := g 0 -g 0 S is continuous on (0, ∞).

(ii) Under H 0 , T n d -→ T 0 and sup t≥0 |P (T n ≤ t) -P ( T 0 ≤ t)| → 0, as n → ∞.

Calibrating by estimating the limiting distribution

Here, we build a random variable that weakly converges to T 0 and can be approximated via Monte-Carlo simulations. Let S n = max{X 1 , . . . , X n }, and let g n be a random vector which, conditionally on (X 1 , . . . , X n ), is a S n +2 centered Gaussian vector whose dispersion matrix Γ n has component (i + 1, j + 1) equal to 1 {i=j} p n (i) -p n (i)p n (j) for all i, j = 0, . . . , S n + 1. Now, let g n = argmin g∈Kn g -g n Sn , the least squares projection of g n on K n , where K n "approaches" K 0 as n → ∞:

K n = g = (g(0), . . . , g(S n + 1)) ∈ R Sn+2 : ∆g(k) ≥ 0 for all k ∈ {1, . . . , S n } such that ∆ p n (k) ≤ v n
for an appropriate positive sequence (v n ) n . Choosing v n = 0 would make K n to be the largest possible and hence g n -g n Sn the smallest possible; this distance would be stochas-tically smaller than the actual limit of T n , yielding a large probability of rejection. In fact, choosing v n = 0 amounts to estimate the knots of p 0 by those of p n , which is not desirable since p n has typically more knots than p 0 . The conditions required on v n are given below.

Theorem 2.1. Let g n , K n , and g n be as above, and take v n > 0 such that v n = o(1) and v n n -1/2 . (i) Then, g n uniquely exists, both g n and T n := g n -g n Sn are measurable.

(ii) Under H 0 , conditionally on X 1 , . . . , X n we have T n d -→ T 0 in probability as n → ∞.

By (i) in Theorem 1.1, the conditional convergence T n d -→ T 0 in probability means that sup t∈R P ( T n ≤ t|X 1 , . . . , X n ) -P ( T 0 ≤ t) = o p (1). (2.1)
In Balabdaoui et al. (2017, Theorem 3.3) it is shown that lim n→∞ P (K n = K 0 ) = 0 for any (v n ) n satisfying the conditions of the theorem. The intuition behind is as follows: when k is a knot of p 0 and ∆ p n (k) ≤ v n , then √ n(∆ p n (k) -∆p 0 (k)) < -√ n 0 for some 0 > 0 and n large enough. Weak convergence of p n to p 0 implies that this happens with probability converging to zero. In case k is not a knot; i.e., ∆p

0 (k) = 0 such that ∆ p n (k) > v n then √ n∆ p n (k) > √ nv n → ∞
, which again happens with decreasing probability. We now state the main result of the section, which is proven is the supplement.

Theorem 2.2. Let T n as in Theorem 2.1. Let α ∈ (0, 1) and t α,n the conditional

(1 -α)- quantile of T n given X 1 , . . . , X n . Under H 0 , lim sup n→∞ P T n > t α,n ≤ α.
Hence, the test is asymptotically calibrated. In fact, it can be shown that the asymptotic Type I error is precisely equal to α for appropriate range of α, i.e. lim n→∞ P (T n > t α,n ) = α, see the supplementary material. An approximative value of t α,n can be computed using Monte-Carlo simulations as follows. Having observed X 1 , . . . , X n , draw independent sequences (Z

(b) i ) 0≤i≤Sn+1 for b ∈ {1, . . . , B}, of i.i.d. N (0, 1) variables Z (b)
i , where B > 0 is the number of Monte-Carlo runs. For all b, compute g

(b) n = Γ 1/2 n (Z (b) 0 , . . . , Z (b) Sn+1 ) T and g (b)
n the minimizer of g

(b)
n -g Sn over K n using Dykstra's algorithm, see [START_REF] Balabdaoui | On asymptotics of the discrete convex LSE of a p[END_REF]. If (s j ) j is the sequence of successive knots of p n such that ∆ p n (s j ) > v n , then the algorithm performs iterative projections on the cones {g ∈ R Sn+2 : ∆g(k) ≥ 0 for all k ∈ {s j + 1, . . . , s j+1 -1}}, whose intersection is precisely K n . See [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF] for more details and a proof of convergence. Then, t α,n can be approximated by the (1 -α)-quantile of the empirical distribution corresponding to g

(b) n -g (b) n
Sn , with b ∈ {1, . . . , B}.

Calibrating under the least favorable hypothesis

We consider below an alternative calibration that is easier to implement than the first one since it does not involve a sequence (v n ). In what follows we denote by T a the triangular pmf supported on {0, . . . , a -1} for a given integer a ≥ 1; i.e., T a (i) = 2(a -i) + [a(a + 1)] -1 .

Consider K 0 the set of all vectors g = (g(0), . . . , g(S+1)) ∈ R S+2 such that ∆g(k) ≥ 0 for all k ∈ {1, . . . , S}. Similarly, let K n be the set of all vectors g = (g(0), . . . , g(S n + 1)) ∈ R Sn+2 such that ∆g(k) ≥ 0 for all k ∈ {1, . . . , S n }. Let g0 be the least squares projection of g 0 onto K 0 and gn that of g n onto K n , with g 0 and g n as in Sections 1 and 2. Finally, let tα,n be the conditional (1 -α)-quantile of T n := gn -g n Sn given (X 1 , . . . , X n ). Recall that T 0 = g 0 -g 0 S is the weak limit of T n . We start with the following proposition.

Proposition 3.1. We have g0 -g 0 S ≥ T 0 . Moreover under H 0 , conditionally on

X 1 , . . . , X n , T n d -→ g0 -g 0 S in probability, as n → ∞.
Hence, T n is typically larger that T n , which provides the intuition for the following theorem. A proof of the theorem is given in the Supplement.

Theorem 3.2. Let α ∈ (0, 1). Under H 0 , lim sup n→∞ P T n > tα,n ≤ α with equality if p 0 is T S+1 and α < α0 , with α0 = P ( g0

-g 0 S = 0) ∈ [1/2, 1).
The test is asymptotically calibrated since the Type I error does not exceed α. It reaches precisely α when p 0 is triangular, which can be viewed as the least favorable hypothesis for testing convexity. The test is thus refereed to as the LFH test in the sequel. The theorem above does not exclude existence of other least favorable hypotheses. Since with an appropriate sequence (v n ), the test of Section 2 has the exact asymptotic Type I error α, it is typically less conservative and more powerful than the LFH test. The performance could be however poor with a bad choice of (v n ). On the other hand, the LFH test has the great advantage of not requiring the choice of such a sequence, which makes it easier to implement.

Simulations

To illustrate the theory, we have considered four pmf's supported on {0, . . . , 5}. It follows from Theorem 7 in [START_REF] Durot | Least-squares estimation of a convex discrete distribution[END_REF] that any convex pmf on N can be written as k≥1 π k T k where π k ∈ [0, 1], k≥1 π k = 1 and T k is the triangular pmf supported on {0, . . . , k -1} already defined above. Under H 0 , we considered the triangular pmf p

(1) 0 = T 6 and p

(2) 0 = 6 k=1 π k T k with π 1 = 0, π 2 = π 3 = 1/6, π 4 = 0 and π 5 = π 6 = 1/3, which has knots at 2, 3, 5 and 6. Under H 1 we considered p (1) 0 (1) -0.008). To investigate the asymptotic type I error and power of our tests, we have drawn n ∈ {50, 100, 500} rv's from the aforementioned pmf's. Here, α = 5%, t α,n and tα,n were estimated for each drawn sample using B = 1000 i.i.d. copies of g n . The rejection probability was estimated using N = 500 replications of the whole procedure. For the first convexity test, we considered the sequences v n = log(log n)n -1/2 and n -1/4 . We also added the sequence v n = 0 to compare our approach PMF n = 50 n = 100 n = 500 0 (2) 0 and power for p

(log(log n)) 1/2 n 1/2 n -1/4 0 (log(log n)) 1/2 n 1/2 n -1/4 0 (log(log n)) 1/2 n 1/2 n -1/4 p (
(1) 1 and p

(2)

1 of the test based on the least favorable hypothesis. The asymptotic level is 5%.

with the naive one where no knot extraction is attempted. The results are reported in Tables 1 and2. Some more results for pmfs with a larger support are given in the Supplement.

The simulations confirm that the first test depends on the choice of v n and that choosing v n = 0 does not yield a valid test. The choice v n = log(log n)n -1/2 is appropriate for p

(2) 0 whereas v n = n -1/4 is appropriate for p

(1) 0 . Small sequences can have the same effect as choosing v n = 0. On the contrary, large sequences can yield a conservative test in case the true pmf has small slope changes so that the estimated conditional quantile is too large to be exceeded: in such a case it may happen that no knot is extracted up so that the distance g n -g n Sn is too large. The question is then open as to how to choose such a sequence.

To select among possible candidates of sequences (v n ) n , it would be possible to resort to a cross-validation approach. This is however beyond the main scope of this article.

The second testing approach is, as expected, conservative when the true pmf is not triangular. For p (1) 1 , which strongly violates the convexity constraint, the power is equal to 1 for n = 500. For p

(2)
1 , which has only a small flaw at 2 with a change of slope equal to -0.008, the power is comparable to that obtained with the first testing method and v n ≡ n -1/4 . This is somehow expected as it is the largest sequence among the ones considered, yielding the largest distance g n -g n Sn . As for any result based on asymptotic theory, it is difficult to know in advance the range of the sample size which would guarantee the applicability of this theory, since this would require a theoretical statement on the convergence of the involved processes at the second order. Our numerical investigation above shows that n = 500 is enough to obtain the correct type I error for both tests.

Proofs

We begin with two preparatory lemmas and then go to the proof of the main results.

Lemma 5.1. Let s ∈ N\{0} and K ⊂ R s+2 a non-empty closed convex set. Then, for u = (u(0), . . . , u(s + 1)) ∈ R s+2 , the minimizer of g -u s over g ∈ K uniquely exists.

Denoting by Φ(u) the minimizer, the applications Φ and u → u -Φ(u) s are measurable.

Proof: It follows from standard results on convex optimization that Φ(u) uniquely exists for all u, and Φ(u)-Φ(v) s ≤ u-v s for all u and v in R s+2 . Hence, Φ is continuous, whence it is measurable. Now, the function u → (u, Φ(u)) is continuous, whence measurable. By continuity of the norm, this ensures that u → u -Φ(u) s is measurable.

Lemma 5.2. If p 0 has support {0, . . . , S}, with g 0 = (g 0 (0), . . . , g 0 (S + 1)) as in Section 1, (∆g 0 (1), . . . , ∆g 0 (S)) is a centered Gaussian vector with invertible dispersion matrix.

Proof: For notational convenience, we assume in the sequel that S ≥ 3. The case S ≤ 2 can be handled likewise. Let B be the S × (S + 1)-matrix which j-th line has components j,j + 1 and j + 2 equal to 1, -2 and 1 respectively while the other components are zero, for j = 1, . . . , S -1, and S-th line has components equal to zero except the penultimate and the last one, which are equal respectively to 1 and -2. We have p n (S + 1) = p 0 (S + 1) = 0 almost surely so that in the limit, g 0 (S + 1) = 0 almost surely and (∆g 0 (1), . . . , ∆g 0 (S)) T = B (g 0 (0), . . . , g 0 (S)) T . Hence, (∆g 0 (1), . . . , ∆g 0 (S)) is a centered Gaussian vector with dispersion matrix V = BΣ 0 B T where Σ 0 is the dispersion matrix of (g 0 (0), . . . , g 0 (S)) T , i.e. Σ 0 is obtained by deleting a line and a column of zeros in Γ 0 , the dispersion matrix of g 0 . It remains to prove that V is invertible.

Let

√ p 0 be the column vector with components p 0 (0), . . . , p 0 (S) and let diag( √ p 0 ) be the (S + 1) × (S + 1) diagonal matrix with diagonal vector √ p 0 . Denoting by I the identity matrix on R S+1 , the matrix (in the canonical basis) associated with the orthogonal projection from R S+1 onto the orthogonal supplement of the linear space generated by √ p 0 is given by I -Π 0 = I -√ p 0 √ p T 0 . The linear subspace of R S+1 generated by √ p 0 has dimension 1, so its orthogonal supplement in R S+1 has dimension S, whence rank(I -

Π 0 ) = S. Now, Σ 0 = diag( √ p 0 )(I -Π 0 )diag( √ p 0 ) where diag( √ p 0 ) is invertible, so that rank(Σ 1/2 0 ) = rank(Σ 0 ) = S.
Hence, the kernel of Σ 1/2 0 is a linear subspace of R S+1 of dimension 1. Let λ be the column vector in R S+1 whose components are all equal to 1.

Using that S k=0 p 0 (k) = 1, it is easy to see that Σ 0 λ is the null vector in R S+1 , whence λ T Σ 0 λ = 0. This means that Σ 1/2 0 λ 2 = 0 with . the Euclidean norm in R S+1 .

Hence, Σ 1/2 0 λ is the null vector in R S+1 , so the kernel of Σ 1/2 0 is the linear subspace of R S+1 generated by λ.

Next, let us determine the kernel of V . Let µ = (µ 1 , . . . , µ S ) ∈ R S with V µ = 0. Then,

µ T V µ = 0 where V = BΣ 0 B T , so that Σ 1/2 0 B T µ 2 = µ T BΣ 0 B T µ = 0. Hence, Σ 1/2 0 B T µ = 0. Since the kernel of Σ 1/2 0
is the linear subspace of R S+1 spanned by λ, we conclude that

B T µ = aλ for some a ∈ R. Hence, µ 1 = a, µ 2 -2µ 1 = a, µ k-2 -2µ k-1 + µ k = a for all k ∈
and v ∈ A t . By definition, both u and v belong to the convex set K 0 and therefore,

λ u + (1 -λ) v ∈ K 0 . Since x minimizes g -x S over g ∈ K 0 , we conclude that x -x S ≤ λ u + (1 -λ) v -x S ≤ λ( u -u) + (1 -λ)( v -v) S , using that x = λu + (1 -λ)v. It then follows from the triangle inequality that x -x S ≤ λ u -u S + (1 -λ) v -v S . Since u ∈ A t and v ∈ A t , we have u -u S ≤ t and v -v S ≤ t , which implies that x -x S ≤ λt + (1 -λ)t . Hence, x ∈ A λt+(1-λ)t and λA t + (1 -λ)A t ⊂ A λt+(1-λ)t .
(5.3)

Now, µ = P • g -1 0 is a Gaussian probability measure on B(R S+2 ), so it follows from Lemma 1.1 in [START_REF] Gaenssler | On continuity and strict increase of the cdf for the sup-functional of a gaussian process with applications to statistics[END_REF] that µ is log-concave in the sense that µ (λA + (1 -λ)B) ≥ µ(A) λ µ(B) 1-λ for all λ ∈ (0, 1) and A, B ∈ B(R S+2 ), with µ the inner measure pertaining to µ. Applying this with A = A t and B = A t , and combining with (5.3) yields µ (A λt+(1-λ)t ) ≥ µ(A t ) λ µ(A t ) 1-λ for all t, t ∈ [0, ∞) and λ ∈ (0, 1). The same inequality remains true with µ replaced by µ. Using the equality F (t) = µ(A t ), t ∈ [0, ∞), and taking the logarithm on both sides of the inequality, we conclude that log(F (λt

+ (1 -λ)t )) ≥ λ log(F (t)) + (1 -λ) log(F (t )
) for all λ ∈ (0, 1), and t, t ∈ [0, ∞). This means that the function log(F ) is concave on [0, ∞). Recalling (5.2), this implies (i).

Now, with probability one, p n is supported on {0, . . . , S + 1} for sufficiently large n (see Balabdaoui et al. (2017, Proposition 3.5)), and p n also is supported on that set by definition, so T n = √ n p n -p n S . It can be proved in the same manner as (Balabdaoui et al., 2017, Theorem 3.2) that under H 0 , √ n ( p n -p 0 , p n -p 0 ) d -→ ( g 0 , g 0 ) as n → ∞ as a joint weak convergence on {0, . . . , S + 1}, so T n d -→ T 0 as n → ∞. The last assertion follows from the latter weak convergence together with (i), since T n ≥ 0 and T 0 ≥ 0.

Proof of Theorem 2.1: Clearly, K n is a non-empty closed convex subset of R Sn+2 .

Hence, we consider the specific case of s = S n and K = K n in Lemma 5.1. In the notation of the lemma, we have g n := Φ(g n ), so that g n is uniquely defined. Moreover, since both Φ and g n are measurable, we conclude that g n = Φ(g n ) is measurable. Likewise, g n -g n Sn is measurable. This proves the fist two assertions in Theorem 2.1. Next, similar to Balabdaoui et al. (2017, Theorem 3.3), the following joint weak convergence on {0, . . . , S + 1} can be proved under H 0 : conditionally on X 1 , . . . , X n , ( g n , g n ) d -→ ( g 0 , g 0 ) in probability as n → ∞. The result follows, since S n = S with provability that tends to one.

Proof of Proposition 3.1: The first assertion follows from the inclusion K 0 ⊂ K 0 . Next, similar to Balabdaoui et al. (2017, Theorem 3.3), the following joint weak convergence on {0, . . . , S + 1} can be proved under H 0 : conditionally on X 1 , . . . , X n , (g n , g n ) d -→ (g 0 , g 0 ) in probability as n → ∞. The second assertion follows since lim n→∞ P (S n = S) = 1.

  of a truncated Poisson on {0, . . . , 5} with rate λ = 1.5 and p

{3, . . . , S} and µ S-1 -2µ S = a. This is equivalent to 2µ k = ak(k + 1) for all k ∈ {1, . . . , S} and 2µ S = µ S-1 -a. Combining the first equation with k = S, S -1 to the second equation yields aS(S + 1) = -a + a(S -1)S/2. Therefore, a(2 + S 2 + 3S)/2 = 0, whence a = 0.

Hence, µ is the null vector in R S and therefore, rank(V ) = S, and hence V is invertible.

Proof of Theorem 1.1: Let F be the c.d.f. of T 0 := g 0 -g 0 S . Note that T 0 is a properly defined random variable by the measurability proved in Lemma 5.1. The proof of Theorem 1.1 is inspired by that of Lemma 1.2 in [START_REF] Gaenssler | On continuity and strict increase of the cdf for the sup-functional of a gaussian process with applications to statistics[END_REF]. First, we will prove that F (t) > 0 for all t ≥ 0.

(5.2)

To this end, note that F (0) = P (g 0 ∈ K 0 ), hence F (0) ≥ P (δ ∈ A) with A the set of all vectors (u 1 , . . . , u S ) ∈ R S such that u k ≥ 0 for all k ∈ {1, . . . , S}, and δ = (∆g 0 (1), . . . , ∆g 0 (S)). From Lemma 5.2, the dispersion matrix of the centered Gaussian vector δ is invertible, so δ has a strictly positive density with respect to the Lebesgue measure on R S . This implies that the probability that δ belongs to a Borel set whose Lebesgue measure is not equal to zero, is strictly positive, whence P (δ ∈ A) > 0. Thus,

Combining this with the monotonicity of F yields (5.2).

Next, we prove that the function log(F ) (which is well defined on [0, ∞) thanks to (5.2))

is concave on [0, ∞). For u = (u(0), . . . , u(S + 1)) ∈ R S+2 , let us write u the minimizer of g -u S over g ∈ K 0 . For all t ∈ [0, ∞), letA t be the set of all u = (u(0), . . . , u(S + 1)) ∈ R S+2 such that u-u S ≤ t. Note that A t is a Borel set in B(R S+2 ) for all t since, according to Lemma 5.1, the application u → u -u S is measurable. Finally, let µ = P • g -1 0 be the distribution of g 0 on R S+2 endowed with the Borel σ-algebra B(R S+2 ). This means that F (t) = µ(A t ). Fix λ ∈ (0, 1), t, t ∈ [0, ∞), and consider an arbitrary x ∈ λA t + (1 -λ)A t .

Then, x takes the form x = λu + (1 -λ)v for some (non necessarily unique) u ∈ A t