Jean-Luc Marichal

Ivan Kojadinovic

Distribution functions of linear combinations of lattice polynomials from the uniform distribution

Keywords: Lovász extension, discrete Choquet integral, lattice polynomial, order statistic, distribution function, moment, B-Spline, divided difference

We give the distribution functions, the expected values, and the moments of linear combinations of lattice polynomials from the uniform distribution. Linear combinations of lattice polynomials, which include weighted sums, linear combinations of order statistics, and lattice polynomials, are actually those continuous functions that reduce to linear functions on each simplex of the standard triangulation of the unit cube. They are mainly used in aggregation theory, combinatorial optimization, and game theory, where they are known as discrete Choquet integrals and Lovász extensions.

Introduction

Let h : [0, 1] n → R be an aggregation function and let X be a random vector uniformly distributed on [0, 1] n . An interesting but generally difficult problem Email addresses: jean-luc.marichal [at]uni.lu (Jean-Luc Marichal), ivan[at]stat.auckland.ac.nz (Ivan Kojadinovic). 1 Formerly at LINA CNRS FRE 2729, École Polytechnique de l'Université de Nantes, France.

Preprint submitted to Elsevier Science Submitted 18 Jul 06, Accepted 15 Sep 07

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT is to provide explicit expressions for the distribution function and the moments of the aggregated random variable Y = h(X).

This problem has been completely solved for certain aggregation functions (see for instance [21, §7.2]), especially piecewise linear functions such as weighted sums [START_REF] Barrow | Spline notation applied to a volume problem[END_REF] (see also [START_REF] Marichal | Slices, slabs, and sections of the unit hypercube[END_REF]), linear combinations of order statistics [START_REF] Agarwal | Linear functions of uniform order statistics and B-splines[END_REF][START_REF] Matsunawa | The exact and approximate distributions of linear combinations of selected order statistics from a uniform distribution[END_REF][START_REF] Weisberg | The distribution of linear combinations of order statistics from the uniform distribution[END_REF] (see also [7, §6.5] for an overview), and lattice polynomials [START_REF] Marichal | Cumulative distribution functions and moments of lattice polynomials[END_REF], which are max-min combinations of the variables.

In this note we solve the case of linear combinations of lattice polynomials, which include the three above-mentioned cases. Actually, linear combinations of lattice polynomials are exactly those continuous functions that reduce to linear functions on each simplex of the standard triangulation of [0, 1] n . In particular, these functions are completely determined by their values at the 2 n vertices of [0, 1] n .

The concept of linear combination of lattice polynomials is known in combinatorial optimization and game theory as the Lovász extension [START_REF] Algaba | The Lovász extension of market games[END_REF][START_REF] Grabisch | Equivalent representations of set functions[END_REF][START_REF] Lovász | Submodular functions and convexity[END_REF][START_REF] Singer | Extensions of functions of 0-1 variables and applications to combinatorial optimization[END_REF] of a pseudo-Boolean function (recall that a pseudo-Boolean function is a realvalued function of 0-1 variables). When it is nondecreasing in each variable, it is known in the area of nonlinear aggregation and integration as the discrete Choquet integral [START_REF] Denneberg | Non-additive measure and integral[END_REF][START_REF] Grabisch | Fuzzy measures and integrals[END_REF][START_REF] Marichal | Aggregation of interacting criteria by means of the discrete Choquet integral[END_REF], which is an extension of the discrete Lebesgue integral (weighted mean) to non-additive measures. The equivalence between the Lovász extension and the Choquet integral is discussed in [START_REF] Marichal | Aggregation of interacting criteria by means of the discrete Choquet integral[END_REF]. This note is set out as follows. In Section 2 we elaborate on the definition of linear combinations of lattice polynomials and we show how to concisely represent them. In Section 3 we provide formulas for the distribution function and the moments of any linear combination of lattice polynomials from the uniform distribution. Finally, in Section 4 we provide an application of our results to aggregation theory. Throughout we will use the notation [n] := {1, . . . , n}. Also, for any subset A ⊆ [n], 1 A will denote the characteristic vector of A in {0, 1} n . Finally, for any function h : [0, 1] n → R, we define the set function v h : 2

[n] → R as v h (A) := h(1 A ) for all A ⊆ [n].

Linear combinations of lattice polynomials

In the present section we recall the definition of lattice polynomials and we show how an arbitrary combination of lattice polynomials can be represented.

Basically an n-place lattice polynomial p : [0, 1] n → [0, 1] is a function defined from any well-formed expression involving n real variables x 1 , . . . , x n linked by the lattice operations ∧ = min and ∨ = max in an arbitrary combination of parentheses (see e.g. Birkhoff [6,§II.2]). For instance,

p(x 1 , x 2 , x 3 ) = (x 1 ∧ x 2 ) ∨ x 3 is a 3-place lattice polynomial.
Consider the standard triangulation of [0, 1] n into the canonical simplices

S σ := {x ∈ [0, 1] n | x σ(1) • • • x σ(n) } (σ ∈ S n ), (1)
where S n is the set of all permutations on [n]. Clearly, any linear combination of n-place lattice polynomials

h(x) = m i=1 c i p i (x)
is a continuous function whose restriction to any canonical simplex is a linear function. According to Singer [23, §2], h is then the Lovász extension of the pseudo-Boolean function h| {0,1} n , that is, the continuous function defined on each canonical simplex S σ as the unique linear function that coincides with h| {0,1} n at the n + 1 vertices

ε σ i := 1 {σ(1),...,σ(i)} (i = 0, . . . , n) of S σ . It can be written as [23, §2] h(x) = n i=1 h σ i -h σ i-1 x σ(i) (x ∈ S σ ), (2)
where 1), . . . , σ(i)} for all i = 0, . . . , n. In particular, h σ 0 = 0.

