A goodness-of-fit test for the multivariate Poisson distribution

Statistics & Operations Research Transactions SORT 40 (1) January-June 2016, 113-138 Statistics & Operations Research Transactions © Institut d’Estad´ıstica de Catalunya [email protected] ISSN: 1696-2281 eISSN: 2013-8830 www.idescat.cat/sort/ A goodness-of-fit test for the multivariate Poisson distribution F. Novoa-Mu˜noz1,∗and M.D. Jim´enez-Gamero2 Abstract Bivariate count data arise in several different disciplines and the bivariate Poisson distribution is commonly used to model them. This paper proposes and studies a computationally convenient goodness-of-fit test for this distribution, which is based on an empirical counterpart of a system of equations. The test is consistent against fixed alternatives. The null distribution of the test can be consistently approximated by a parametric bootstrap and by a weighted bootstrap. The goodness of these bootstrap estimators and the power for finite sample sizes are numerically studied. It is shown that the proposed test can be naturally extended to the multivariate Poisson distribution. MSC: 62F03; 62F05; 62F12; 62F40. Keywords: Bivariate Poisson distribution, goodness-of-fit, empirical probability generating function, parametric bootstrap, weighted bootstrap, multivariate Poisson distribution. 1. Introduction Univariate count data appear in many real life situations and the univariate Poisson distribution is frequently used to model this kind of data (see for example Haight, 1967; Johnson and Kotz 1969; Sahai and Khurshid, 1993). G¨urtler and Henze (2000) present a wide variety of procedures for testing goodness-of-fit (gof) for the univariate Poisson distribution. In practice, bivariate count data appear in different areas of knowledge and the bivariate Poisson distribution (BPD), being a generalization of the Poisson distribution, plays a key role in modelling them, provided that such data present a positive correlation. ∗Corresponding author: [email protected] 1Departamento de Estad´ıstica, Universidad del B´ıo-B´ıo, Chile. 2Departamento de Estad´ıstica e I.O., Universidad de Sevilla, Spain. Received: January 2015 Accepted: April 2016 114 A goodness-of-fit test for the multivariate Poisson distribution Different authors have given a definition for the BPD (see for example Kocherlakota and Kocherlakota, 1992). In this article we will work with the one that has received more attention (see for example Holgate, 1964; Johnson, Kotz and Balakrishnan, 1997). Let X1=Y1+Y3and X2=Y2+Y3, where Y1,Y2and Y3are mutually independent Poisson random variables with means θ′ 1=θ1−θ3>0, θ′ 2=θ2−θ3>0 and θ3>0, respectively. The joint distribution of the vector (X1,X2)is called BPD with parameter θ= (θ1,θ2,θ3),(X1,X2)∼BP(θ)for short. In the statistical literature on gof tests for the BPD, which is not so rich as in the univariate case, we found the following: the tests given by Crockett (1979), Loukas and Kemp (1986), Rayner and Best (1995) – these three tests are not consistent against all fixed alternatives – and, more recently, the tests in Novoa-Mu˜noz and Jim´enez-Gamero (2014) (hereafter abbreviated to NJ). The two tests in NJ are consistent against all fixed alternatives. The results in Janssen (2000) assert that the global power function of any nonparametric test is flat on balls of alternatives except for alternatives coming from a finite dimensional subspace. Because of this reason, it is interesting to propose new gof tests able to detect different sets of alternatives. This paper presents a consistent gof test for the BPD. It is based on the following: since the probability generating function (pgf) of the BPD is the unique pgf satisfying certain system of partial differential equations, and the empirical probability generating function (epgf) consistently estimates the pgf, the epgf should approximately satisfy such system. The proposed test statistic is a function of the coefficients of the polynomials of an empirical version of that system. The asymptotic behaviour of the proposed test under alternatives is shared with the ones in NJ. An advantage of the test proposed in this paper over those in NJ is that its application does not entail the choice of a weight function, which is rather arbitrary. The null distribution of the test statistic can be consistently approximated by a parametric bootstrap as well as by means of a weighted bootstrap. The finite sample performance of the proposed test is investigated by means of a simulation study, where the goodness of the proposed approximations is numerically studied and the test is compared, in terms of power, to the tests cited above. The numerical power study reveals that, as expected from the results in Janssen (2000), there is no test yielding the highest power against all considered alternatives. In most cases, the power of the proposed test is quite close to the highest one; in other cases, the proposed test is the most powerful. In addition, from a computational point of view, the test proposed in this paper is more efficient than its competitors. The work is organized as follows. Section 2 introduces the test statistic and derives its asymptotic null distribution. Since the asymptotic null distribution does not provide a useful means of approximating the null distribution of the test statistic, Section 3 stud- F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 115 ies two bootstrap estimators. Specifically, it is shown that the parametric bootstrap and a conveniently defined weighted bootstrap estimators produce consistent null distribution estimators. This Section also studies the power of the resulting tests against fixed alternatives. Section 4 deals with the practical implementation of the bootstrap null distribution estimators as well as other related issues. Section 5 reports a summary of the results of a simulation study carried out to examine the finite sample performance of the tests and to compare them with the existing ones. All stated results are valid for θ3>0. Section 6 deals with the case θ3=0. Section 7 shows how the proposed technique can be applied to the general multivariate case. All proofs are relegated to the last section. Hereinafter we shall use the following notation: all vectors are row vectors and vT is the transposed of the row vector v; for any vector v,vkdenotes its kth coordinate, and kvkits Euclidean norm; N0={0,1,2,3,...};I{A}denotes the indicator function of the set A;Pθdenotes the probability law of the BPD with parameter θ;Pdenotes the probability law of the data; Eθdenotes expectation with respect to the probability function Pθ;Edenotes expectation with respect to the true probability function of the data; P∗denote the probability law, given the data; all limits in this work are taken as n→ ∞;L −→ denotes convergence in distribution; P −→ denotes convergence in probability; a.s. −→ denotes almost sure (a.s.) convergence; for any function h:S⊂Rm→R, for some fixed m∈N, we will denote Da1···amh(u) = ∂ k ∂ ua1 1··· ∂ uam m h(u), ∀a1,...,am∈N0such that k=a1+···+am. 2. The test statistic and its asymptotic null distribution Let X1= (X11,X12),X2= (X21,X22),...,Xn= (Xn1,Xn2)be independent identically distributed (iid) from a random vector X= (X1,X2)∈N2 0. Based on the sample X1,X2,...,Xn, the objective is to test the hypothesis H0:(X1,X2)∼BP(θ1,θ2,θ3),for some (θ1,θ2,θ3)∈Θ, against the alternative H1:(X1,X2)≁BP(θ1,θ2,θ3),∀(θ1,θ2,θ3)∈Θ, 116 A goodness-of-fit test for the multivariate Poisson distribution where Θ=(θ1,θ2,θ3)∈R3:θ1>θ3,θ2>θ3,θ3>0. Since the distribution of a random vector X= (X1,X2)∈N2 0is determined by its pgf g(u) = EuX1 1uX2 2,u= (u1,u2)∈[0,1]2, and the joint pgf of a random vector X∼BP(θ)is g(u;θ) = Eθ(uX1 1uX2 2) = expθ1(u1−1) +θ2(u2−1)+θ3(u1−1)(u2−1),(1) testing H0vs H1is equivalent to testing H0:g(u) = g(u;θ),∀u∈[0,1]2,for some (θ1,θ2,θ3)∈Θ, versus H1:g(u)6=g(u;θ),for some u∈[0,1]2,∀(θ1,θ2,θ3)∈Θ. Proposition 2 in NJ shows that g(u1,u2;θ)is the only pgf in G2={g:[0,1]2→R, such that gis a pgf and ∂ ∂ u1g(u1,u2)and ∂ ∂ u2g(u1,u2)exist ∀(u1,u2)∈[0,1]2}satisfying the following system, Di(u;θ) = 0,i=1,2,∀u∈[0,1]2, where D1(u;θ) = ∂ ∂ u1 g(u1,u2)−θ1+θ3(u2−1)g(u1,u2), D2(u;θ) = ∂ ∂ u2 g(u1,u2)−θ2+θ3(u1−1)g(u1,u2). Now we consider the following empirical versions of the functions Di(u;θ),i=1,2, D1n(u;ˆ θ) = ∂ ∂ u1 gn(u1,u2)−ˆ θ1+ˆ θ3(u2−1)gn(u1,u2), D2n(u;ˆ θ) = ∂ ∂ u2 gn(u1,u2)−ˆ θ2+ˆ θ3(u1−1)gn(u1,u2), where ˆ θ= ( ˆ θ1,ˆ θ2,ˆ θ3)is a consistent estimator of θand gn(u1,u2)is the epgf associated to the data, gn(u1,u2) = 1 n n X i=1 uXi1 1uXi2 2. F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 117 Proposition 1 in NJ shows that g(u)and its derivatives can be consistently estimated by the epgf and the derivatives of the epgf, respectively. Thus, if H0is true then D1n(u;ˆ θ) and D2n(u;ˆ θ)should be close to 0, ∀u∈[0,1]2. This proximity to 0 can be interpreted in several ways. For example, NJ interpreted this proximity as Sn,w(ˆ θ) = nZ{D1n(u;ˆ θ)2+D2n(u;ˆ θ)2}w(u)du ≈0,(2) where w(u)is a non-negative function on [0,1]2. Here we present another interpretation, reasoning as in Nakamura and P´erez-Abreu (1993) for the univariate case. With this aim, observe that Din(u;ˆ θ) = X r≥0X s≥0 di(r,s;ˆ θ)ur 1us 2,i=1,2,(3) where d1(r,s;ˆ θ) = (r+1)pn(r+1,s)−(ˆ θ1−ˆ θ3)pn(r,s)−ˆ θ3pn(r,s−1), d2(r,s;ˆ θ) = (s+1)pn(r,s+1)−(ˆ θ2−ˆ θ3)pn(r,s)−ˆ θ3pn(r−1,s), and pn(r,s) = 1 n n X k=1 I(Xk1=r,Xk2=s) is the relative frequency of the pair (r,s). Thus, Din(u;ˆ θ) = 0, ∀u∈[0,1]2,i=1,2, if and only if the coefficient of ur 1us 2in the right hand side of (3) is null, ∀r,s≥0, i=1,2. This leads us to consider the following statistic for testing H0, Wn(ˆ θ) = X r≥0X s≥0{d1(r,s;ˆ θ)2+d2(r,s;ˆ θ)2}= M X r,s=0{d1(r,s;ˆ θ)2+d2(r,s;ˆ θ)2},(4) where M=max{X(n)1,X(n)2},X(n)k=max1≤i≤nXik,k=1,2. Taking into account that dk(r,s;ˆ θ) = 1 n n X i=1 φkrs(Xi;ˆ θ),k=1,2, 118 A goodness-of-fit test for the multivariate Poisson distribution with φ1rs(x;θ) = (r+1)I(x1=r+1,x2=s)−(θ1−θ3)I(x1=r,x2=s)−θ3I(x1=r,x2=s−1), φ2rs(x;θ) = (s+1)I(x1=r,x2=s+1)−(θ2−θ3)I(x1=r,x2=s)−θ3I(x1=r−1,x2=s), where x= (x1,x2), the statistic Wn(ˆ θ)can be expressed as follows, Wn(ˆ θ) = 1 n2 n X i,j=1 h(Xi,Xj;ˆ θ), with h(x,y;θ) = h1(x,y;θ) +h2(x,y;θ), h1(x,y;θ) = Pr≥0Ps≥0φ1rs(x;θ)φ1rs(y;θ) ={x2 1+(θ1−θ3)2+θ2 3}I(x1=y1,x2=y2)−(θ1−θ3)x1I(x1=y1+1,x2=y2) −θ3x1I(x1=y1+1,x2=y2+1)+ (θ1−θ3)θ3I(x1=y1,x2=y2+1) −(θ1−θ3)y1I(y1=x1+1,y2=x2)−θ3y1I(y1=x1+1,y2=x2+1) +(θ1−θ3)θ3I(y1=x1,y2=x2+1), h2(x,y;θ) = Pr≥0Ps≥0φ2rs(x;θ)φ2rs(y;θ) ={x2 2+(θ2−θ3)2+θ2 3}I(x1=y1,x2=y2)−(θ2−θ3)x2I(x1=y1,x2=y2+1) −θ3x2I(x1=y1+1,x2=y2+1)+ (θ2−θ3)θ3I(x1=y1+1,x2=y2) −(θ2−θ3)y2I(y1=x1,y2=x2+1)−θ3y2I(y1=x1+1,y2=x2+1) +(θ2−θ3)θ3I(y1=x1+1,y2=x2), where x= (x1,x2)and y= (y1,y2). In order to give a sound justification of Wn(ˆ θ)as a test statistic for testing H0we next derive its a.s. limit. Theorem 1 Let X1,X2,...,Xnbe iid from X= (X1,X2)∈N2 0with E(X2 k)<∞, k =1,2. Let p(r,s) = P(X1=r,X2=s). If ˆ θa.s. −→θ, for some θ∈R3, then Wn(ˆ θ)a.s. −→X r,s≥0a1(r,s;θ)2+a2(r,s;θ)2=η(P;θ), where a1(r,s;θ) = (r+1)p(r+1,s)−(θ1−θ3)p(r,s)−θ3p(r,s−1), a2(r,s;θ) = (s+1)p(r,s+1)−(θ2−θ3)p(r,s)−θ3p(r−1,s). F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 119 Note that η(P;θ)≥0 and, taking into account that Dk(u;θ) = X r≥0X s≥0 ak(r,s;θ)ur 1us 2,k=1,2, it follows that η(P;θ) = 0 if and only if H0is true. Thus, a reasonable test for testing H0 should reject the null hypothesis for large values of Wn(ˆ θ). Now, to determine what are large values we must calculate its null distribution, or at least an approximation to it. We first try to estimate the null distribution of Wn(ˆ θ)by means of its asymptotic null distribution. In order to derive it, it will be assumed that the estimator ˆ θis asymptotically linear, as expressed in the next assumption. Assumption 1 Under H0, if θ= (θ1,θ2,θ3)∈Θdenotes the true parameter value, then √nˆ θ−θ=1 √n n X i=1 ℓ(Xi;θ)+oP(1), where ℓ:N2 0×Θ−→R3is such that Eθ{ℓ(X1;θ)}=0and J(θ) = Eθℓ(X1;θ)Tℓ(X1;θ) <∞. Assumption 1 is not restrictive at all since it is fulfilled by some commonly used estimators such as the moment estimator, the maximum likelihood estimator, the double zero estimator, the even points estimator and the conditional even points estimator (see Kocherlakota and Kocherlakota, 1992, and Papageorgiou and Loukas, 1988). The next result gives the asymptotic null distribution of Wn(ˆ θ). Theorem 2 Let X1,X2,...,Xnbe iid from X= (X1,X2)∼BP(θ1,θ2,θ3). Suppose that Assumption 1 holds. Then nWn(ˆ θ)L −→ X j≥1 λjχ2 1j, where χ2 11,χ2 12,... are independent χ2variates with one degree of freedom and the set {λj}are the non-null eigenvalues of the operator C(θ)defined on the function space {τ:N2 0→R,such that Eθτ2(X)<∞,∀θ∈Θ}, as follows C(θ)τ(x) = Eθ{K(x,X;θ)τ(X)}, with K(x,y;θ) = h(x,y;θ)+ℓ(x;θ)µ(y;θ)T+ℓ(y;θ)µ(x;θ)T+ℓ(x;θ)S(θ)ℓ(y;θ)T ,µ(x;θ) = (µ1(x;θ),µ2(x;θ),µ3(x;θ)), 120 A goodness-of-fit test for the multivariate Poisson distribution µ1(x;θ) = −x1Pθ(x1−1,x2) +θ3Pθ(x1,x2+1) +(θ1−θ3)Pθ(x1,x2), µ2(x;θ) = −x2Pθ(x1,x2−1) +θ3Pθ(x1+1,x2) +(θ2−θ3)Pθ(x1,x2), µ3(x;θ) = −µ1(x;θ)−x1Pθ(x1−1,x2−1) +θ3Pθ(x1,x2) +(θ1−θ3)Pθ(x1,x2−1) −µ2(x;θ)−x2Pθ(x1−1,x2−1) +θ3Pθ(x1,x2) +(θ2−θ3)Pθ(x1−1,x2), S(θ) = Pr,s≥0Srs(θ), Srs(θ) =   a20a(b−a) 0a2a(c−a) a(b−a)a(c−a) (b−a)2+(c−a)2 , a=Pθ(r,s), b =Pθ(r,s−1), c =Pθ(r−1,s), The asymptotic null distribution of Wn(ˆ θ)does not provide a useful approximation to its null distribution since it depends on the unknown true value of θ. Even if θwere known or replaced by an appropriate estimator, to determine the eigenvalues of an operator is a rather hard problem. So, we next study two further ways of approximating it: a parametric bootstrap (PB) estimator and a weighted bootstrap (WB) estimator. 3. Approximating the null distribution 3.1. Parametric bootstrap Let X1,X2,...,Xnbe iid taking values in N2 0such that ˆ θ=ˆ θ(X1,X2,...,Xn)∈Θ. Let X∗ 1,X∗ 2,...,X∗ nbe iid from a population with distribution BPˆ θ, given X1,X2, ...,Xn, and let W∗ n(ˆ θ∗)be the bootstrap version of Wn(ˆ θ)obtained by replacing X1,X2, ...,Xn and ˆ θ=ˆ θX1,X2,...,Xnby X∗ 1,X∗ 2,...,X∗ nand ˆ θ∗=ˆ θX∗ 1,X∗ 2,...,X∗ n, respectively, in the expression of Wn(ˆ θ). To prove that the PB can be used to consistently approximate the null distribution of Wn(ˆ θ), we will assume the following, which is a bit stronger than Assumption 1. Assumption 2 Assumption 1 holds and the functions ℓand J satisfy (1) supϑ∈Θ0Eϑkℓ(X;ϑ)k2I{kℓ(X;ϑ)k>γ}→0, as γ→∞, where Θ0⊆Θis an open neighborhood of θ. (2) ℓ(X;ϑ)and J(ϑ)are continuous as functions of ϑat ϑ=θand J(ϑ)is finite ∀ϑ∈Θ0. F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 121 Theorem 3 Let X1,X2,...,Xnbe iid from X= (X1,X2)∈N2 0. Suppose that Assumption 2 holds and that ˆ θa.s. −→θ, for some θ∈Θ. Then sup x∈RP∗nW ∗ n(ˆ θ∗)≤x−PθnWn(ˆ θ)≤x a.s. −→ 0. Let w∗ n,α=inf{x:P∗(W∗ n(ˆ θ∗)≥x)≤α}be the αupper percentile of the PB distribution of Wn(ˆ θ)and let Wobs be the observed value of the test statistic. From Theorem 3, the test function Ψ∗ PB =(1,if Wn(ˆ θ)≥w∗ n,α, 0,otherwise, or equivalently, the test that rejects H0when p∗=P∗W∗ n(ˆ θ∗)≥Wobs≤α, is asymptotically correct in the sense that Pθ(Ψ∗ PB =1)→α. 3.2. Weighted bootstrap From the proof of Theorem 2, when H0is true, we have that nWn(ˆ θ) = nW1n(θ)+oP(1), where nW1n(θ) = 1 n n X i,j=1 K(Xi,Xj;θ), which converges in law to W0=Pj≥1λjχ2 1j. As observed before, the greatest difficulty with W0is to determine the set {λj}. Nevertheless, Delhing and Mikosch (1994) have shown that the eigenvalues {λj}can be consistently (a.s.) approximated by the eigenvalues of the matrix Hn=1 nK(Xi,Xj;θ)1≤i,j≤n , say ˆ λ1,...,ˆ λn. Therefore, we could approximate the null distribution of nW1n(ˆ θ)(and thus that of nWn(ˆ θ)) through the conditional distribution, given X1,...,Xn, of nW ∗ 1n= n X j=1 ˆ λjχ2 1j. This is tantamount to approximate the null distribution of nW1n(ˆ θ)by means of the conditional distribution, given X1,...,Xn, of 128 A goodness-of-fit test for the multivariate Poisson distribution algorithm 1 with Imhof’s method) were, as expected, quite close. As for raw multipliers versus centered multipliers, a bit better results are obtained when using the centered multipliers. Table 2 displays the fraction of estimated p-values less than or equal to 0.05 and 0.10, which are the estimated type I error probabilities for α=0.05 and 0.10, respectively by using PB and WB with centered multipliers. From the results in this table it can be concluded that both approximations give rise to conservative tests for small sample sizes. As the values of θ1and θ2increase, the tests become more conservative, specially the one based on the WB approximation. For example, when θ1=θ2=50 and ρ=0.25, the sample size required to get empirical levels close to the nominal values is n=4000. For θ1=θ2=50 and ρ=0.75, n=3000 is enough. In general, better results (in the sense of closeness to the nominal values) are obtained for ρ=0.75 than for ρ=0.25. Finally, it is also observed a bit better results when the parameter is estimated by the maximum likelihood estimator. To study the power we repeated the above experiment for samples with size n= 50 from the following alternatives: bivariate binomial distribution BB(m;p1,p2,p3), where p1+p2−p3≤1,p1≥p3,p2≥p3and p3>0; bivariate Hermite distribution BH(µ,σ2;λ1,λ2,λ3), where µ>σ2(λ1+λ2+λ3); bivariate logarithmic series distribution BLS(λ1,λ2,λ3), where 0 <λ1+λ2+λ3<1; bivariate Neyman type A distribution BNTA(λ;λ1,λ2,λ3), where 0 <λ1+λ2+λ3≤1; bivariate Poisson distribution mixtures of the form pBP(θ) + (1−p)BP(λ),0<p<1, denoted by BPP(p;θ,λ); and (X1,X2) with X1=max{Y1,Y3}and X2=|Y1−Y3|(type 1), X1=max{Y2,Y3}and X2=|Y2−Y3| (type 2), X1=max{Y1,Y3}and X2=min{Y2,Y3}(type 3), X1=max{Y2,Y3}and X2= min{Y1,Y3}(type 4), X1=max{Y1,Y3}and X2=max{Y2,Y3}(type 5), where Y1,Y2,Y3 are independent variables taking values in N0whose distribution are binomial B(m;p), negative binomial BN(m;p), Poisson P(λ)and uniform on 1,2,...,m,U(m). The values of the parameters were chosen so that the expectations E(X1)and E(X2)are small for the PB and the WB not to be excessively conservative. In this part of the simulation experiment we only considered the maximum likelihood estimator of the parameter. In addition to the tests proposed in this paper, Ψ∗ PB and Ψ∗ WB, we also considered the tests given in Crockett (1979) (denoted by T), Loukas and Kemp (1986) (denoted by IB), Rayner and Best (1995) (denoted by NIB) and NJ (denoted by Rnand Sn, with weight function w(u) = 1). Table 3 displays the alternatives considered and the estimated power for nominal significance level α=0.05. Looking at this table we conclude that the tests Ψ∗ PB,Ψ∗ WB,Rnand Snare able to detect all considered alternatives while, as expected, the other tests cannot, specially the tests based on IBand NIB. For the alternatives in the first half of Table 3 we see that the powers of the new tests, Rnand Snare quite close; while for the other alternatives the tests proposed in this paper are more powerful than Rnand Sn. We also compared these tests from a computational point of view. From the results in Table 1 we saw that, in this respect, Ψ∗ WB is more efficient than Ψ∗ PB. Since Rn and Snare both based on a PB, for the comparisons to be fair, we compared Ψ∗ PB,Rnand Sn. Table 4 reports the ratio of the average CPU to get a p-value. Clearly, regarding the required computing time, Ψ∗ PB is more efficient than Rnand Sn. F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 129 Table 3: Simulation results for the power. Alternative E(X1)var(X1) E(X1)E(X2)var(X2) E(X2)ρRnSnΨ∗ PB Ψ∗ W B T IBNIB BB(1,0.45,0.02,0.01)0.450 0.550 0.020 0.980 0.014 0.953 0.961 0.944 1.000 0.263 0.000 0.000 BB(1,0.55,0.03,0.02)0.550 0.450 0.030 0.970 0.041 0.998 0.999 0.998 1.000 0.768 0.000 0.000 BB(2,0.71,0.04,0.03)1.420 0.290 0.080 0.960 0.018 0.997 0.993 1.000 1.000 0.999 0.000 0.000 BH(0.99,1,0.66,0.10,0.10)0.752 1.768 0.198 1.202 0.446 0.963 0.985 0.996 0.982 0.747 0.809 0.842 BH(1.40,1,1.00,0.26,0.12)1.568 1.800 0.532 1.271 0.430 0.967 0.983 0.999 0.999 0.795 0.849 0.899 BH(1.50,1,1.00,0.38,0.10)1.650 1.733 0.720 1.320 0.411 0.939 0.971 0.989 0.968 0.771 0.849 0.887 BLS(0.30,0.01,0.11)1.298 0.409 0.380 0.827 0.303 1.000 1.000 1.000 1.000 0.830 0.012 0.006 BLS(0.40,0.01,0.02)1.311 0.426 0.094 0.959 0.039 1.000 1.000 1.000 1.000 0.763 0.008 0.007 BLS(0.50,0.01,0.02)1.465 0.641 0.085 0.083 0.093 1.000 1.000 1.000 1.000 0.446 0.027 0.025 BNTA(0.1,0.01,0.01,0.93)0.094 1.940 0.094 1.940 0.995 0.797 0.793 0.860 0.225 0.565 0.007 0.599 BNTA(0.1,0.01,0.01,0.92)0.093 1.930 0.093 1.930 0.994 0.810 0.810 0.869 0.222 0.580 0.005 0.617 BNTA(0.1,0.01,0.01,0.95)0.096 1.960 0.096 1.960 0.995 0.835 0.833 0.884 0.188 0.604 0.003 0.629 BPP(0.30; (0.2,0.2,0.1),(0.9,0.9,0.5)) 0.690 1.149 0.690 1.149 0.665 0.846 0.844 0.623 0.729 0.548 0.003 0.001 BPP(0.31; (0.2,0.2,0.1),(1.0,1.2,0.9)) 0.752 1.182 0.890 1.240 0.909 0.726 0.652 0.673 0.755 0.527 0.004 0.001 BPP(0.35; (0.2,0.2,0.1),(0.9,0.9,0.6)) 0.655 1.170 0.655 1.170 0.778 0.805 0.799 0.618 0.722 0.516 0.001 0.000 B B P(2,0.1;2,0.1;0.9), type 1 0.983 0.854 0.870 0.935 0.937 0.229 0.486 0.992 0.938 0.054 0.045 0.058 B P P(2,0.1;0.9;0.9), type 1 0.981 0.850 0.867 0.930 0.936 0.224 0.481 0.991 0.926 0.081 0.06 0.058 B U U(2,0.1;2;2), type 1 1.070 0.566 0.941 0.643 0.902 0.066 0.265 1.000 1.000 0.483 0.000 0.000 B B P(2,0.1;2,0.1;2.1), type 2 2.126 0.950 1.953 1.063 0.962 0.329 0.455 0.965 0.773 0.060 0.075 0.158 B B U(2,0.1;2,0.1;2), type 2 1.069 0.566 0.940 0.641 0.901 0.073 0.289 1.000 1.000 0.491 0.000 0.000 B P U(2,0.1;0.1;2), type 2 1.032 0.618 0.966 0.659 0.949 0.095 0.448 1.000 1.000 0.350 0.000 0.001 B BN P(2,0.1;2,0.1;2.1), type 3 2.129 0.950 2.032 0.997 0.951 0.370 0.516 1.000 0.995 0.028 0.080 0.134 P BN B(0.1;2,0.1;2,0.1), type 3 0.281 0.825 0.196 0.906 0.796 0.120 0.350 0.838 0.851 0.007 0.051 0.018 P BN P(0.6;2,0.1;0.6), type 3 0.976 0.725 0.586 1.003 0.644 0.123 0.205 0.848 0.825 0.177 0.066 0.026 BN B B(2,0.1;2,0.1;2,0.15), type 4 0.446 0.695 0.296 0.848 0.750 0.033 0.146 0.855 0.914 0.063 0.006 0.001 BN B P(2,0.1;2,0.1;0.4), type 4 0.537 0.790 0.394 0.993 0.833 0.176 0.355 0.927 0.850 0.085 0.087 0.053 BN B P(2,0.1;2,0.1;0.5), type 4 0.623 0.809 0.493 1.002 0.872 0.180 0.407 0.971 0.899 0.114 0.087 0.053 B U B(2,0.1;1;2,0.15), type 5 0.445 0.700 0.662 0.406 0.428 0.994 0.996 0.982 0.995 0.940 0.000 0.000 B U U(2,0.1;1;1), type 5 0.605 0.427 0.752 0.248 0.462 1.000 1.000 1.000 1.000 1.000 0.000 0.000 B U U(2,0.1;3;3), type 5 1.557 0.741 2.126 0.404 0.597 0.300 0.187 0.953 0.952 0.580 0.000 0.000 130 A goodness-of-fit test for the multivariate Poisson distribution Table 4: Ratio of average CPU time (in seconds). n=30 n=50 n=70 n=100 n=200 n=300 Rn/Ψ∗ PB 73.50 75.71 77.44 80.19 79.69 79.20 Sn/Ψ∗ PB 5.01 11.07 20.20 43.73 145.28 303.92 Table 5: Results for the real data sets. Plants Health Rn,(0,0)0.003 0.000 Rn,(1,0)0.005 0.000 Rn,(0,1)0.010 0.000 Sn,(0,0)0.005 0.002 Sn,(1,0)0.009 0.000 Sn,(0,1)0.011 0.000 Ψ∗ PB 0.049 0.000 ˆ θn(0.64000, 0.94000, 0.19852) (0.30173, 1.21830, 0.12518) To end this section, Ψ∗ PB is applied to two real data sets. The first one were first given and analysed by Holgate (1966), and refers to the number of plants of the species Lacistema aggregatum and Protium guianense in each of 100 contiguous quadrats. Crockett (1979), Loukas and Kemp (1986), Rayner and Best (1995) and NJ tested the data for agreement with the bivariate Poisson model, they all concluded the data were not well modelled by a BPD. The second data set were analysed in Karlis and Tsiamyrtzis (2008), who used two variables, the number of consultations with a doctor or a specialist (X1) and the total number of prescribed and non-prescribed medications used in past 2 days (X2), from the Australian Health survey for 1977–1978. The sample size was quite large (n=5190). These authors assumed that (X1,X2)has a BPD. NJ tested these data sets for agreement with the bivariate Poisson model, concluding that they were not well modelled by a BPD. The p-values obtained by applying the test proposed in this paper to these two real data sets are 0.049 and 0.000, respectively, in agreement with the previous analyses. Table 6: Simulations results for the type I error probabilities when θ3=0. θ1=θ2=1θ1=θ2=3θ1=θ2=10 ML MM ML MM ML MM n5% 10% 5% 10% 5% 10% 5% 10% 5% 10% 5% 10% 100 3.4 7.4 3.5 7.9 3.1 7.3 3.0 7.3 1.0 3.4 0.9 3.4 200 4.2 8.0 4.3 9.1 3.4 8.0 3.3 7.9 2.2 6.6 2.3 6.6 300 4.4 8.7 4.6 9.4 3.7 8.5 3.7 8.5 3.7 7.9 3.7 8.0 F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 131 6. Case θ3=0 θ3=0 θ3=0 The case θ3=0 has been excluded from H0because it is a boundary point. It is wellknown (see, for example Andrews, 1999, Self and Liang, 1987, and the references therein) that in such a case the MLE is not asymptotically normally distributed and thus Assumption 1 is not satisfied. Moreover, Andrews (2000) have proven that the bootstrap does not provides a consistent estimator of the distribution. Therefore, the theory so far developed is not valid for θ3=0. Next we give two possible ways of dealing with this case. A first way consist in applying the method in Feng and McCulloch (1992), which proposed to enlarge the parametric space to θ∈R3, so that negative values for ˆ θ3are allowed. With this approach all required assumptions in our theory are satisfied. The only problem with this solution is how to apply in practice the PB approximation because it implies the generation of samples from a BP(θ1,θ2,θ3)distribution with θ3<0. Nevertheless, the WB approximation can be applied. Table 6 gives the result of a small simulation that studies the goodness of this solution. Observe that the results are quite close to those obtained for θ3>0. Another possible way of dealing with this case is to adapt the alternatives to the usual bootstrap proposed in Andrews (2000). Two of them consists in subsampling, while the other two are based on testing if the parameter is in the boundary. For the later methods we could calculate a confidence interval for θ3and look if it contains 0 by applying, for example, the method in Feng and McCulloch (1992) but, as recognized by the authors, it requires rather large sample sizes. Note that testing for θ3=0 is tantamount to having two independent Poisson variables. Another way of investigating the independence of the marginal distributions is by applying the classical χ2-test. Nevertheless, such test requires the data to be grouped in classes, and the decision could depend on the grouping. In our view, there is a need of a test for independence of variables taking values on N0, which will be the topic of a future research. If it can be reasonably assumed that the variables are independent, then by using Raikov’s theorem (which states that the sum of two independent non-negative random variables has a Poisson distribution if and only if both random variables have the Poisson distribution), testing gof for an independent Poisson model is equivalent to testing gof to the sum of the components to a univariate Poisson model. In the statistical literature there is a variety of test for testing gof to a univariate Poisson model (see, for example, the review in G¨urtler and Henze, 2000). 132 A goodness-of-fit test for the multivariate Poisson distribution 7. The general m m m-variate case This section shows that the proposed test can be extended to the general m-variate case, for any m≥2. Let X1=Y1+Ym+1,X2=Y2+Ym+1, ..., Xm=Ym+Ym+1, where Y1,Y2,...,Ym+1are mutually independent Poisson random variables with means θ′ 1=θ1−θm+1>0,...,θ′ m=θm−θm+1>0 and θm+1>0, respectively. The joint distribution of the vector (X1,X2,...,Xm)is called a m-variate Poisson distribution with parameter θ= (θ1,θ2,...,θm+1)(see Johnson, Kotz and Balakrishnan, 1997). The joint pgf of (X1,X2,...,Xm)is g(u;θ) = exp(m X i=1 θi(ui−1)+ θm+1 m ∏ i=1 ui− m X i=1 ui+m−1!),∀u∈Rm.(7) Now, the objective is to test the hypothesis H0m:(X1,X2,...,Xm)has a m-variate Poisson distribution. In order to extend the proposed test to the general m-variate case we will use the following result in Proposition 3 in NJ which states that g(u;θ)is the only pgf in Gm= {g:[0,1]m→R, such that gis a pgf and ∂ ∂ uig(u1,u2,...,um)exists ∀u∈[0,1]m, 1 ≤i≤ m}satisfying the following system, Di(u;θ) = 0,1≤i≤m,(8) ∀u∈[0,1]m, where Di(u;θ) = ∂ ∂ ui g(u)−(θi+θm+1 ∏ j6=i uj−1!)g(u),1≤i≤m. Let (X1,X2,...,Xm)∈Nm 0be a random vector and let g(u1,u2,...,um) = EuX1 1uX2 2···uXm m its pgf. Then, taking into account that g(u) = X r1,r2,...,rm≥0 ur1 1ur2 2···urm mp(r1,r2,...,rm), where p(r1,r2,...,rm) = P(X1=r1,X2=r2,...,Xm=rm), we can write Di(u;θ) = X r1,r2,...,rm≥0n(ri+1)p(r1,...,ri−1,ri+1,ri+1,...,rm)−(θi−θm+1)p(r1,r2,...,rm) −θm+1p(r1−1,...,ri−1−1,ri,ri+1−1,...,rm−1)our1 1ur2 2···urm m,1≤i≤m. F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 133 Let Din(u;ˆ θ)denote the empirical counterpart of Di(u;θ)obtained by replacing the pgf gby the epgf gnand θby a consistent estimator ˆ θ, 1 ≤i≤m. If H0mis true then the functions Din(u;ˆ θ), 1 ≤i≤m, should be close to 0, ∀u∈[0,1]m. This proximity to zero can be interpreted as we did in Section 2, for the bivariate case. Observe that Din(u;ˆ θ) = X r1,r2,...,rm≥0 di(r1,r2,...,rm;ˆ θ)ur1 1ur2 2···urm m,1≤i≤m, where di(r1,r2,...,rm;ˆ θ) = (ri+1)pn(r1,...,ri−1,ri+1,ri+1,...,rm) −(ˆ θi−ˆ θm+1)pn(r1,r2,...,rm) −ˆ θm+1pn(r1−1,...,ri−1−1,ri,ri+1−1,...,rm−1),1≤i≤m, and pn(r1,r2,...,rm) = 1 nPn k=1I(Xk1=r1,Xk2=r2,...,Xkm =rm)is the relative frequency of (r1,r2,...,rm). Therefore, Din(u;ˆ θ) = 0, ∀u∈[0,1]m, 1 ≤i≤m, if and only if the coefficients of ur1 1ur2 2···urm min the previous expansions are null, ∀r1,r2,...,rm≥0. This leads us to consider the following statistic for testing H0m, Wm,n(ˆ θ) = X r1,r2,...,rm≥0(m X i=1 di(r1,r2,...,rm;ˆ θ)2)= M X r1,r2,...,rm=0(m X i=1 di(r1,r2,...,rm;ˆ θ)2), where M=max{X(n)1,X(n)2,...,X(n)m},X(n)k=max1≤i≤nXik, 1 ≤k≤m. Similar results to those stated in Sections 2 and 3 for the bivariate case can be established for Wm,n(ˆ θ). 8. Proofs Here we give a sketch of the proofs of the results in Sections 2 and 3. A detailed derivation of the results can be obtained from the authors upon request. Proof of Theorem 1 Observe that d1(r,s;ˆ θ) = d1(r,s;θ)−(ˆ θ1−θ1)pn(r,s) +( ˆ θ3−θ3){pn(r,s)−pn(r,s−1)} and X r,s≥0 d1(r,s;θ)2=1 n2X i6=j h1(Xi,Xj;θ)+ 1 n2 n X i=1 h1(Xi,Xi;θ). 134 A goodness-of-fit test for the multivariate Poisson distribution By the SLLN, 1 n n X i=1 h1(Xi,Xi;θ)a.s. −→E  X r,s≥0 φ1rs(X1;θ)2   <∞. By the SLLN for U-statistics (Theorem 5.4 in Serfling, 1980), 1 n2X i6=j h1(Xi,Xj;θ)a.s. −→E{h1(X1,X2;θ)}=X r,s≥0 a1(r,s;θ)2. Therefore, X r,s≥0 d1(r,s;θ)2a.s. −→X r,s≥0 a1(r,s;θ)2. Since pn(r,s)2≤pn(r,s),∀r,s≥0, and Pr,s≥0pn(r,s) = 1, we have (ˆ θ1−θ1)2X r,s≥0 pn(r,s)2≤(ˆ θ1−θ1)2=o(1), and analogously, (ˆ θ3−θ3)2X r,s≥0{pn(r,s)−pn(r,s−1)}2=o(1). Thus, X r,s≥0 d1(r,s;ˆ θ)2a.s. −→X r,s≥0 a1(r,s;θ)2.(9) Following similar steps we get X r,s≥0 d2(r,s;ˆ θ)2a.s. −→X r,s≥0 a2(r,s;θ)2.(10) Finally, the result is obtained from (9) and (10). Proof of Theorem 2 Let us consider the separable Hilbert space of functions H={g: N0→R,so thatkgk2 H=Pr≥0Ps≥0g(r,s)2<∞}. We have that F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 135 √ndk(r,s;ˆ θ) = 1 √n n X i=1 φkrs(Xi;θ)+ √n(ˆ θ−θ)ˆvk(r,s)T,k=1,2, with ˆv1(r,s) = (−pn(r,s),0,pn(r,s)−pn(r,s−1)) and ˆv2(r,s) = (0,−pn(r,s),pn(r,s)− pn(r−1,s)). From Assumption 1 and the SLLN, we get that √ndk(r,s;ˆ θ) = √nd1k(r,s;θ)+Rk(r,s),k=1,2, with d1k(r,s;θ) = 1 n n X i=1φkrs(Xi;θ)+ℓ(Xi;θ)vk(r,s;θ)T,k=1,2, v1(r,s;θ) = (−Pθ(r,s),0,Pθ(r,s)−Pθ(r,s−1)), v2(r,s;θ) = (0,−Pθ(r,s),Pθ(r,s)−Pθ(r−1,s)), and kRkkH=oP(1),k=1,2. From the CLT in Hilbert spaces (see, for example, van der Vaart and Wellner, 1996, pp. 50–51), it follows that k√nd1kk2 H=OP(1),k=1,2, and therefore nWn(ˆ θ) = k√nd1kk2 H+k√nd12k2 H+oP(1). Routine calculations show that k√nd1kk2 H+k√nd12k2 H=1 n n X i,j=1 K(Xi,Xj;θ). The result is achieved by applying Theorem 6.4.1.B in Serfling (1980) to 1 nPn i,j=1K(Xi,Xj;θ). Proof of Theorem 3 Following similar steps to those given in the proof of Theorem 2 but instead of applying the CLT for iid random elements taking values in H, we apply a CLT for triangular arrays, such as Theorem 1.1 in Kundu et al. (2000). Proof of Theorem 4 nW ∗ 2n(ˆ θ)can be expressed as nW∗ 2n(ˆ θ) = W∗ 1+W∗ 2+2W∗ 3+W∗ 4, where W∗ 1=1 n n X i,j=1 K(Xi,Xj;θ)ξiξj, 136 A goodness-of-fit test for the multivariate Poisson distribution W∗ 2=1 n n X i,j=1{h(Xi,Xj;ˆ θ)−h(Xi,Xj;θ)}ξiξj, W∗ 3=1 n n X i,j=1{ˆ ℓ(Xi;ˆ θ)µ(Xj;ˆ θ)T−ℓ1(Xi;θ)µ(Xj;θ)T}ξiξj, W∗ 4=1 n n X i,j=1{ˆ ℓ(Xi;ˆ θ)S(ˆ θ)ˆ ℓ(Xj;ˆ θ)T−ℓ1(Xi;θ)S(θ)ℓ1(Xj;θ)T}ξiξj. From the results in Delhing and Mikosch (1994), sup x|P∗{W∗ 1≤x}−P{W1≤x}| a.s. −→0. Thus, to show the result it suffices to see that W∗ k=oP∗(1)in probability, k=2,3,4. We first deal with W∗ 2. Observe that E∗(W∗2 2)≤M1 n2 n X i,j=1{h(Xi,Xj;ˆ θ)−h(Xi,Xj;θ)}2, for some positive M>0. From the assumptions made, the right-hand side of the above expression is oP(1). Therefore, W∗ 2=oP∗(1)in probability. As for W∗ 3, we have that W∗ 3=W∗ 31W∗⊤ 32 +W∗ 33W∗⊤ 34 , with W∗ 31 =1 √n n X i=1{ˆ ℓ(Xi;ˆ θ)−ℓ1(Xi;θ)}ξi, W∗ 32 =1 √n n X i=1 µ(Xi;ˆ θ)ξi, W∗ 33 =1 √n n X i=1 ℓ1(Xi;θ)ξi, W∗ 34 =1 √n n X i=1{µ(Xi;ˆ θ)−µ(Xi;θ)}ξi. From the assumptions made, E∗(W∗2 31 ) = oP(1),E∗(W∗2 32 )is bounded in probability and E∗(W∗2 33 )is bounded a.s.. Now taking into account that F. Novoa-Mu˜ noz and M.D. Jim´ enez-Gamero 137 ∂ ∂ θ1 Pθ(r,s) = Pθ(r−1,s)−Pθ(r,s), ∂ ∂ θ2 Pθ(r,s) = Pθ(r,s−1)−Pθ(r,s), ∂ ∂ θ3 Pθ(r,s) = Pθ(r−1,s−1)−Pθ(r−1,s)−Pθ(r,s−1)+Pθ(r,s), it follows that sup r,s∈N0|Pˆ θ(r,s)−Pθ(r,s)|≤Mkˆ θ−θk,(11) for some positive M>0. This implies that E∗(W∗2 34 ) = oP(1). Therefore, W∗ 3=oP∗(1)in probability. By using (11) and the assumptions made, it readily follows that W∗ 4=oP∗(1) in probability. This concludes the proof. Acknowledgements The authors thank the anonymous referees for their valuable time and careful comments, which improved the quality of this paper. F. Novoa-Mu˜noz wishes to thank his institution, University of B´ıo-B´ıo (Chile) and the scholarship given by the Chilean Ministry of Education through the Superior Education MECESUP Program 2, which make his doctorate studies possible. M.D. Jim´enez-Gamero acknowledges financial support from grant MTM2014-55966-P of the Spanish Ministry of Economy and Competitiveness. References Andrews, D.W.K. (1999) Estimation when a parameter is on a boundary. Econometrica, 67, 1341–1383. Andrews, D.W.K. (2000) Inconsistency of the bootstrap when a parameter is on the boundary of the parameter space. Econometrica, 68, 399–405. Burke, M.D. (2000). Multivariate tests-of-fit and uniform confidence bands using a weighted bootstrap. Statistics & Probability Letters, 46, 13–20. Crockett, N. G. (1979). A quick test of fit of a bivariate distribution. In Interactive Statistics, D. McNeil (ed.), 185–191. Amsterdam: North-Holland. Delhing, H. and Mikosch, T. (1994). Random quadratic forms and the bootstrap for U-statistics. Journal of Multivariate Analysis, 51, 392–413. Duchesne, P. and Lafaye De Micheaux, P. (2010). Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods. Computational Statistics and Data Analysis, 54, 858–862.