Table Of ContentElectronic Journal of Statistics
Vol.0(2009)
ISSN:1935-7524
Testing distribution in deconvolution
problems
9 Denys Pommeret
0
0 Institut de Math´ematiques de Luminy - Case 907 - Universit´e de la M´editerran´ee - 13288
2 Marseille cedex 9, France e-mail: [email protected]
n Abstract: InthispaperweconsiderarandomvariableY contaminedby
a anindependentadditivenoiseZ.WeassumethatZhasknowndistribution.
J Ourpurposeistotestthe distributionof theunobserved randomvariable
Y.Weproposeadatadrivenstatisticbasedonadevelopmentofthedensity
7
of Y +Z, which is valid in the discrete case and in the continuous case.
2
Thetestisillustratedinbothcases.
] AMS 2000 subject classifications:Primary62G10; secondary62F05.
T
Keywordsandphrases:contaminateddata,Laguerrepolynomials,Meixner
S polynomials,Legendrepolynomials.
.
h
at Contents
m
[ 1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1
3 Some particular cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
v
3.1 Continuous positive case . . . . . . . . . . . . . . . . . . . . . . . 5
6
8 3.2 Continuous bounded case . . . . . . . . . . . . . . . . . . . . . . 5
1 3.3 Discrete case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4 3.4 Dependent case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
.
1 4 Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
0 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
9
5 Appendix: Orthogonalpolynomials . . . . . . . . . . . . . . . . . . . . 8
0
:
v
i 1. Introduction
X
r
a Consider the convolution
X =Y +Z, (1.1)
where Y and Z are independent r.v. with densities f and g w.r.t. a reference
measure µ. According to the nature of the variable, µ can be the Lebesgue
measure on R (or on interval), or the counting measure δ , where δ (x)=1
n n
if x=n and0 otherwise.The errortermZ is assumedtoPhavea knowndensity
g. Observing X instead of Y and Z yields to a deconvolution problem which
consistsindistinguishingthe twocomponentsofthe variable.Ourpurposeisto
test the distribution of Y.
1
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 2
Different authors have considered the nonparametric deconvolution problem
of estimating f or its associated distribution function from i.i.d. observations;
that is, to recover the distribution of Y using the contaminated measurements
X.Nonparametrickernel-typeestimationoff basedoncontaminateddatawas
studiedby(2),(6)and(3)amongothers.Thisproblemisalsorelatedtomixture
problem since (1.1) may be seen as a location mixture f(x m)Π(dm), with
−
mthelocationparameterandΠtheknownmixingdistrRibution(seeforinstance
(4) and (11)).
Recently, (7) proposed a new procedure for testing the density of Y. Their
methodisbasedonFouriertransformsofvariablesY,Z andX (viaitsnonpara-
metric estimator). In this paper we also use a nonparametric representation of
the density h of X to test the distribution of Y. But our approach is based on
a polynomial expansion of the density h under the null hypothesis. The com-
parison of this density with the empirical one allows to build a test statistic.
A data driven approach permits to select automatically the number of com-
ponents of the statistic. In fact, we consider the same problem than (7), but
instead of considering the problem as a L2 distance between a nonparametric
deconvolutionestimatorandasmoothedversionofthedensityf,werestrictour
attention to a known reference measure µ and its associatedL2 basis where we
represent all densities. In (7), the authors need to make different assumptions
onFouriertransformsandthe use ofkerneldensity estimatorsrequiresclassical
supplementaryassumptionsonthebandwidth.Intheprocedureweproposedin
this paper we did not need such parameter thanks to the use of a data driven
technic.Thenouradditionalassumptionisthatdensitiesaresquaredintegrable
with respect to µ. Also, for asymptotic results we make technical conditions on
eigenvaluesofthevariancematrix.Butinpractice,thetestmaybeusedwithout
data driven procedure for a large enough number of statistics components.
More recently, such data driven approach has been used for deconvolution
problem: in (9), the problem of testing the density of Y is considered.The idea
is based on the convolution formula obtained in the continuous case; that is,
when the reference measure µ is the Lebesgue one. In that case the author
constructedascoreteststatistic whichis combinedwith amodelselectionrule.
In our paper, the data driven method is close to the work of (9). Also our test
statisticcanberelatedtothescoreone(seeRemark2.1).Howeverwedidnotuse
theclassicalconvolutionformula,whatallowsustoconsiderdiscreteconvolution
problems, as that illustrated in our simulations. Also it permits to introduce
a dependance between Y and Z as we underline it then. But the extension
in a multivariate setting, as done in (7), requires here boring calculus because
orthogonalpolynomialsarenotwelldevelopedinthisframe.Finally,thenovelty
of our test is that we can also test discrete convolution while in our knowledge
it is essentially the continuous case which is studied in the literature. Also, we
extend our study to the dependent case, testing the conditional distribution of
Y Z.
|
Thepaperisorganisedasfollows:inSection2weintroducethemethodbased
onthe polynomialexpansionofthe (null) densityandweproposeadata driven
procedure for testing the distribution of the contaminated density. In Section 3
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 3
a simulation study shows some comportments of the proposed test.
2. Method
Let X , ,X be i.i.d. random variables with density h, satisfying the convo-
1 n
···
lution formula (1.1). We assume that g is known and we want to test
H :f =f against f =f ,
0 0 0
6
where f is a fixed density. Let µ be a probability measure on R with density
0
m andsuchthatbothdistributions of Y andX aredominatedby µ.Denote by
= Q ;i = 0,1, an associated basis of dense orthogonal functions with
i
B { ···}
respect to µ. When h belongs to L2(µ) (what we shall suppose) we have the
following expansion:
h(x) = a Q (x), (2.1)
i i
nXN
∈
where a =E(Q (X)m(X)). Our method consists in comparing the coefficients
i i
a’sbasedonthisexpansionwiththoseobtainedunderthenullhypothesis.Under
H , equality (2.1) may be rewritten as
0
h(x) = E(Q (Y +Z)m(Y +Z))Q (x)
i i
Xi N
∈
= α Q (x).
i i
Xi N
∈
Testing H is equivalent to test
0
H :a =α i=1,2, ,
0 i i
∀ ···
Under H these coefficients α’s can be easy to calculate with a good choice of
0
the reference measure and its associated polynomials (see illustrations below).
Thus,anaturalstatisticcanbe constructedasfollows:forsomeintegerk,write
B = (b , ,b ),
k 1 k
···
b b
where
n
1
b = Q (X )m(X ) α .
j j i i j
√n −
(cid:0)Xi=1 (cid:1)
b
By the Central Limit Theorem, we have the following convergence under H :
0
Uk =Σk−1/2BkT −→L N(0,I),
where Σ is the k k covariance matrix Var(B ).
k k
×
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 4
Under H , T = U 2 is asymptotically Chi-squared distributed with k
0 k k
k k
degreesoffreedomandthenwe rejectthe nullhypothesis forlargevalues ofT .
k
But the arbitrary choice of k is the weakness of this method and we adapt a
data driven method to select the number of components in the test statistics.
Imitating the work of (10), we consider an increasing sequence of number of
components k(n) such that lim k(n) = . Our selection rule is based on
n
→∞ ∞
the Schwarz criteria. Write
S =min argmax(T klog(n)) . (2.2)
n k
−
(cid:8)1 k k(n) (cid:9)
≤ ≤
Then our data driven test statistic is T .
Sn
Proposition 2.1. Let λ be the smallest eigenvalue of Σ . Assume that
k(n) k(n)
loλgkk(n(n)) = oP(logn). Then, under H0, TSn converges in distribution to a Chi-
squared random variable with 1 degree of freedom.
Proof. Under H , T converges to a Chi-squared random variable with one de-
0 1
gree of freedom and then we have to show that P(S =1) tends to zero.
n
Since (S = k) implies (T klog(n) > T log(n)) we have P(S = k)
n k 1 n
− − ≤
P(T >(k 1)log(n)). Then we have
k
−
k(n) k(n)
P(S >2)= P(S =k) P T >(k 1)log(n) . (2.3)
n n k
≤ (cid:18) − (cid:19)
kX=2 kX=2
Using the fact that λ1k =supX∈R∗k X′XΣ′−kX1X, we obtain
B 2
Tk =tBkΣ−k1Bk ≤ k λkk (2.4)
k
and
P T >(k 1)log(n) P B 2 >λ (k 1)log(n) . (2.5)
k k k
(cid:18) − (cid:19)≤ (cid:18)k k − (cid:19)
As Σ is an Hilbert Schmidt operator, its trace is bounded by a constant, say
k
M, independently of k, and we have (under H )
0
E( B 2) = Var(B )
k k
k k
k
= E(Q (X)m(X))2 (α )2
i i
−
Xi=1(cid:0) (cid:1)
= Trace(Σ ) M.
k
≤
Combining (2.5) with Markov inequality and the above result yields
M
P T >(k 1)log(n) . (2.6)
k
(cid:18) − (cid:19)≤ (k 1)log(n)λ
k
−
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 5
Finally, using a result on harmonic sum we get n 1 = O (logn) and
k=1 k +∞
we obtain P
M logk(n)
P(S >2) .
n
≤ inf λ (cid:18) logn (cid:19)
1 k k(n) k
≤ ≤
Since matrices Σ are embedded, λ is a decreasing sequence and we get the
k k
result.
Remark2.1. Itis clear thatourhypotheses canberewritten H :θ =0against
0
H : θ = 0, where θ = E(V(k)) with V(k) = (Q (X)m(X) α ) . Then
1 6 i − i i=1,···,k
T coincides with the score statistic if the maximum likelihood estimator θ of
k
θ equals the empirical mean of the sample of the V(k)’s, that is θ = Bk, as it
√n b
is the case for instance when the distribution of V(k) belongs to an exponential
b
family.
3. Some particular cases
3.1. Continuous positive case
Assume that Y and Z are continuous positive random variables. We choose µ
the exponentialdistribution with mean1 and Q =L its associatedLaguerre
i i,1
orthogonalpolynomials(see(1)).InAppendixwerecallsomebasicpropertiesof
LaguerreandgeneralizedLaguerrepolynomials.Equality(2.1)mayberewritten
as
h(x) = E(Q (Y +Z)exp( Y)exp( Z))Q (x)
i i
− −
Xi N
∈
= C E(E (Y)E (Z)) Q (x),
i,s s,u i s,v i
Xi N(cid:0)Xs i − (cid:1)
∈ ≤
where E (x) = L (x)exp( x)), C are coefficients given in Appendix and
s,u s,u i,s
−
u,varearbitrarpositiverealssatisfyingu+v=1.UnderH ,usingtheprevious
0
decomposition, we have:
Σ = E(Q (Y +Z)Q (Y +Z)exp( 2Y)exp( 2Z)) α α
ij i j i j
− − −
= C C E(E (Y)E (Z)) α α ,
s,i t,j (s,t),u (i s,j t),v i j
Xs iXt j − − −
≤ ≤
where E (x)=L (x)L (x)exp( 2x)).
(s,t),u s,u t,u
−
3.2. Continuous bounded case
Assume that Y and Z arebounded, w.l.o.g.taking values in [0;1].We chooseµ
theuniformdistributionandQ itsassociatedLegendreorthogonalpolynomials
i
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 6
described in Appendix. Then equality (2.1) may be rewritten as
h(x) = E(Q (Y +Z))Q (x)
i i
Xi N
∈
= C E(YsZt) Q (x),
i,s,t i
Xi N(cid:0)sX+t i (cid:1)
∈ ≤
where C are coefficients obtained by expansion of Q . Under H :
i,s,t i 0
Σ = E(Q (Y +Z)Q (Y +Z)) α α
ij i j i j
−
= C C E(Ys+u)E(Zt+v) α α .
i,s,t j,u,v i j
−
s+Xt iu+Xv j
≤ ≤
3.3. Discrete case
Assume that Y and Z are discrete random variables taking values in N. We
choose µ the geometric distribution with probability µ(x) = px(1 p) for x =
−
0,1, ,p (0,1)andQ =M itsassociatedMeixnerorthogonalpolynomials
i i,1
··· ∈
described in Appendix. Equality (2.1) may be rewritten as
h(x) = E(Q (Y +Z)pYpZ(1 p))Q (x)
i i
−
Xi N
∈
= C E(E (Y)E (Z)) Q (x),
i,s s,u i s,v i
Xi N(cid:0)Xs i − (cid:1)
∈ ≤
whereE (x)=E(M (x)px(1 p)1/2),C arecoefficientsgiveninAppendix
s,u s,u i,s
−
and u,v are arbitrar positive reals satisfying u+v = 1. Under H , using the
0
previous decomposition, we have :
Σ = E(Q (Y +Z)Q (Y +Z)p2Yp2Z(1 p)2) α α
ij i j i j
− −
= C C E(E (Y)E (Z)) α α ,
s,i t,j (s,t),u (i s,j t),v i j
Xs iXt j − − −
≤ ≤
where E (x)=M (x)M (x)p2x(1 p).
(s,t),u s,u t,u
−
3.4. Dependent case
WhenY andZ arenotindependent,problem(1.1)canbetraitedbycondition-
ningwithrespecttothenoiseZ.Inthatcase,hypothesesconcerntheconditional
distributionofY Z;thatisH :f =f (.,Z).Wecanusethesameapproach,
0 Y Z 0
| |
replacing coefficients α by the quantities
i
E(E(Q (Y +Z)m(Y +Z)Z)).
i
|
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 7
4. Illustrations
In this section we present the results of two simulation studies of our testing
procedure. We consider i.i.d. data X , ,X generated from two convolution
1 n
···
models (1.1) satisfying:
First model (Mod1): Y has exponential distribution with mean 1 and Z is
Chi-squareddistributedwith1degreeoffreedom.Threealternativesarestudied:
Alt1: instead of (1.1), X is a mixture with two components:
50% exponential with mean 2 and 50% Chi squared with 2 degrees of freedom.
Alt2: convolution (1.1) with both Y and Z e−xponential distributed with mean 1.
Alt3: convolution (1.1) with both Y and Z Chi Squared distributed with degree 1.
−
Second model (Mod2): Y has Poisson distribution with mean 1 and Z has
Geometric distributed with mean 1. Three alternatives are proposed:
Alt4: instead of (1.1), X is a mixture with two components:
50% Poisson with mean 2 and 50% Geometric with mean 2.
Alt5: convolution (1.1) with both Y and Z Poisson distributed with mean 1.
Alt6: convolution (1.1) with both Y and Z Geometric distributed with mean 1.
It is clear that for Models 1 and 2, the two convolution’s components have dis-
tributionswithrelativelyclosecharacteristicsandweareinterestedindetecting
a confusion between these components (Alternatives 2-3 and Alternatives 5-6)
or we are interested in detecting a mixture of these two components instead of
a convolution (Alternative 1 and Alternative 4).
For each model and alternative, we compute the test statistic based on a
sample size n = 50,100 and 500 for a theoretical level α = 5%. The empirical
level (resp. power) of the test is defined as the percentage of rejection of the
null hypothesis over 10000 replications of the test statistic under H (resp.
0
under Alternative). We can see that for Alternative 2, the power is weak for
small samples. For Alternative 4 the power is very low. Then the mixture from
Alternative 4 and the convolution from Model 2 are very close and it is very
difficult to distinguish them.
Figure 1 here
Figure 2 here
References
[1] Abramowitz, M. and Stegun, I. A. (Eds.) (1972). Orthogonal Polyno-
mials. Ch. 22 in Handbook of Mathematical Functions, 9th printing. New
York: Dover.
[2] Carroll, R. J. and Hall, P. (1988). Optimal rates of convergence for
deconvolving a density. J. Amer. Statist. Assoc. 83 1184–1186.
[3] Delaigle, A. and Gijbels, I. (2004). Bootstrap bandwidth selection in
kernel density estimation from a contaminated sample. Ann. Inst. Statist.
Math. 56 19–47.
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 8
[4] Devroye, L. (1989). Consistent deconvolution in density estimation.
Canad. J. Statist. 17 235–239.
[5] Eagleson, G. 1964. Polynomial expansion of bivariate distribution. Ann.
Math. Statist. 35 1208–1215.
[6] Fan, J. (1991). On the optimal rates of convergence for nonparametric de-
convolution problems. Ann. Statist. 19 1257–1272
[7] Holzmann, H.,Bissantz, N.,andMunk, A. (2007).Density testing ina
contaminated sample. J. Multiv. Anal. 98 57–75
[8] Koudou, A.E. (1996). Probabilit´es de Lancaster, Exp. Math. 14 247–275.
[9] Langovoy,M.A.(2008).Data-drivenefficientscoretestsfordeconvolution
hypotheses. Inverse Problem 24 025028.
[10] Ledwina, T. (1994). Data-driven version of neymans smooth test of Fit.
J. Amer. Statist. Assoc.. 89 1000–1005.
[11] Zhang, C.H.(1990).Fouriermethods forestimating mixing densities and
distributions. Ann. Statist. 18 806–831.
5. Appendix: Orthogonal polynomials
We follow the notation used in (1).
Laguerre polynomials Laguerre polynomials L ;n N are defined by
n,1
{ ∈ }
their recurrence relation
L =1 L (x)=1 x
0,1 1,1
−
(i+1)L (x) = (2i+1 x)L (x) iL (x)
i+1,1 i,1 i 1,1
− − −
They are orthogonal w.r.t. the exponential distribution with density f(x) =
exp( x)onR+.Forα>0,GeneralizedLaguerrepolynomials L ;n N are
n,α
− { ∈ }
orthogonalw.r.t.theGammadistributionwithdensityf(x,α)=exp( x)xα−1Γ(α)−1.
−
They satisfy
L =1 L (x)=α x
0,α 1,α
−
(i+1)L (x) = (2i+α x)L (x) (i+α 1)L (x)
i+1,α i,α i 1,α
− − − −
Their norms are given by L 2 = α 2i(i!)Γ(i+α)/Γ(α). In the exponential
i,α −
k k
case we simply get L 2 = (i!)2. These polynomials satisfy the following
i,1
k k
relation that we used in our simulations study (see (5) or (8) for an idea of the
proof) for α=u+v, u>0, v >0:
L (y+z) = C usL (y)vn sL (z) (5.1)
n,α s,n s,u − n s,v
sXn −
≤
where C =n!/(s!(n s)!).
s,n
−
Legendre polynomials For µ the uniform distribution on (0,1), associated
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 9
orthogonal polynomials are (shifted) Legendre polynomials defined by the re-
currence relation :
P =1 P (x)=2x 1
0 1
−
(n+1)P (x) = (2n+1)xP (x) nP (x),
n+1 n n 1
− −
and sastifying P 2 = 1P (x)2dx=(2n+1) 1.
k nk 0 n −
R
Meixner polynomials For µ the Pascal distribution, µ(x) = (1 p)px, p
(0,1), x N, the associatedorthogonalpolynomials satisfy the follo−wing recur∈-
∈
rence formula:
M =1 M (x)=1 p x/p
0 1
− −
p/(1 p)M (x) = ((p 1)x+(1+p)n+p)M (x) (1 p)n2M (x),
n+1 n n 1
− − − − −
with M 2 = p n(1 p)2n(n!)2. They are particular case of Meixner polyno-
n −
k k −
mials M with b = 1 and c = p. These polynomials also satisfy a relation
n;b,c
similar to (5.1) that we used for Illustrations:
M (y+z) = C M (y)M (z) (5.2)
n;b,c s,n s,u,c n s,v,c
sXn −
≤
where u+v =b and C =n!/(s!(n s)!).
s,n
−
100
80
Mod1HlevelL
values 60 Alt1HpowerL
Empirical 40 Alt2HpowerL
Alt3HpowerL
20
100 200 300 400 500
Samplesize
Fig 1. Empirical level for Model 1 and empirical powers for Alternatives 1-3, with α=5%,
n=50,100,500
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
/ 10
80
Mod2HlevelL
values60 Alt4HpowerL
Empirical40 Alt5HpowerL
20 Alt6HpowerL
100 200 300 400 500
Samplesize
Fig 2. Empirical level for Model 2 and empirical powers for Alternatives 4-6, with α=5%,
n=50,100,500
imsart-ejs ver. 2008/08/29 file: ejs_2009_364.tex date: January 27, 2009