h σ i := h ε σ i = v h {σ(
Conversely any continuous function h : [0, 1] n → R that reduces to a linear function on each canonical simplex is a linear combination of lattice polynomials:

h(x) = A⊆[n] m h (A) i∈A x i (x ∈ [0, 1] n ), (3)
where

m h : 2 [n] → R is the Möbius transform of v h , defined as m h (A) := B⊆A (-1) |A|-|B| v h (B).
Indeed, expression (3) reduces to a linear function on each canonical simplex and agrees with h(1 B ) at

1 B for each B ⊆ [n].
Eq. ( 2) thus provides a concise expression for linear combinations of lattice polynomials. We will use it in the next section to calculate their distribution functions and their moments.

Remark 1 As we have already mentioned, the class of linear combinations of lattice polynomials covers three interesting particular cases, namely: lattice polynomials, linear combinations of order statistics, and weighted sums. These are characterized as follows. Let h : [0, 1] n → R be a linear combination of lattice polynomials.

(1) The function h reduces to a lattice polynomial if and only if the set function v h is monotone, {0, 1}-valued, and such that v h ([n]) = 1.

(2) As the order statistics are exactly the symmetric lattice polynomials (see [START_REF] Marichal | On order invariant synthesizing functions[END_REF]), the function h reduces to a linear combination of order statistics if and only if the set function

v h is cardinality-based, that is, such that v h (A) = v h (A ) whenever |A| = |A |. (3

) The function h reduces to a weighted sum if and only if the set function

v h is additive, that is, v h (A) = i∈A v h ({i}).

Distribution functions and moments

Before yielding the main results, let us recall some basic material related to divided differences. See for instance [START_REF] Davis | Interpolation and approximation[END_REF][START_REF] Devore | Constructive approximation[END_REF][START_REF] Powell | Approximation theory and methods[END_REF] for further details.

Consider the plus (resp. minus) truncated power function x n + (resp. x n -), defined to be x n if x > 0 (resp. x < 0) and zero otherwise. Let A (n) be the set of n -1 times differentiable one-place functions g such that g (n-1) is absolutely continuous. The nth divided difference of a function g ∈ A (n) is the symmetric function of n + 1 arguments defined inductively by ∆[g : a 0 ] := g(a 0 ) and ∆[g : a 0 , . . . , a n ] :=

         ∆[g : a 1 , . . . , a n ] -∆[g : a 0 , . . . , a n-1 ] a n -a 0 , if a 0 = a n , ∂ ∂a 0 ∆[g : a 0 , . . . , a n-1 ], if a 0 = a n .
The Peano representation of the divided differences, which can be obtained by a Taylor expansion of g, is given by

∆[g : a 0 , . . . , a n ] = 1 n! R g (n) (t) M (t | a 0 , . . . , a n ) dt, (4)
where M (t | a 0 , . . . , a n ) is the B-spline of order n, with knots {a 0 , . . . , a n }, defined as

M (t | a 0 , . . . , a n ) := n ∆[(• -t) n-1 + : a 0 , . . . , a n ]. (5)
We also recall the Hermite-Genocchi formula: For any function g ∈ A (n) , we have

∆[g : a 0 , . . . , a n ] = S id g (n) a 0 + n i=1 (a i -a i-1 )x i dx, (6)
where S id is the simplex defined in (1) when σ is the identity permutation.

For distinct arguments a 0 , . . . , a n , we also have the following formula, which can be verified by induction,

∆[g : a 0 , . . . , a n ] = n i=0 g(a i ) j =i (a i -a j ) . ( 7
)
Now, consider a random vector X uniformly distributed on [0, 1] n and set Y h := h(X), where the function h : [0, 1] n → R is a linear combination of lattice polynomials as given in formula (2). We then have the following result.

Theorem 2 For any function g ∈ A (n) , we have

E[g (n) (Y h )] = σ∈S n ∆[g : h σ 0 , . . . , h σ n ]. (8)
Proof. Using (2), we simply have

E[g (n) (Y h )] = [0,1] n g (n) [h(x)] dx = σ∈S n S σ g (n) n i=1 h σ i -h σ i-1 x σ(i) dx.
Finally, after an elementary change of variables, we conclude by the Hermite-Genocchi formula (6). 2

Theorem 2 provides the expectation E[g (n) (Y h )] in terms of the divided differences of g with arguments h σ 0 , . . . , h σ n (σ ∈ S n ). An explicit formula can be obtained by [START_REF] David | Wiley Series in Probability and Statistics[END_REF] whenever the arguments are distinct for every σ ∈ S n .

Clearly, the special cases

g(x) = r! (n + r)! x n+r , r! (n + r)! [x -E(Y h )] n+r
, and e tx t n (9)

give, respectively, the raw moments, the central moments, and the momentgenerating function of Y h . As far as the raw moments are concerned, we have the following result.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Proposition 3 For any integer r 1, setting A 0 := [n], we have,

E[Y r h ] = 1 n+r r A 1 ⊆[n] A 2 ⊆A 1 ••• Ar ⊆A r-1 r i=1 1 |A i-1 | |A i | v h (A i ).
Proof. Let r 1. It can be shown [START_REF] Ali | Content of the frustum of a simplex[END_REF] that ∆[(•) n+r : a 0 , . . . , a n ] = r 0 ,...,rn 0

r 0 +•••+r n =r a r 0 0 • • • a r n n = 0 i 1 ••• ir n a i 1 • • • a ir .
Hence, from ( 8) and ( 9) it follows that

E[Y r h ] = r! (n + r)! 0 i 1 ••• i r n σ∈S n h σ i 1 • • • h σ i r = r! (n + r)! 0 i 1 ••• ir n m∈M n v h (m i 1 ) • • • v h (m ir ),
where M n is the set of the n! maximal chains of the lattice (2 [n] , ⊆), and where, for any m ∈ M n , m i is the unique element of m of cardinality i.

For any

B 1 ⊆ • • • ⊆ B r ⊆ [n], let M B 1 ,...,B r n
denote the subset of maximal chains of (2 [n] , ⊆) containing B 1 , . . . , B r . It is then easy to see that, for any fixed 0 i 1 • • • i r n, the following identity holds:

B 1 ⊆•••⊆B r ⊆[n] |B 1 |=i 1 ,...,|Br|=ir M B 1 ,...,B r n = M n
and the union is disjoint. Therefore, we have

E[Y r h ] = r! (n + r)! 0 i 1 ••• ir n B 1 ⊆•••⊆B r ⊆[n] |B 1 |=i 1 ,...,|Br|=ir m∈M B 1 ,...,Br n v h (B 1 ) • • • v h (B r ) = r! (n + r)! B 1 ⊆•••⊆Br⊆[n] |M B 1 ,...,Br n | r i=1 v h (B i ),
where

|M B 1 ,...,B r n | = |B 1 |! (|B 2 | -|B 1 |)! (|B 3 | -|B 2 |)! • • • (n -|B r |)!.
Finally, we get the result by setting A i := B r+1-i for all i = 1, . . . , r. 2

Proposition 3 provides an explicit expression for the rth raw moment of Y h as a sum of (r + 1) n terms. For instance, the first two moments are

E[Y h ] = 1 n + 1 A⊆[n] 1 n |A| v h (A), E[Y 2 h ] = 2 (n + 1)(n + 2) A 1 ⊆[n] 1 n |A 1 | v h (A 1 ) A 2 ⊆A 1 1 |A 1 | |A 2 | v h (A 2 ).
We now yield a formula for the distribution function

F h (y) := Pr[Y h y] of Y h .
Theorem 4 There holds

F h (y) = 1 - 1 n! σ∈S n ∆[(• -y) n + : h σ 0 , . . . , h σ n ]. (10)
Proof. We have

F h (y) = 1 -Pr[h(X) > y] = 1 -E (Y h -y) 0 + .
Then, using [START_REF] Davis | Interpolation and approximation[END_REF] with

g(x) = 1 n! (x -y) n + leads to the result. 2
It follows from [START_REF] Denneberg | Non-additive measure and integral[END_REF] that the distribution function of Y h is absolutely continuous and hence its probability density function is simply given by

f h (y) = 1 (n -1)! σ∈Sn ∆[(• -y) n-1 + : h σ 0 , . . . , h σ n ] (11)
or, using the B-spline notation [START_REF] Barrow | Spline notation applied to a volume problem[END_REF],

f h (y) = 1 n! σ∈Sn M (y | h σ 0 , . . . , h σ n ).
Remark 5 (i) It is easy to see that [START_REF] Denneberg | Non-additive measure and integral[END_REF] can be rewritten by means of the minus truncated power function as

F h (y) = 1 n! σ∈Sn ∆[(• -y) n -: h σ 0 , . . . , h σ n ].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

(ii) When the arguments h σ 0 , . . . , h σ n are distinct for every σ ∈ S n , then combining [START_REF] David | Wiley Series in Probability and Statistics[END_REF] with [START_REF] Denneberg | Non-additive measure and integral[END_REF]

immediately yields the following explicit expression

F h (y) = 1 - 1 n! σ∈Sn n i=0 (h σ i -y) n + j =i (h σ i -h σ j )
, or, using the minus truncated power function,

F h (y) = 1 n! σ∈S n n i=0 (h σ i -y) n - j =i (h σ i -h σ j )
.

(iii) The knowledge of f h (y) immediately gives an alternative proof of [START_REF] Davis | Interpolation and approximation[END_REF]. Indeed, using Peano's representation ( 4), we simply have

E[g (n) (Y h )] = R g (n) (y) f h (y) dy = 1 n! σ∈Sn R g (n) (y) M (y | h σ 0 , . . . , h σ n ) dy = σ∈S n ∆[g : h σ 0 , . . . , h σ n ].
(iv) The case of linear combinations of order statistics is of particular interest. In this case, each h σ i is independent of σ (see Remark 1), so that we can write h i := h σ i . The main formulas then reduce to (see for instance [START_REF] Adell | Error bounds in divided difference expansions. A probabilistic perspective[END_REF] and [START_REF] Agarwal | Linear functions of uniform order statistics and B-splines[END_REF])

E[g (n) (Y h )] = n! ∆[g : h 0 , . . . , h n ], F h (y) = ∆[(• -y) n -: h 0 , . . . , h n ], f h (y) = M (y | h 0 , . . . , h n ).
We also note that the Hermite-Genocchi formula [START_REF] Birkhoff | Lattice theory[END_REF] provides nice geometric interpretations of F h (y) and f h (y) in terms of volumes of slices and sections of canonical simplices (see also [START_REF] Ali | Content of the frustum of a simplex[END_REF] and [START_REF] Gerber | The volume cut off a simplex by a half-space[END_REF]).

Both functions F h (y) and f h (y) require the computation of divided differences of truncated power functions. On this issue, we recall a recurrence equation, due to de Boor [START_REF] De Boor | On calculating with B-splines[END_REF] and rediscovered independently by Varsi [START_REF] Varsi | The multidimensional content of the frustum of the simplex[END_REF] (see also [START_REF] Ali | Content of the frustum of a simplex[END_REF]), which allows to compute ∆[(• -y) n-1

+ : a 0 , . . . , a n ] in O(n 2 ) time.
Rename as b 1 , . . . , b r the elements a i such that a i < y and as c 1 , . . . , c s the elements a i such that a i y so that r + s = n + 1. Then, the unique solution of the recurrence equation

α k,l = (c l -y)α k-1,l + (y -b k )α k,l-1 c l -b k (k r, l s), (12)
with initial values α 1,1 = (c 1 -b 1 ) -1 and α 0,l = α k,0 = 0 for all l, k 2, is given by

α k,l := ∆[(• -y) k+l-2 + : b 1 , . . . , b k , c 1 , . . . , c l ] (k + l 2).
In order to compute ∆[(• -y) n-1 + : a 0 , . . . , a n ] = α r,s , it suffices therefore to compute the sequence α k,l for k + l 2, k r, l s, by means of two nested loops, one on k, the other on l.

We can compute ∆[(• -y) n -: a 0 , . . . , a n ] similarly. Indeed, the same recurrence equation applied to the initial values α 0,l = 0 for all l 1 and α k,0 = 1 for all k 1, produces the solution

α k,l := ∆[(• -y) k+l-1 - : b 1 , . . . , b k , c 1 , . . . , c l ] (k + l 1).
See for instance [START_REF] Ali | Content of the frustum of a simplex[END_REF] and [START_REF] Varsi | The multidimensional content of the frustum of the simplex[END_REF] for further details.

Application to aggregation theory

As we have already mentioned, the concept of linear combination of lattice polynomials, when it is nondecreasing in each variable, is known in aggregation theory as the discrete Choquet integral, which is extensively used in nonadditive expected utility theory, cooperative game theory, complexity analysis, measure theory, etc. (see [START_REF] Grabisch | Fuzzy measures and integrals[END_REF] for an overview.) For instance, when a discrete Choquet integral is used as an aggregation tool in a given decision making problem, it is then very informative for the decision maker to know its distribution. In that context, the most natural a priori density on [0, 1] n is the uniform one, which makes the results derived here of particular interest.

Example 6 Let h : [0, 1] 3 → R be a discrete Choquet integral defined by v h ({1}) = 0.1, v h ({2}) = 0.6, v h ({3}) = v h ({1, 2}) = v h ({1, 3}) = v h ({2, 3}) = 0.9, and v h ({1, 2, 3}) = 1. According to [START_REF] Algaba | The Lovász extension of market games[END_REF], it can be written as h(x) = 0.1 x 1 + 0.6 x 2 + 0.9 x 3 + 0.2(x 1 ∧ x 2 ) -0.1(x 1 ∧ x 3 ) -0.6(x 2 ∧ x 3 ) -0.1(x 1 ∧ x 2 ∧ x 3 ).

Its density, which can be computed through [START_REF] Devore | Constructive approximation[END_REF] and the recurrence equation [START_REF] Gerber | The volume cut off a simplex by a half-space[END_REF], is represented in Figure 1 by the solid line. The dotted line represents the density estimated by the kernel method from 10 000 randomly generated realizations. The typical value and standard deviation can also be calculated through the raw moments: we have

E[Y h ] ≈ 0.608 and E[Y 2 h ] -E[Y h ] 2 ≈ 0.204.

ACCEPTED MANUSCRIPT

Fig. 1. Density of a discrete Choquet integral (solid line).


Feike C Drost

Ramon Van Den Akker

• Bas

J M Werker

Note on integer-valued bilinear time series models

Keywords: count data, integer-valued time series, bilinear model

.

Introduction

In many sciences one encounters nonnegative discrete valued time series, often as counts of events or objects at consecutive points in time. Especially in economics and medicine many interesting variables are (nonnegative) integer-valued. For example: the number of transactions in IBM during each minute, the number of patients in a hospital at the end of the day, the number of claims an insurance company receives during each day, the number of epileptic seizures a patient suffers each day, etcetera. Until the mid seventies modeling discrete valued time series did not attract much attention since most traditional representations of dependence become either impossible or impractical. The last two decades there were many developments in the literature on integer-valued time series; see [START_REF] Mckenzie | Discrete variate time series[END_REF] for a detailed review; for example, INteger-valued Moving Average processes (Al-Osh and [START_REF] Al-Osh | Binomial autoregressive moving average models[END_REF]), INteger-valued AutoRegressive processes (Al-Osh and [START_REF] Al-Osh | First-order integer-valued autoregressive (INAR(1)) processes[END_REF], Al-Osh and [START_REF] Alzaid | An integer-valued pth-order autoregressive structure (INAR(p)) process[END_REF], and [START_REF] Du | The integer-valued autoregressive (INAR(p)) model[END_REF]), and Generalized INteger-valued AutoRegressive processes [START_REF] Latour | Existence and stochastic structure of a non-negative integer-valued autoregressive process[END_REF]) are common choices. [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF] introduced the class of nonnegative INteger-valued BiLinear time series, which contains the three aforementioned classes. Higher order processes turn out to be highly releveant see [START_REF] Latour | Existence and stochastic structure of a non-negative integer-valued autoregressive process[END_REF]) for the class of GINAR models. For extensions to regression based count data time series models see, for example, [START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF] and [START_REF] Kedem | Regression Models for Time Series Analysis[END_REF].

Recall the definition of an INBL(p, q, m, n) process: let (X -1 , . . . , X -p∨m , ε -1 , . . . , ε -q∨n ) be the (nonnegative) integer-valued initial points generated by some probability distribution ν, and define X t , t = 0, 1, 2, . . ., recursively by,

X t = p i=1 α i • X t-i + q j=1 β j • ε t-j + m k=1 n =1 γ k, • (X t-k ε t-) + ε t , (1)
where (ε t ) t≥0 is a collection of i.i.d. nonnegative integer-valued variables with distribution G with mean µ G ∈ [0, ∞] and variance σ 2 G ∈ [0, ∞], and where the Steutel-van Harn operators η• applied to some nonnegative integer-valued random variable Z are defined as

η • Z = Z s=1 U (η) s ,
where (U

(η)
s ) s≥1 is a collection of i.i.d. nonnegative integer-valued variables with distribution function F η with mean η. * In proofs we suppose, without loss of generality, p = q = m = n since we may introduce additional lags in (1) with η = 0 or equivalently F η {0} = 1. The class of GINAR(p) processes arises by taking q = m = n = 0. Specializing further to the situation where the F η -distributions are Bernoulli, one gets the class of INAR(p) models. Similarly INMA(q) processes are obtained when p = m = n = 0.

In this note we consider stationarity and existence of moments of superdiagonal INBL(p, q, m, n) processes with γ k, = 0 if k < . Using Markov chain techniques we obtain the existence and uniqueness of a strictly stationary solution of (1) when

i α i + µ G k, γ k, < 1 (see Theorem 2.1
). Theorem 2.2 provides sufficient parameter restrictions ensuring the existence of (higher-order) moments under the stationary distribution. Along completely different lines our results generalize Sections 2 and 3 in [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF] who restricted attention to the INBL(1, 0, 1, 1) process. Moreover our parameter restrictions in Theorem 2.2 are less severe. Estimation of superdiagonal INBL processes can be performed along the same lines as in [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF]; The details are beyond the scope of this note.

Results

First we discuss the existence of a stationary solution. The insightful proof of Theorem 2.1 in [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF] gives an explicit construction of approximating processes for the INBL(1, 0, 1, 1) process. For the general superdiagonal INBL(p, q, m, n) case we provide a compact proof using Markov chain techniques, yielding also an alternative proof for the INBL(1, 0, 1, 1) process.

Let F * G denote the convolution of the probability measures F and G and introduce the process Z = (Z t ) t≥0 by

Z t = (X t-1 , . . . , X t-p∨m , ε t-1 , . . . , ε t-q∨n ).
(2)

Notice that Z 0 ∼ ν. It is easy to see that Z is a Markov chain. For t ≥ 0 and z t = (x t-1 , . . . , x t-p∨m , e t-1 , . . . , e t-q∨n ) we have

P(Z t+1 = z t+1 | Z t = z t ) = G{e} * p i=1 ( * x t-i s=1 F α i ) * q j=1 ( * e t-j s=1 F β j ) * m k=1 * n =1 ( * x t-k e t- s=1
F γ k, ) {z -e}, * Implicitly it is understood that distinct Steutel-van Harn operators are independent and also that Steutel-van Harn operators applied to different random variables are independent; see also [START_REF] Du | The integer-valued autoregressive (INAR(p)) model[END_REF].

if z t+1 = (z, x t-1 , . . . , x t+1-p∨m , e, e t-1 , . . . , e t+1-q∨n ) for some e, z, and with

P(Z t+1 = z t+1 | Z t = z t ) = 0 otherwise. From (1) it is immediate that, for t ≥ 0, E[X t | Z t ] = µ G + p i=1 α i X t-i + q j=1 β j ε t-j + m k=1 n =1 γ k, X t-k ε t-∈ [0, ∞].
(3)

The following proposition provides sufficient conditions for Z to be an irreducible aperiodic Markov chain. The proof is immediate from the definition of the state space.

Proposition 2.1 Let Z be the time series defined in (2) linked to the INBL X defined in (1) and assume, for all relevant i, j, k, ,

0 < G{0}, 0 < F αi {0}, 0 < F βj {0}, 0 < F γ k, {0}. ( 4
)
Then Z is an irreducible aperiodic Markov chain on the state space

S = {z | ∃t ≥ 1 : P(Z t = z | Z 0 = 0) > 0}.
Remark Condition (4) allows the process X to arrive at 0 and therefore (4) seems to be harmless in real-life applications. Condition (4) can be easily verified from a sufficiently long data series just by checking whether the process has hit zero. At the cost of lengthy derivations, condition (4) might be weakened but this is outside the scope of this short note.

Theorem 2.1 Let X be the superdiagonal INBL process as defined in (1) with γ k, = 0 if k < and assume (4), σ 2 G < ∞, q j=1 β j < ∞, and

p i=1 α i + µ G m k=1 n =1 γ k, < 1. ( 5
)
Then there exists a unique initial distribution ν * such that X is strictly stationary under P ν * . Furthermore the first moment of X is finite under P ν * .

Proof. We prove that there exists a unique initial distribution ν * such that Z is stationary under P ν * , which yields the result. Since Z is irreducible and aperiodic on S (and since, for countable Markov chains, finite sets are petite) it suffices to prove, see Theorem 15.0.1 in [START_REF] Meyn | Markov chains and stochastic stability[END_REF], that there exists a mapping V : z → [1, ∞) such that the following Foster-Lyapunov drift criterium holds: there exists δ > 0 and a finite set such that for all Z t outside this finite set we have

W (Z t ) = (1 -δ)V (Z t ) -E[V (Z t+1 ) | Z t ] ≥ 0. ( 6
)
Take, for notational simplicity, p = q = m = n. To verify (6) choose δ > 0 sufficiently small such that

i α i + µ G k≥ γ k, < (1 -δ) n -nδ - 1 2 n(n + 1)µ G δ,
and define V : z → [1, ∞) by V (Z t ) = 1 + i αi X t-i + j βj ε t-j + k≥ γk, X t-k ε t-,
where the nonnegative tilde parameters are recursively defined from αn+1 = βn+1 = γn+1, = 0 and

(1 -δ)α i = αi+1 + µ G γi+1,1 + α i + δ, (1 -δ) βj = βj+1 + β j + δ, (1 -δ)γ k, = γk+1, +1 + γ k, + δ.
Note, by the choice of δ,

θ ≡ α1 + µ G γ1,1 = i (α i + δ)/(1 -δ) i + µ G k≥ (γ k, + δ)/(1 -δ) k < 1.
Hence we obtain, using ( 1), ( 3), and the independence of ε t and X t -ε t ,

W (Z t ) = -δ -µ G (α 1 + β1 ) -(µ 2 G + σ 2 G )γ 1,1 + i [(1 -δ)α i -(α 1 + µ G γ1,1 )α i -αi+1 -µ G γi+1,1 ]X t-i + j [(1 -δ) βj -(α 1 + µ G γ1,1 )β j -βj+1 ]ε t-j + k≥ [(1 -δ)γ k, -(α 1 + µ G γ1,1 )γ k, -γk+1, +1 ]X t-k ε t- = -δ -µ G (α 1 + β1 ) -(µ 2 G + σ 2 G )γ 1,1 + (1 -θ) ×   i (α i + δ 1 -θ )X t-i + j (β j + δ 1 -θ )ε t-j + k≥ (γ k, + δ 1 -θ )X t-k ε t-   .
Conclude W (Z t ) ≥ 0 outside a finite set. Hence the drift condition (6) holds, which concludes the proof of the stationarity. The existence of the first moment under ν * is immediate from Theorem 2.2 below. 2

Remark For p = m = n = 1 and q = 0, (5) reduces to the condition of Theorem 2.1 in [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF].

The next theorem gives a sufficient condition for the existence of higher order moments of X t under P ν * .

Theorem 2.2 Let X be the superdiagonal INBL process as defined in (1) with γ k, = 0 if k < . Suppose, for some K ∈ N, the existence of the K-th order moments of F αi , F βj , and F γ k, , and assume ( 4),

E G ε 2K 0 < ∞, q j=1 β j < ∞, and p i=1 α i + E G ε K 0 1/K m k=1 n =1 γ k, < 1. (7)
Then, with ν * the stationary distribution of Theorem 2.2,

E ν * X K 0 < ∞.
Proof. Recall from Proposition 2.1 that Z is an irreducible aperiodic Markov chain on S and 0 ∈ S. Hence, since Z has stationary distribution ν * , we have Z t d -→ ν * , under P δ 0 as t → ∞. Thus the convergence also holds for marginals. By the Portmanteau theorem, and nonnegativity of X, this implies for L ∈ N,

E ν * X L-1 0 ≤ lim inf t→∞ E δ0 X L-1 t ≤ sup 0≤t<∞ E δ0 X L-1 t . ( 8
)
To prove that the right hand side of ( 8) is bounded for L ≤ K + 1, we use induction. Clearly the statement holds true for L = 1. To show the boundedness for L = 2 we make some preliminary remarks. Denote, for notational convenience, E δ 0 by E 0 and let, for P ≥ 1,

Z P = E 0 |Z| P 1/P denote the L P (P δ 0 ) norm. Recall • P ≤ • Q for 1 ≤ P ≤ Q. Hence (7) implies that i α i + ε 0 v k, γ k, < 1 for all 1 ≤ v ≤ K.
Since, for t ≤ s, X t -ε t and ε s are independent, we have, for P ≥ 1,

X t ε s P ≤ (X t -ε t )ε s P + ε t ε s P ≤ ε 0 P X t P + ε 2 0 P . ( 9
)
Using ( 3), ( 9) with P = 1, the zero starting values, and ( 7), the boundedness for L = 2 is obtained from,

E 0 X t = µ G + i α i E 0 X t-i + j β j µ G + k≥ γ k, E 0 X t-k ε t- ≤ µ G + µ G j β j + (µ 2 G + σ 2 G ) k≥ γ k, + i α i E 0 X t-i + µ G k≥ γ k, E 0 X t-k ≤   µ G + µ G j β j + (µ 2 G + σ 2 G ) k≥ γ k,   ∞ s=0   i α i + µ G k≥ γ k,   s < ∞.
To complete the induction, we assume that the right hand side of ( 8) is bounded for some 2 ≤ L ≤ K and we prove that this implies sup 0≤t<∞ E δ0 X L t < ∞. We use the following result (see, for example, [START_REF] Dharmadhikari | Bounds on the moments of martingales[END_REF]).

Lemma 2.1 If Z 1 , . . . , Z n are i.i.d. variables with mean zero and a finite k-th moment, k ≥ 2, then we have the bound

E n s=1 Z s k ≤ C k n k/2 E|Z 1 | k ,
where the constant C k > 0 only depends on k (and not on the distribution of Z 1 ).

Using that the Steutel-van Harn operator

η • Z, conditional on Z, follows a * Z s=1 F η distri- bution, Lemma 2.1 implies the following inequality for L ≥ 2, E|η • Z -ηZ| L = EE |η • Z -ηZ| L | Z ≤ C L EZ L/2 E|U (η) -η| L ≤ C L,η Z L/2 L-1 . (10)
We have for t ≥ 0, using (9), van den Akker, and Werker

X t L ≤ X t - i α i X t-i - j β j ε t-j - k≥ γ k, X t-k ε t-L Drost,
+ i α i X t-i + j β j ε t-j + k≥ γ k, X t-k ε t-L ≤ i α i • X t-i -α i X t-i L + j β j • ε t-j -β j ε t-j L + k≥ γ k, • (X t-k ε t-) -γ k, X t-k ε t-L + ε t L + i α i X t-i L + j β j ε t-j L + k≥ γ k, ( ε 0 L X t-k L + ε 2 0 L ) ≤ M + i α i X t-i L + ε 0 L k≥ γ k, X t-k L ,
where the appropriately chosen constant M < ∞, not depending on t, can be obtained from the assumed moment conditions and the induction hypothesis applied to (10) and ( 9). Condition (7) completes the induction argument just as for the case L = 2,

X t L ≤ M ∞ s=0   i α i + ε 0 L k≥ γ k,   s .
This completes the proof. 2

Remark For the INBL(1,0,1,1) process [START_REF] Doukhan | A simple integer-valued bilinear time series model[END_REF] showed that the condition U (α) K + ε 0 K U (γ) K < 1 suffices to obtain the existence of the K-th moment of X t under P ν * . Condition ( 7) is weaker since, for K ≥ 1, η = U (η) 1 ≤ U (η) K .

Drost, van den Akker, and Werker

Note on integer-valued bilinear time series models


L Baringhaus

N Henze

A new weighted integral goodness of fit statistic for exponentiality

Keywords: Exponential distribution, mean residual life function, Brownian bridge, Bessel functions AMS 1991 subject classifications: Primary 62G10, Secondary 62E20

We propose a new weighted integral goodness of fit statistic for exponentiality. The statistic is motivated by a characterization of the exponential distribution via the mean residual life function. Its limit null distribution is seen to be the same as that of a certain weighted integral of the squared Brownian bridge. The Laplace transform and cumulants of the latter are expressible in terms of Bessel functions.

Introduction

Let X 1 , . . . , X n be independent and identically distributed non negative random variables with some unknown distribution function F (x) = P (X 1 ≤ x), x ≥ 0. It is well-known that, under the assumption 0 < E(X 1 ) < ∞, the distribution of X 1 is exponential, i.e. F (x) = 1 -exp(-λx), x ≥ 0, for some λ > 0, if and only if the mean residual life function is constant, i.e. we have E(X 1 -z|X 1 > z) = E(X 1 ) for each z > 0.

(1) [START_REF] Baringhaus | Tests of fit for exponentiality based on a characterization via the mean residual life function[END_REF] noted that (1) is equivalent to E(min(X 1 , z)) = E(X 1 )F (z) for each z > 0.

(2)

Arguing that, under the assumptions X 1 ≥ 0 and 0 < E(X 1 ) < ∞, (2) is a characteristic property of the class {Exp(λ) : λ > 0} of exponential distributions, they suggested a new approach to assess exponentiality. In particular, they proposed the Cramér-von Mises type statistic

G n = n ∞ 0 1 n n k=1 min(U k , z) - 1 n n k=1 I(U k ≤ z) 2 e -z dz, (3)
where X = 1 n n k=1 X k and U k = X k /X, k = 1, . . . , n. Interestingly, the limit null distribution of G n is the same as the limit null distribution of the classical Cramér-von Mises statistic when testing the simple hypothesis of uniformity on the interval [0, 1]. The goodness of fit test based on G n rejecting the hypothesis of exponentiality if G n > g n where g n is the (1 -α)-quantile of the null distribution of G n for some given α ∈ (0, 1), is consistent against any fixed alternative distribution. The test is discussed in the recent review paper on goodness of fit tests for exponentiality by [START_REF] Henze | Recent and classical tests for exponentiality: a partial review with comparisons[END_REF]. From the comparative simulation study given there one may conclude that it is a a serious competitor, although there are other procedures showing a better power performance for special alternatives. However, the favorable behavior of these procedures is recognized by varying some weight parameter built in the corresponding test statistic. So, to obtain a possible gain in power performance we suggest to generalize (3) in a natural way by using the more general weight function e -az , z ≥ 0, where a > -1 is some real parameter. Then we have the test statistic

G n,a = n ∞ 0 1 n n k=1 min(U k , z) - 1 n n k=1 I(U k ≤ z) 2 e -az dz. (4)
By some elementary calculations we see that in the case a = 0 the statistic can be written in the form

G n,a = 1 n n k, =1 2 a 3 - 1 a 2 + 1 a min(U k , U ) exp(-aU k ) + exp(-aU ) - 2 a 3 + 1 a 2 exp(-a min(U k , U )) + 1 a 2 + 1 a exp(-a max(U k , U )) .

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

In the special case a = 0 the statistic has the alternative representation

G n,0 = 1 n n k, =1 1 3 min(U k , U ) 3 + min(U k , U ) -1 max(U k , U ) - 1 2 U 2 k - 1 2 U 2 .
Using the central limit theorem in Hilbert spaces as a method of proof different to that given in [START_REF] Baringhaus | Tests of fit for exponentiality based on a characterization via the mean residual life function[END_REF] we show in the next section, that the limit null distribution of G n,a is the same as that of

1 0 B(t) 2 (1 -t) a-1 dt, where (B(t), 0 ≤ t ≤ 1
) is a Brownian bridge. Using this representation we shall be able to express the Laplace transform of the limit null distribution in terms of the Bessel function I 1 a+1 , and, additionally, its cumulants via sums of the even powers of the reciprocal zeros of the Bessel function J 1

a+1

. Critical values and empirical power values of the test obtained by simulation are shown in the last section.

The limit null distribution of G n,a

The null distribution of G n,a does not depend on the parameter λ of the underlying exponential distribution. We therefore assume that λ = 1, i.e. that X k has the distribution function

1 -exp(-x), x ≥ 0. Putting R n (z) = √ n 1 n n k=1 min(X k , z) -XI(X k ≤ z) , z ≥ 0, it follows that G n,a = 1 X 3 ∞ 0 R 2 n (z) exp - 1 X -1 az -1 e -az dz + 1 X 3 ∞ 0 R 2 n (z)e -az dz. (5)
We shall show below that for each a > -1 the statistic ∞ 0 R 2 n (z) exp(-az) dz has a limit distribution. Using this in advance we first prove that

∞ 0 R 2 n (z) exp - 1 X -1 az -1 e -az dz = o P (1). (6)
For, let us treat the cases a ≥ 0 and -1 < a < 0 separately. If a ≥ 0 we put η = 1 2 , a * = -η and 0 = 1. For given 0 < ≤ 0 we consider the case where

| 1 X -1| < . From exp - 1 X -1 az -1 ≤ za exp( az), z ≥ 0, we get using z ≤ e ηz /η for z ≥ 0 that ∞ 0 R 2 n (z) exp - 1 X -1 az -1 e -az dz ≤ a ∞ 0 R 2 n (z)z exp(-a(1 -)z) dz ≤ a η ∞ 0 R 2 n (z) exp(-[a(1 -) -η]z) dz. ≤ a η ∞ 0 R 2 n (z) exp(-a * z) dz.
If -1 < a < 0 we put η = a+1 4 , a * = a-3 4 and 0 = 1 2 -1 a -1 . For given 0 < ≤ 0 we consider the case where

| 1 X -1| < . From exp - 1 X -1 az -1 ≤ z|a| exp(-az), z ≥ 0, we obtain ∞ 0 R 2 n (z) exp -( 1 X -1)az -1 e -az dz ≤ |a| ∞ 0 R 2 n (z)z exp(-a(1 + )z) dz ≤ |a| η ∞ 0 R 2 n (z) exp(-[a(1 + ) -η]z) dz. ≤ |a| η ∞ 0 R 2 n (z) exp(-a * z) dz.
Thus in any case there is some a * > -1 such for each 0 and1 X 3 = 1 + o P (1) we obtain (6) and, additionally,

< ≤ 0 ∞ 0 R 2 n (z) exp -( 1 X -1)az -1 exp(-az) dz ≤ |a| η ∞ 0 R 2 n (z) exp(-a * z) dz if | 1 X -1| < . Using that ∞ 0 R 2 n (z) exp(-a * z) dz has a limit distribution, 1 X = 1 + o P (1)
G n,a = ∞ 0 R 2 n (z) exp(-az) dz + o P (1).
It remains to derive the limit distribution of ∞ 0 R 2 n (z)e -az dz. To this end note that

R n (z) = H n (z) + (X -1)L n (z), z ≥ 0,
where

H n (z) = √ n 1 n n k=1 [min(X k , z) -I(X k ≤ z) -(X k -1)(1 -e -z )] , z ≥ 0,
and

L n (z) = √ n 1 n n k=1 1 -e -z -I(X k ≤ z) , z ≥ 0.
H n and L n are random elements in the Hilbert space L 2 (R + , B + , µ a ), where R + = [0, ∞), B + is the Borel σ-field on R + , and µ a is the σ-finite measure having density exp(-ax), x ≥ 0, with respect to Lebesgue measure on (R

+ , B + ). Since ∞ 0 Var min(X 1 , z)-I(X 1 ≤ z)-(X 1 -1)(1-e -z ) dµ a (z) = ∞ 0 [1-e -z ]e -(a+1)z dz < ∞ A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT and ∞ 0 Var 1 -e -z -I(X 1 ≤ z) dµ a (z) = ∞ 0 [1 -e -z ]e -(a+1)z dz < ∞
the central limit theorem for Hilbert space-valued random variables applies (see [START_REF] Ledoux | Probability in Banach Spaces[END_REF], Corollary 10.9). Recognizing that

ρ(w, z) = min(1 -e -w , 1 -e -z ) -(1 -e -w )(1 -e -z ), w, z ≥ 0, (7)
is the covariance function of the process (H n (z), z ≥ 0) and also that of the process (L n (z), z ≥ 0), there is a zero mean Gaussian process (H(z), z ≥ 0) with sample paths in L 2 (R + , B + , µ a ) and covariance kernel ( 7) such that

H 2 n (z) dµ a (z) D -→ H 2 (z) dµ a (z) and L 2 n (z) dµ a (z) D -→ H 2 (z) dµ a (z),
where "

D -→" means convergence in distribution. Due to X -1 = o P (1) this implies that ∞ 0 R 2 n (z) exp(-az) dz D -→ ∞ 0 H 2 (z) exp(-az) dz.
The result just proved is summarized as follows.

Theorem 1. The limit null distribution of the test statistic G n,a is that of H 2 (z) dµ a (z) where (H(z), z ≥ 0) is a zero mean Gaussian process with sample paths in L 2 (R + , B + , µ a ) and covariance function (7).

The Gaussian process (H(z), z ≥ 0) has the same covariance function as the process B(1 -e -z ), z ≥ 0 where (B(t), 0 ≤ t ≤ 1) is the Brownian bridge. Consequently,

H 2 (z) dµ a (z) D = 1 0 B 2 (1 -t)t a-1 dt D = 1 0 B 2 (t)(1 -t) a-1 dt
where " D =" denotes equality in distribution. Since the limit distribution may be of independent interest, it is studied further. Let

G a = 1 0 B(t) 2 (1 -t) a-1 dt. Theorem 2. a) The Laplace transform of G a is ψ a (z) = ∞ j=1 1 + 2 a + 1 2 2z γ 2 j -1/2 = Γ 1 + 1 a + 1 1 1 a+1 √ 2z 1 a+1 I 1 a+1 2 a + 1 √ 2z -1/2 , z ≥ 0, A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

where I ν is the modified Bessel function of the first kind of order ν.

b) The first four cumulants of G a are

κ 1 = E(G a ) = 1 (a + 1)(a + 2) , κ 2 = Var(G a ) = 2 (a + 2) 2 (a + 1)(2a + 3) , κ 3 = 16 (a + 2) 3 (a + 1)(2a + 3)(3a + 4) , κ 4 = 48(11a + 16) (a + 2) 4 (2a + 3) 2 (a + 1)(3a + 4)(4a + 5) .
Proof. Although part a) of the theorem can be obtained from the material presented in [START_REF] Deheuvels | Karhunen-Loéve expansions for weighted Wiener processes and Brownian bridges via Bessel functions[END_REF], we give a derivation for the special case considered here. Note that

G a = B 2 (t) dν a (t),
where ν a is the σ-finite measure with density (1 -t) a-1 with respect to Lebesgue measure on the Borel sets

B [0,1] of the interval [0, 1]. Let A be the integral operator on L 2 ([0, 1], B [0,1] , ν a ) associated with the covariance kernel k(s, t) = min(s, t) -st, 0 ≤ s, t ≤ 1, of the Brownian bridge. Thus Af (t) = 1 0 k(t, s)f (s) dν a (s), 0 ≤ t ≤ 1, for f ∈ L 2 ([0, 1], B [0,1] , ν a ).
The integral operator A is positive definite. It is well known (see, e.g. [START_REF] Vakhania | Probability Distributions on Linear Spaces[END_REF], p. 58) that

G a D = k≥1 λ k Z 2 k , (8)
where the Z k are independent unit normal random variables, and the λ k are the eigenvalues of A. To obtain these eigenvalues, assume that an eigenfunction f of A with associated positive eigenvalue λ is smooth enough so that, starting with the equation

k(t, s)f (s) dν a (s) = λf (t), 0 ≤ t ≤ 1, (9)
we may differentiate twice on both sides of (9). This leads to the differential equation

λf (t) = -(1 -t) a-1 f (t), 0 < t < 1.
From (9) we also infer the boundary conditions f (0) = 0 and f (1) = 0. Putting ϕ(t) = f (1 -t), ϕ satisfies the differential equation

λϕ (t) = -t a-1 ϕ(t), 0 < t < 1, (10)
A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

subject to the boundary conditions

ϕ(0) = ϕ(1) = 0. ( 11
)
The general solution of (10) was found by Lommel (1868); for another reference see [START_REF] Nielsen | Handbuch der Theorie der Zylinderfunktionen[END_REF], p. 130. In view of (11), the general solution is

ct 1/2 J 1 a+1 2 a + 1 λ -1/2 t (a+1)/2 , 0 ≤ t ≤ 1,
where J ν is the Bessel function of the first kind of order ν, and c is some constant. In what follows, 0 < γ 1 < γ 2 < . . . are the positive zeros of J 1 a+1

. Putting

λ j = 2 a + 1 2 1 γ 2 j , j = 1, 2, . . . , and
f j (t) = (a + 1) 1/2 J a+2 a+1 γ j ) J 1 a+1 (γ j (1 -t) (a+1)/2 (1 -t) 1/2 , (0 ≤ t ≤ 1, j = 1, 2, . . . ), it follows from f j (t)f k (t) dν a (t) = a + 1 J a+2 a+1 (γ j )J a+2 a+1 (γ k ) 1 0 (1 -t) a J 1 a+1 γ j (1 -t) (a+1)/2 J 1 a+1 γ k (1 -t) (a+1)/2 dt = 0, j = k, 1, j = k,
(see [START_REF] Erdélyi | Higher transcendental functions[END_REF], p. 70) that {f j , j ≥ 1} is an orthonormal set in

L 2 ([0, 1], B [0,1] , ν a ). For f ∈ L 2 ([0, 1], B [0,1] , ν a ) the Fourier-Bessel series j≥1 α j f j with α j = f (t)f j (t) dν a (t), j ≥ 1, converges in L 2 ([0, 1], B [0,1] , ν a ), i.e. lim n→∞ n j=1 α j f j (t) -f (t) 2 dν a (t) = 0.
This result can be proved in a rather elementary way by using the work of [START_REF] Hochstadt | The mean convergence of Fourier-Bessel series[END_REF], for example. Thus

{f j , j ≥ 1} is a complete orthonormal set in L 2 ([0, 1], B [0,1] , ν a ). Using -f j (t)(1 -t) a-1 = λ j f j (t), 0 < t < 1, A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

integration by parts gives

(1 -t) t 0 sf j (s)(1 -s) a-1 ds + t 1 t f (s)(1 -s) a ds = λ j f j (t), 0 ≤ t ≤ 1, or, equivalently, k(t, s)f j (s) dν a (s) = λ j f j (t), 0 ≤ t ≤ 1.
Thus the f j , j ≥ 1, form a complete orthonormal system of eigenfunctions of A with associated eigenvalues λ j , j ≥ 1. In view of (8) the Laplace transform of G a is

ψ a (z) = ∞ j=1 1 + 2 a + 1 2 2z γ 2 j -1/2 , z ≥ 0, = Γ 1 + 1 a + 1 1 1 a+1 √ 2z 1 a+1 I 1 a+1 2 a + 1 √ 2z -1/2 , z ≥ 0,
where I ν is the modified Bessel function of the first kind of order ν. For the second equality see [START_REF] Nielsen | Handbuch der Theorie der Zylinderfunktionen[END_REF], p.358. This proves part a) of the theorem. To prove part b) which asserts an interesting connection between the cumulants cumulants κ ν of G a , the latter being an integral involving the Brownian bridge, and the sums of the even powers of the reciprocal zeros γ j of the Bessel functions J 1 a+1 we note that

κ ν = 2 ν-1 (ν -1)! ∞ j=1 λ ν j = 2 ν-1 (ν -1)! 2 a + 1 2ν ∞ j=1 1 γ j 2ν , ν ≥ 1.
We can derive the first four cumulants of G a directly by using Fubini's theorem and calculating mixed moments of a zero mean multivariate normal vector with covariance structure specified by the Brownian bridge. Alternatively, we can use that [START_REF] Nielsen | Handbuch der Theorie der Zylinderfunktionen[END_REF], p. 360) to obtain the first four cumulants stated in part b) of the theorem.

∞ j=1 1 γ 2 j = 2 -2 1 + 1 a+1 ∞ j=1 1 γ 4 j = 2 -4 (1 + 1 a+1 ) 2 ( 1 a+1 + 2) ∞ j=1 1 γ 6 j = 2 -6 • 2 (1 + 1 a+1 ) 3 ( 1 a+1 + 2)( 1 a+1 + 3) ∞ j=1 1 γ 8 j = 2 -8 (5 1 a+1 + 11) (1 + 1 a+1 ) 4 ( 1 a+1 + 2) 2 ( 1 a+1 + 3)( 1 a+1 + 4) (see

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

It is also possible to derive 'limit test statistics' in the cases a → -1 and a → ∞. For, it is immediately seen that

G n,-1 = lim a↓-1 G n,a = -2n + 1 n n k, exp min(U k , U ) .
Of course, one can use G n,-1 as test statistic for testing the hypothesis of exponentiality.

However, presently we are unable to present any asymptotic theory for G n,-1 as n → ∞.

In the second case, a → ∞, putting

U 1:n = min(U 1 , . . . , U n ) we have lim a→∞ a -2 exp(aU 1:n ) a 3 G n,a n -2 = 1 n -2U 1:n .
Thus one may suggest U 1:n as test statistic rejecting hypothesis of exponentiality for small values of U 1:n . We remark that a test based on this statistic is easily done because U 1:n has the beta distribution B(1, n -1) with the density (n -1)(1 -t) n-2 , 0 ≤ t ≤ 1. Finally, we remark that for each -1 < a < ∞ and given level α the test obtained by rejecting the hypothesis of exponentiality if G n,a > g n,a (α), where g n,a (α) is the (1 -α)quantile of G n,a in the case where the hypothesis of exponentiality is true, is consistent against any fixed alternative distribution. A proof of this assertion is easily done by adapting the arguments used by [START_REF] Baringhaus | Tests of fit for exponentiality based on a characterization via the mean residual life function[END_REF].

Empirical Results

For sample sizes n = 10, 20,30,40,50,100,200, parameter values a = -0.99, -0.9, -0.5, 0, 0.5, 1, 1.5, 2, 5, 10 and levels α = 0.05, 0.1 we got approximations of the critical values g n,a (α) by simulation with 100000 replications. The results are shown in TABLE 1 and TABLE 2 of the appendix. An empirical power study based on simulations with 10000 replications was done for the sample sizes n = 20, 50, 100 and level α = 0.05. The alternative distributions were chosen from distribution families considered also by Baringhaus andHenze (1991), (2000). The distributions included are the Gamma distributions (G), Weibull distributions (W), Lognormal distributions (LN ) with scale parameter 1 and shape parameter θ, the uniform distribution U[0, 1], the Half-Normal distribution (HN ), the Half-Cauchy distribution (HC), the χ 2 1 -distribution, the Power distributions (PW) with density θ -1 x 1/θ-1 , 0 < x < 1, the linear increasing failure rate distributions (LIFR) with density (1 + θx) exp(-(x + θx 2 /2)), x > 0, and the JSHAPE distributions (JS) with density (1 + θx) -1/θ-1 , x > 0. The empirical power values (rounded to the nearest integer) are shown in TABLE 3, TABLE 4 and TABLE 5 of the appendix. For estimated power values of various other competitive procedures we refer to [START_REF] Baringhaus | A class of consistent tests for exponentiality based on the empirical Laplace transform[END_REF] and also to [START_REF] Henze | Recent and classical tests for exponentiality: a partial review with comparisons[END_REF], although there is merely a partial overlap with the alternative distributions considered in the latter paper. We conclude, that taking a weight parameter 1 ≤ a ≤ 2 is a rather good choice. For some alternative distributions, the Gamma distributions or lognormal distributions with shape parameter θ < 1 for example, a choice of a > 2 can be recommended. In these cases the power performance is also better than that of the top ranked test in the comparative study of [START_REF] Henze | Recent and classical tests for exponentiality: a partial review with comparisons[END_REF].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

A Appendix: Tables of critical values and power values

TABLE 1

 1 Critical values of G n,a for α = 0.05

	Sample					Parameter a
	size n	-0.99	-0.9	-0.5	0	0.5	1	1.5	2	5	10
	10	5.056	3.944 1.846 1.017 0.634 0.427 0.305 0.226 0.065 0.018
	20	8.069	5.983 2.336 1.147 0.678 0.449 0.313 0.230 0.067 0.021
	30	9.886	7.282 2.597 1.193 0.691 0.450 0.316 0.231 0.068 0.021
	40	11.299 8.107 2.723 1.224 0.702 0.446 0.316 0.234 0.069 0.022
	50	12.433 8.793 2.846 1.249 0.699 0.453 0.320 0.232 0.069 0.022
	100	15.848 10.654 3.072 1.279 0.713 0.457 0.316 0.233 0.069 0.022
	200	18.338 11.900 3.233 1.296 0.711 0.464 0.320 0.235 0.069 0.022

TABLE 2

 2 Critical values of G n,a for α = 0.1

	Sample					Parameter a
	size n	-0.99	-0.9	-0.5	0	0.5	1	1.5	2	5	10
	10	3.063 2.580 1.406 0.793 0.490 0.332 0.237 0.174 0.049 0.014
	20	4.395 3.526 1.705 0.875 0.515 0.340 0.238 0.176 0.050 0.016
	30	5.341 4.180 1.847 0.908 0.528 0.341 0.240 0.176 0.050 0.016
	40	6.033 4.581 1.945 0.927 0.531 0.343 0.240 0.176 0.051 0.016
	50	6.569 4.946 2.042 0.948 0.536 0.346 0.238 0.175 0.051 0.016
	100	8.200 5.929 2.226 0.977 0.539 0.346 0.240 0.177 0.051 0.016
	200	9.917 6.961 2.359 0.991 0.544 0.345 0.239 0.176 0.051 0.017

TABLE 3

 3

	Percentage of 10000 Monte Carlo samples declared significant;
		level α = 0.05; sample size n = 20
	Distribution				Parameters a
		-0.99 -0.9 -0.5 0 0.5 1 1.5 2	5 10
	W(0.6)	59	61	65	69 69 70 72 73 78 80
	W(1.2)	1	1	4	9 12 13 14 14 12 8
	W(1.4)	0	0	11	26 33 35 37 37 32 20
	W(1.6)	0	0	26	50 60 63 64 65 57 39
	χ 2 1	38	39	44	47 49 52 55 57 66 72
	PW(0.8)	1	13	87	94 94 93 92 89 69 40
	PW(1.2)	0	1	33	48 48 44 40 35 17 7
	PW(1.4)	0	0	16	28 27 24 21 18 9	7
	PW(2.0)	1	1	4	8 10 11 12 14 25 37
	PW(3.0)	13	17	28	38 46 53 59 64 79 87
	LIFR(1)	0	0	6	14 18 19 19 19 13 7
	LIFR(2)	0	0	10	24 29 30 30 30 20 11
	LIFR(4)	0	0	17	36 42 43 43 42 30 16
	LIFR(6)	0	0	23	43 50 51 52 50 37 20
	LIFR(10)	0	0	30	53 60 61 61 60 47 26
	HN (0, 1)	0	0	7	17 21 21 22 21 14 8
	LN (0.7)	2	3	16	34 46 53 59 64 73 70
	LN (0.8)	6	6	10	19 23 27 31 34 43 39
	LN (1.0)	21	21	20	18 16 15 14 13 11 8
	LN (1.5)	66	67	68	67 67 66 65 64 59 53
	HC	72	72	73	71 69 67 66 65 59 52
	JS(0.5)	50	50	51	49 47 45 44 43 39 35
	JS(1.0)	81	82	84	84 83 83 83 82 80 77
	U[0, 1]	0	3	58	74 74 71 67 62 38 18
	G(0.4)	56	58	65	68 72 75 77 80 86 90
	G(0.6)	26	28	30	30 32 33 35 37 44 50
	G(0.8)	11	11	12	11 10 10 11 11 13 16
	G(1.4)	1	1	4	9 13 14 15 16 14 9
	G(1.6)	0	0	6	16 21 24 26 27 25 17
	G(1.8)	0	0	10	24 31 34 38 39 38 26
	G(2.0)	0	0	15	33 42 46 50 52 49 36
	G(2.4)	0	0	28	53 64 69 72 74 72 57
	G(3.0)	0	1	49	77 85 89 91 92 91 80

TABLE 4

 4

	Percentage of 10000 Monte Carlo samples declared significant;
			level α = 0.05; sample size n = 50
	Distribution				Parameters a
		-0.99 -0.9 -0.5	0	0.5	1	1.5	2	5	10
	W(0.6)	86	88	94	96	97	98	98	98	99	99
	W(1.2)	0	0	8	21	28	30	31	32	29	22
	W(1.4)	0	1	41	66	75	77	78	79	73	61
	W(1.6)	0	5	80	95	97	98	98	98	96	90
	χ 2 1	58	62	75	82	87	89	90	92	95	97
	PW(0.8)	50	94	100 100 100 100 100 100 99	89
	PW(1.2)	1	20	92	96	94	90	84	78	40	16
	PW(1.4)	0	5	72	80	75	65	54	45	15	8
	PW(2.0)	0	0	19	30	30	28	29	32	50	66
	PW(3.0)	12	19	59	78	87	91	94	95	99 100
	LIFR(1)	0	0	20	39	45	45	43	42	29	17
	LIFR(2)	0	1	41	64	69	69	67	65	49	31
	LIFR(4)	0	2	64	84	87	87	86	85	71	49
	LIFR(6)	0	4	75	90	93	93	92	91	81	60
	LIFR(10)	0	8	84	95	97	97	96	96	89	72
	HN	0	0	28	49	54	53	51	50	34	21
	LN (0.7)	4	5	43	75	88	94	96	98 100 100
	LN (0.8)	10	10	20	37	51	60	67	73	89	95
	LN (1.0)	36	37	36	34	29	25	22	21	22	29
	LN (1.5)	91	92	94	95	95	95	94	94	91	84
	HC	95	95	96	96	95	94	94	93	88	80
	JS(0.5)	78	79	82	82	81	80	79	78	70	61
	JS(1.0)	98	98	99	99	99	99	99	99	99	98
	U[0, 1]	9	56	99	100 100 99	98	97	79	49
	G(0.4)	80	83	93	97	98	99	99	99 100 100
	G(0.6)	38	40	51	58	63	67	69	72	78	81
	G(0.8)	13	13	16	17	17	18	18	19	23	27
	G(1.4)	1	1	7	21	28	32	33	35	35	30
	G(1.6)	0	0	18	40	51	56	59	62	62	56
	G(1.8)	0	0	33	61	72	77	79	81	82	77
	G(2.0)	0	1	49	78	86	90	91	93	93	90
	G(2.4)	0	5	79	95	98	99	99	99	99	99
	G(3.0)	1	20	97	100 100 100 100 100 100 100

TABLE 5

 5 Percentage of 10000 Monte Carlo samples declared significant; level α = 0.05; sample size n = 100

	Distribution				Parameters a
		-0.99 -0.9 -0.5 0.1 0.5	1	1.5	2	5	10
	W(0.6)	98	99	100 100 100 100 100 100 100 100
	W(1.2)	0	0	23	44	51	54	56	56	52	44
	W(1.4)	0	8	84	96	98	98	98	98	97	93
	W(1.6)	8	53	100 100 100 100 100 100 100 100
	χ 2 1	78	83	95	98	99	99 100 100 100 100
	PW(0.8)	100	100	100 100 100 100 100 100 100 100
	PW(1.2)	51	94	100 100 100 100 100 98	72	30
	PW(1.4)	14	66	100 100 99	96	90	82	28	10
	PW(2.0)	0	6	76	81	75	69	66	65	77	89
	PW(3.0)	23	47	95	99 100 100 100 100 100 100
	LIFR(1)	0	2	57	76	78	77	76	73	56	37
	LIFR(2)	0	10	85	95	95	95	95	93	82	61
	LIFR(4)	2	31	98	100 100 100 100 99	96	84
	LIFR(6)	6	47	99	100 100 100 100 100 99	91
	LIFR(10)	14	66	100 100 100 100 100 100 100 97
	HN (0, 1)	0	4	70	85	87	86	84	82	64	43
	LN (0.7)	6	12	83	99 100 100 100 100 100 100
	LN (0.8)	15	16	39	71	85	92	96	97 100 100
	LN (1.0)	53	54	56	53	47	44	42	41	45	61
	LN (1.5)	99	99	100 100 100 100 100 100 99	98
	HC	100	100	100 100 100 100 100 100 99	96
	JS(0.5)	94	95	97	97	97	97	96	96	92	85
	JS(1.0)	100	100	100 100 100 100 100 100 100 100
	U[0, 1]	94	100	100 100 100 100 100 100 98	84
	G(0.4)	96	98	100 100 100 100 100 100 100 100
	G(0.6)	52	58	77	86	90	92	94	95	97	97
	G(0.8)	16	18	24	27	29	31	34	35	40	43
	G(1.4)	1	1	22	43	54	59	62	64	66	61
	G(1.6)	0	1	49	76	84	88	90	91	93	91
	G(1.8)	0	5	75	93	97	98	98	99	99	99
	G(2.0)	1	13	91	99	99 100 100 100 100 100
	G(2.4)	6	45	99	100 100 100 100 100 100 100
	G(3.0)	39	90	100 100 100 100 100 100 100 100

Acknowledgement: The authors are indebted to Nora Gürtler for assistance in computation.


Dip Daniele De Martini
email: demartin@eco.unipmn.it

Semeq

Reproducibility probability estimation for testing statistical hypotheses

Keywords: Noncentrality parameter estimation, plug-in power estimation, statistical tests definitions

The application of statistical tests often generates problems, because of their randomness. Here, stemming from the variability of statistical tests, the reproducibility probability of a statistically significant result (RP ) is estimated first. Then, statistical tests are defined on the basis of RP estimation. Lower bounds for RP , and then pointwise RP estimates are used. It is emphasized that the threshold for the RP pointwise estimator to define statistical tests turns out to be 1/2. Hence, taking the perspective of RP estimation, one can compute statistical tests and have suitable interpretations of their results.

Introduction

The applications of statistical tests in various research fields often generates problems, because the results obtained through statistical hypothesis testing have been, and are, often confuted in further studies, for example in the biomedical field. Indeed, [START_REF] Goodman | A comment on replication, p-values and evidence[END_REF] showed that the estimated probability of repeating a statistically significant result is substantially lower than expected, and he complained that the p-value might lead to too optimistic interpretations of the test's result.

Moreover, many institutional drug agencies, which use statistical tests as experimental proofs for evaluating the efficacy and safety of new drugs and medical systems in general, are aware of the randomness of statistical test results, and ask for the reproduction of statistically significant confirmatory trials (e.g. the United States FDA). [START_REF] Shao | Reproducibility Probability in Clinical Trials[END_REF] estimated the reproducibility probability of a statistically significant result (RP ) for the t-test, and argued that a single clinical trial is sufficient if the statistically significant result of the first experiment is evaluated to be highly reproducible.

On the theoretical side, statisticians seldom consider the Bernoulli randomness nature of statistical tests, and therefore their variability.

Here, it is shown that RP estimation can be used not only to evaluate the results of statistical tests, but even to define statistical tests themselves. Then, taking the perspective of RP estimation, one can compute statistical tests and have suitable interpretations of their results.

In Section 2 the theoretical framework adopted is introduced, which first considers only onesided alternatives. Then, in Section 3, stemming from the variability of statistical tests, the RP is estimated, and it is shown that statistical tests can be defined on the basis of the lower bound of a one-directional confidence interval for the RP . Moreover, even the pointwise RP estimator can be used for the definition of statistical tests, providing a suitable interpretation of test's results. Some examples, which cover many widely used testing procedures, are shown in Section 4. Finally, Section 5 contains the conclusions. Furthermore, the perspective of RP estimation is broadened to two-sided alternatives in the Appendix, in order to present a fluent exposition.

Theoretical Framework

In this section, the theoretical framework adopted in order to present the perspective of RP estimation is introduced. Here, only one-sided alternatives are considered. The results on two-sided alternatives are given in the Appendix.

Let X n ∈ S be a random sample of size n, with distribution F ∈ F , where F is a family of distributions over S. The statistical hypotheses, which should be of scientific interest according with [START_REF] Spiegelhalter | A predictive approach to selecting the size of a clinical trial, based on subjective clinical opinion[END_REF], are:

H 0 : F ∈ F 0 , H 1 : F ∈ F\F 0 . Then, let T n = T (X n )
be the test statistic with parametric continuous distribution G n,λn , where λ n = L(n, F ) ∈ R is the parameter of the distribution of T n . It is assumed that the following two conditions hold:

I) G n,λn is stochastically strictly increasing on λ n , that is G n,λ (t) > G n,λ (t) ∀t, if λ < λ (1)
II) sup H0 {λ n } = 0, and G n,0 is known. Now, let α ∈ (0, 1) be the fixed type I error. Then, (II) implies that t n,1-α = G -1 n,0 (1 -α) is the critical value of the test, which is:

Φ α (X n )=        1 iff T n > t n,1-α 0 iff T n ≤ t n,1-α
The power π(n) of the test Φ α then is:

π(n) = P F (Φ α (X n ) = 1) = P F (T n > t n,1-α ) = 1 -G n,λn (t n,1-α ) (2)
Since G n,λn is stochastically strictly increasing (I), it is easy obtaining that Φ α is strictly unbiased (see also Lehmann, 1997, p.137)

, i.e. π(n) ≤ α iff F ∈ F 0 , and π(n) > α iff F ∈ F\F 0 . So, it is
worth noting that statistical hypotheses can be redefined as follows:

H 0 : π(n) ≤ α, H 1 : π(n) > α.
Consequently, confidence intervals for the power might be used to test statistical hypotheses.

RP estimation for testing one-sided alternatives

In this section, the perspective of RP estimation to define statistical tests is introduced. In practice, the RP is estimated through confidence intervals, at first providing lower bounds for the noncentrality parameter, and then by applying the plug-in principle to the power. Finally, statistical tests are defined on the bases of RP estimates.

In the literature, the pointwise estimate of the power of the test is called the Reproducibility Probability significant results (see for example Goodman, 1992, Shao and[START_REF] Shao | Reproducibility Probability in Clinical Trials[END_REF]. Here, in order to avoid confusion between unknown probabilities and their estimates, and between the power in (2) and the power function (see Lehmann, 1997, p.69), the following definition is introduced.

Definition 1. The Reproducibility Probability of a statistically significant result (RP) of the test Φ α (X n ) is the power π(n).

Lemma 1. Lower bound for the noncentrality parameter. Under condition (I), let λn,γ be the solution of the equation G n, λn,γ (T n ) = 1 -γ. Then, λn,γ is the (1 -γ)-lower bound for λ n , that is:

P F ( λn,γ ≤ λ n ) = 1 -γ (3) Proof. Exploit (1) and the uniform distribution on [0, 1] of G n,λn (T n ).
Theorem 1. Lower bound for the RP. Under conditions (I) and (II), πn,γ (n

) = 1-G n, λn,γ (t n,1-α )
is the (1 -γ)-lower bound for the RP, that is:

P F (π n,γ (n) ≤ π(n)) = 1 -γ
Proof. Make use of Lemma (1), of the plug-in principle with (2), and of (1).

Corollary 1.

(1 -α)-lower bound for the RP for testing one-sided alternatives. Under conditions (I) and (II) the following result is obtained:

Φ α (X n )=        1 iff πn,α (n) > α 0 iff πn,α (n) ≤ α
Proof. As a direct consequence of Theorem 1 and (3

): πn,α (n) = α ⇔ T n = t n,1-α . Moreover, using (1): πn,α (n) > α ⇔ T n > t n,1-α .
Remark 1. From the mathematical point of view, the statistical test is a Bernoulli random variable whose parameter is its RP . Thus, the RP is the only parameter of the testing problem.

Consequently, it appears natural that RP estimation can be used to test statistical hypotheses.

It is shown here below that even the pointwise estimator of the RP can be used for testing.

To estimate the RP pointwise, the (50%)-lower bound for the RP is considered, that is πn (n) = πn,50% (n).

Corollary 2. (50%)-lower bound for the RP for testing one-sided alternatives. Under conditions (I) and (II) the following result is obtained:

Φ α (X n )=        1 iff πn (n) > 1/2 0 iff πn (n) ≤ 1/2
Proof. Analogous to Cor.1.

Remark 2. In practice, when the probability to reject H 0 is estimated to be greater than that of accept H 0 , then H 0 is rejected.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Examples

In this Section some examples of widely used testing procedures are provided, which refer to test statistics with Gaussian, t, F , and χ 2 distributions. For those examples, conditions (I) and

(II) are satisfied, and so the results of Sec.3 hold.

Test statistics referring to the Normal distribution

Let us consider a testing situation where the distribution G n,λn of the test statistic T n is a

Normal distribution, with mean λ n and unitary variance (i.e. φ λn,1 ). Then, condition (I) is fulfilled. Moreover, consider that sup H0 {λ n } = 0, which is encountered, for example, testing the superiority of a mean with respect to another one. Then, G n,0 = φ 0,1 is known, and even condition (II) is satisfied.

The critical value of the test is

φ -1 0,1 (1 -α) = z 1-α , and the test is Φ α (X n ) = 1 ⇔ T n > z 1-α .
Hence, the power is π

(n) = 1 -φ λn,1 (z 1-α ) = 1 -φ 0,1 (z 1-α -λ n )
, and according to (3), the

(1 -γ)-lower bound of λ n (i.e. λn,γ ) results λn,γ = T n -z 1-γ .

Thus, the 1 -α lower bound for the power turns out to be πn,α (n) = 1 -φ 0,1 (2z 1-α -T n ), and from Cor.1 the statistical test is:

Φ α (X n ) = 1 ⇔ 1 -φ 0,1 (2z 1-α -T n ) > α.
As far as the pointwise RP estimation is concerned, the estimate of λ n turns out to be the test statistic, that is λn,50% = T n , and consequently πn

(n) = 1 -φ Tn,1 (z 1-α ) = 1 -φ 0,1 (z 1-α -T n ). Hence, from Cor.2: Φ α (X n ) = 1 ⇔ 1 -φ 0,1 (z 1-α -T n ) > 1/2.

Test statistics referring to the t distribution

Here, assume that T n has a t distribution with ν(n) degrees of freedom and noncentrality parameter λ n . It is then easy to verify condition (I). Now, consider the problem of testing the positive values of the regression parameters in the linear model, where the errors are Normally distributed (see for example [START_REF] Searle | Linear models[END_REF]. Then, the test statistic has a t distribution, sup H0 {λ n } = 0, and G n,0 = t ν(n) , so that even condition (II) is satisfied. Thus, the results of Sec.3 hold, with particular emphasis to Cor.2.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Discussion

In Sec.4 it is shown that RP estimation can be applied to many statistical tests. In particu-

(T n ∼ χ 2 ).
RP estimation is also useful for experimental planning, and in particular for sample size estimation. Indeed, lower bounds for RP (i.e., lower bounds for the power) can be used for conservative sample size estimation, once the alternative hypothesis is of scientific relevance, in accordance with [START_REF] Spiegelhalter | A predictive approach to selecting the size of a clinical trial, based on subjective clinical opinion[END_REF].

Let m(1 -β) = min{m s.t. π(m) > 1 -β} be the ideal sample size for the experiment, where 1 -β is the targeted power level, and π(m) is defined in (2). Now, assume that a pilot sample of size n is given, and that λm,γ,n is the lower bound for the noncentrality parameter λ m , i.e.,

P ( λm,γ,n ≤ λ m ) = 1 -γ. Then, the (1 -γ)-conservative sample size estimator is m γ,n (1 -β) = min{m s.t. πm,γ,n (m) > 1 -β}, where πm,γ,n (m) = 1 -G m, λm,γ,n (t m,1-α )
, according with the definition of πm,γ (m) in Theorem 1. By adopting m γ,n (1 -β), one has a probability of γ of planning an underpowered experiment. Often, the pointwise estimator of the sample size is used, that is m 50%,n (1 -β). In order to reduce to, for example, 1/4 the probability of an underpowered experiment, γ = 25% is needed. Recently, [START_REF] Wang | Adapting the sample size planning of a phase III trial based on phase II data[END_REF] have argued the use of conservative sample size estimation with γ = 15.9%, i.e., the one-standard error conservative approach.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

As far as concerns criteria for a single clinical trial, in some circumstances an RP estimate of 90% is considered high enough to justify its use for regulatory approval. Nevertheless, [START_REF] Shao | Reproducibility Probability in Clinical Trials[END_REF] emphasized the variability of the pointwise estimate, and suggest adopting the usual lower bounds for RP (e.g., πn,5% (n)), to be compared with standard thresholds of the power (i.e. 80%, 90%).

As a numerical example, let us consider the results of a test statistic referring to a Gaussian distribution, with α = 2.5%: if T n = 3.24 is observed, then the RP pointwise estimate is πn (n) 90%; but the 95% lower bound of RP , i.e., πn,5% (n), is just 48%; moreover, T n = 4.13 is needed to have πn,5% (n) 80%, and T n = 4.57 to obtain πn,5% (n) 90%.

From the results in Sec.3, since the test is statistically significant if πn,α (n) > α, or, equivalently, if πn (n) > 1/2, it follows that a conservative result on the statistical significance is πn,α (n) > 1/2.

Continuing the numerical example, the criterion proposed here above is fulfilled when πn,2.5% (n) > 1/2, which corresponds to T n > 3.92. For the Gaussian distributed test statistics, this turns out to be T n > 2z 1-α . In our opinion, this criterion might be used as a sufficient statistical condition for a single clinical trial.

Conclusions

Statistical tests contain randomness. Consequently, the estimates of the reproducibility probability of a statistically significant test (RP ) are measures of great interest for the analysis of test results.

From the mathematical point of view, the test is a Bernoulli random variable. The RP is the parameter of the statistical test, and so it is the only parameter to estimate. Indeed, in this paper it is shown that RP confidence intervals can be used, under mild conditions, to compute statistical tests.

Hence, on one hand statistical tests can be computed through RP estimation, and on the other hand RP estimates provide suitable interpretations of the results.

In particular, recall that the threshold for defining statistical tests on the basis of the pointwise estimator of the RP turned out to be 1/2. This means that the null hypothesis H 0 is rejected when it is estimated there is more inclination to reject it than to accept it. In our opinion, this result might bring Frequentists and Bayesians closer to one another.

Appendix. RP estimation for testing two-sided alternatives

In order to provide a fluent exposition of the perspective based on RP estimation for testing statistical hypotheses, in this paper only one-sided alternatives are considered. This Appendix broadens the results of Sec.3 on RP estimation to statistical tests with two-sided alternatives.

Assume that the following two conditions hold: AI) G n,λn is stochastically decreasing on λ n when λ < 0, and stochastically increasing when λ > 0;

AII) under H 0 : λ n = 0, and G n,0 is known.

If α ∈ (0, 0.5), then t n,α/2 = G -1 n,0 (α/2) and t n,1-α/2 = G -1 n,0 (1 -α/2) are the critical values of the test, which is:

Ψ α (X n )=        1 iff T n < t n,α/2 or T n > t n,1-α/2 0 iff t n,α/2 ≤ T n ≤ t n,1-α/2
Then, the power of the test Ψ α is:

π(n) = P F (Ψ α (X n ) = 1) = P F (T n < t n,α/2 ) + P F (T n > t n,1-α/2 )
Under the alternative hypothesis, the power is the sum of the probability of rejecting H 0 choosing the wrong direction, that is π III (n) = P (|T n | < -t n,α/2 ) and is called the type III error [START_REF] Harter | Error Rates and Sample Sizes for Range Tests in Multiple Comparisons[END_REF], and of the probability of rejecting choosing the right one, that is named here the Good Power and is π g (n) = P (|T n | > t n,1-α/2 ).

  lar, RP estimates can be computed for tests for comparing proportions (whose large sample test statistic is Gaussian); for statistical tests for means (whose test statistic T n is t distributed, with two-sample parallel designs as well as with crossover designs); for many nonparametric rank-tests, e.g., the Wilcoxon rank-sum and signed-rank tests (whose test statistic is approximately Gaussian); for tests for comparing more than two means, i.e., multiple-sample one-way ANOVA, and for tests for comparing variabilities (T

n ∼ F ); for tests for comparing time-to-event data, e.g., the log-rank test (which refers to a Gaussian distribution), and for the test on Cox's proportional hazards model

ACCEPTED MANUSCRIPT

Test statistics referring to the F distribution

When T n has a singly non central F distribution (see Stuart and Ord, 1991, Vol.2, Par.23.31) condition (I) is fulfilled.

Consider now the General Linear Model (see for example Stuart and Ord, 1991, Ch.23), where the statistic T n to test general linear hypotheses has a non central F distribution. Under H 0 the distribution of T n reduces to the central F , so that condition (II) is satisfied. Consequently, even in these very widely used tests, the results of Sec.3 hold.

Test statistics referring to the Chi-square distribution

Assume that T n has a Chi-square distribution, with ν(n) degrees of freedom and noncentrality parameter λ n , say χ 2 ν(n),λn = G n,λn . Then, condition (I) on the stochastic ordering is satisfied.

Now, let F belong to the exponential family, and, in the context of the Generalized Linear Model (see [START_REF] Mccullagh | Generalized Linear Models[END_REF], consider the problem of comparing two models specified under H 0 and H 1 respectively. Then, the test statistic T n , which is the difference between the log-likelihood ratio statistic of the two models, is distributed as a Chi-square distribution (in some cases only asymptotically) and under the null hypothesis as a central Chi-square, so that condition (II) is satisfied.

As a further example, consider the problem of testing for independence in a contingency table

(see for example [START_REF] Agresti | Categorical Data Analysis[END_REF]. Then the test statistic has, at least approximately, a Chi-square distribution, which is central under H 0 , and then condition (II) is also fulfilled.

Thus, in both situations the results of the previous Section are available, and mainly to Cor.2.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Now, note that under H 0 : π g (n) = π III (n) = α/2; while under the alternative:

and so statistical hypotheses can redefined as follows:

Definition A1. The Reproducibility Probability of a statistically significant result (RP ) of the test Ψ α (X n ) is the good power π g (n).

Then, analog as to Sec.3, if conditions (AI) and (AII) hold, RP estimation can be used to test statistical hypotheses.

Corollary A1.

(1 -α/2)-lower bound for the RP for testing two-sided alternatives. Under conditions (AI) and (AII) the following result is obtained:

Here, the pointwise estimator of the RP is πg n (n) = πg n,50% (n).

Corollary A2. (50%)-lower bound for the RP for testing two-sided alternatives. Under conditions (AI) and (AII) the following result is obtained: