Chi-Squared Goodness-of-Fit Test

Section 10 Chi-squared goodness-of-fit test. Example. Let us start with a Matlab example. Let us generate a vector X of 100 i.i.d. uniform random variables on [0, 1] : X=rand(100,1). Parameters (100, 1) here mean that we generate a 100×1 matrix or uniform random variables. Let us test if the vector X comes from distribution U [0, 1] using �2 goodness-of-fit test: [H,P,STATS]=chi2gof(X,’cdf’,@(z)unifcdf(z,0,1),’edges’,0:0.2:1) The output is H = 0, P = 0.0953, STATS = chi2stat: 7.9000 df: 4 edges: [0 0.2 0.4 0.6 0.8 1] O: [17 16 24 29 14] E: [20 20 20 20 20] We accept null hypothesis H0 : P = U [0, 1] at the default level of significance � = 0.05 since the p-value 0.0953 is greater than �. The meaning of other parameters will become clear when we explain how this test works. Parameter ’cdf’ takes the handle @ to a fully specified c.d.f. For example, to test if the data comes from N (3, 5) we would use ’@(z)normcdf(z,3,5)’, or to test Poisson distribution �(4) we would use ’@(z)poisscdf(z,4).’ It is important to note that when we use chi-squared test to test, for example, the null hypothesis H0 : P = N (1, 2), the alternative hypothesis is H0 : P ∞= N (1, 2). This is different from the setting of t-tests where we would assume that the data comes from normal distribution and test H0 : µ = 1 vs. H0 : µ ∞= 1. PSfragPearson’s replacements theorem. Chi-squared goodness-of-fit test is based on a probabilistic result that we will prove in this section. �1 �2 B1 p1 B2 p2 �r Br pr ... Figure 10.1: Let us consider r boxes B1 , . . . , Br and throw n balls X1 , . . . , Xn into these boxes inde­ pendently of each other with probabilities P(Xi ≥ B1 ) = p1 , . . . , P(Xi ≥ Br ) = pr , so that p1 + . . . + pr = 1. Let �j be a number of balls in the jth box: �j = #{balls X1 , . . . , Xn in the box Bj } = n � l=1 I(Xl ≥ Bj ). On average, the number of balls in the jth box will be npj since E�j = n � l=1 EI(Xl ≥ Bj ) = n � l=1 P(Xl ≥ Bj ) = npj . We can expect that a random variable �j should be close to npj . For example, we can use a Central Limit Theorem to describe precisely how close �j is to npj . The next result tells us how we can describe the closeness of �j to npj simultaneously for all boxes j ≈ r. The main difficulty in this Thorem comes from the fact that random variables �j for j ≈ r are not independent because the total number of balls is fixed �1 + . . . + �r = n. If we know the counts in n − 1 boxes we automatically know the count in the last box. Theorem.(Pearson) We have that the random variable r � (�j − npj )2 d 2 � �r−1 np j j=1 converges in distribution to �2r−1 -distribution with (r − 1) degrees of freedom. Proof. Let us fix a box Bj . The random variables I(X1 ≥ Bj ), . . . , I(Xn ≥ Bj ) that indicate whether each observation Xi is in the box Bj or not are i.i.d. with Bernoulli distribution B(pj ) with probability of success EI(X1 ≥ Bj ) = P(X1 ≥ Bj ) = pj and variance Var(I(X1 ≥ Bj )) = pj (1 − pj ). Therefore, by Central Limit Theorem the random variable �n �j − npj l=1 I(Xl ≥ Bj ) − npj � � = npj (1 − pj ) npj (1 − pj ) �n l=1 I(Xl ≥ Bj ) − nE ≤ = �d N(0, 1) nVar converges in distribution to N(0, 1). Therefore, the random variable �j − npj d � � 1 − pj N(0, 1) = N (0, 1 − pj ) ≤ npj converges to normal distribution with variance 1 − pj . Let us be a little informal and simply say that �j − npj � Zj ≤ npj where random variable Zj � N (0, 1 − pj ). We know that each Zj has distribution � 2 N (0, 1 − pj ) but, unfortunately, this does not tell Zj will be, because as we mentioned above r.v.s �j us what the distribution of the sum are not independent and their correlation structure will play an important role. To compute the covariance between Zi and Zj let us first compute the covariance between �i − npi �j − npj and ≤ ≤ npi npj which is equal to �i − npj �j − npj 1 E ≤ = ≤ (E�i �j − E�i npj − E�j npi + n2 pi pj ) ≤ npi npj n pi pj 1 1 = ≤ (E�i �j − npi npj − npj npi + n2 pi pj ) = ≤ (E�i �j − n2 pi pj ). n pi pj n pi pj To compute E�i �j we will use the fact that one ball cannot be inside two different boxes simultaneously which means that I(Xl ≥ Bi )I(Xl ≥ Bj ) = 0. (10.0.1) Therefore, E�i �j = E n �� l=1 = E � � l=l� I(Xl ≥ Bi ) n �� � l� =1 � � I(Xl� ≥ Bj ) = E I(Xl ≥ Bi )I(Xl� ≥ Bj ) l,l� I(Xl ≥ Bi )I(Xl� ≥ Bj ) +E � l=l � � I(Xl ≥ Bi )I(Xl� ≥ Bj ) �� � this equals to 0 by (10.0.1) = n(n − 1)EI(Xl ≥ Bj )EI(Xl� ≥ Bj ) = n(n − 1)pi pj . Therefore, the covariance above is equal to � 1 � ≤ 2 n(n − 1)p p − n p p ≤ i j i j = − pi pj . n pi pj To summarize, we showed that the random variable r � (�j − npj )2 j=1 npj � r � Zj2 . j=1 where normal random variables Z1 , . . . , Zn satisfy ≤ EZi2 = 1 − pi and covariance EZi Zj = − pi pj . To prove the Theorem it remains to show that this covariance structure of the sequence of (Zi ) implies that their sum�of squares has �2r−1 -distribution. To show this we will find a Zi2 . different representation for Let g1 , . . . , gr be i.i.d. standard normal random variables. Consider two vectors ≤ ≤ g = (g1 , . . . , gr )T and p = ( p1 , . . . , pr )T ≤ ≤ and consider a vector g − (g · p)p, where g · p = g1 p1 + . . . + gr pr is a scalar product of g and p. We will first prove that g − (g · p)p has the same joint distribution as (Z1 , . . . , Zr ). To show this let us consider two coordinates of the vector g − (g · p)p : r � ≤ ≤ i : gi − gl p l pi th l=1 and j th r � ≤ ≤ : gj − gl pl pj l=1 and compute their covariance: r r � � ≤ ≤ � ≤ ≤ �� gl pl pj E gi − gl pl pi gj − � l=1 l=1 n � ≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤ = − pi pj − pj pi + pl pi pj = −2 pi pj + pi pj = − pi pj . l=1 (10.0.2) Similarly, it is easy to compute that This proves (10.0.2), which provides us with another way to formulate the convergence, namely, we have But this vector has a simple geometric interpretation. Since vector p is a unit vector: vector Vl = (p . g ) p is the projection of vector g on the line along p and, therefore, vector Vz = g (p . g ) p will be the projection of g onto the plane orthogonal to p, as shown in figure 10.2. - Figure 10.2: New coordinate system, Let us consider a new orthonormal coordinate system with the first basis vector (first axis) equal t o p . In this new coordinate system vector g will have coordinates obtained from g by orthogonal transformation V = (p, p2 , . . . , pr ) that maps canonical basis into this new basis. But we proved in Lecure 4 that in that case g1� , . . . , gr� will also be i.i.d. standard normal. From figure 10.2 it is obvious that vector V2 = g − (p · g)p in the new coordinate system has coordinates (0, g2� , . . . , gr� )T and, therefore, |V2 |2 = |g − (p · g)p|2 = (g2� )2 + . . . + (gr� )2 . But this last sum, by definition, has �2r−1 distribution since g2� , · · · , gr� are i.i.d. standard normal. This finishes the proof of Theorem. Chi-squared goodness-of-fit test for simple hypothesis. Suppose that we observe an i.i.d. sample X1 , . . . , Xn of random variables that take a finite number of values B1 , . . . , Br with unknown probabilities p1 , . . . , pr . Consider hypotheses H0 : pi = p�i for all i = 1, . . . , r, H1 : for some i, pi ∞= p�i . If the null hypothesis H0 is true then by Pearson’s theorem T = r � (�i − np� )2 i i=1 np�i �d �2r−1 where �i = #{Xj : Xj = Bi } are the observed counts in each category. On the other hand, if H1 holds then for some index i, pi ∞= p�i and the statistics T will behave differently. If pi is the true probability P(X1 = Bi ) then by CLT �i − npi d � N(0, 1 − pi ). ≤ npi If we rewrite �i − np�i �i − npi + n(pi − p�i ) ≤ � = ≤ � = npi npi � pi �i − npi ≤ pi − p�i + n ≤ � ≤ p�i npi pi then the first term converges to N(0, (1 − pi )pi /p�i ) and the second term diverges to plus or minus → because pi = ∞ p�i . Therefore, (�i − np�i )2 � +→ np�i which, obviously, implies that T � +→. Therefore, as sample size n increases the distri­ bution of T under null hypothesis H0 will approach �2r−1 -distribution and under alternative hypothesis H1 it will shift to +→, as shown in figure 10.3. 0.1 0.09 0.08 H0 : T � �2r−1 0.07 0.06 0.05 0.04 H1 : T � +→ 0.03 0.02 PSfrag replacements � 0.01 0 0 10 c 20 30 40 50 60 Figure 10.3: Behavior of T under H0 and H1 . Therefore, we define the decision rule � H1 : T ≈ c α= H2 : T > c. We choose the threshold c from the condition that the error of type 1 is equal to the level of significance � : � = P1 (α = ∞ H1 ) = P1 (T > c) � �2r−1 (c, →) since under the null hypothesis the distribution of T is approximated by �2r−1 distribution. Therefore, we take c such that � = �2r−1 (c, →). This test α is called the chi-squared goodness­ of-fit test. Example. (Montana outlook poll.) In a 1992 poll 189 Montana residents were asked (among other things) whether their personal financial status was worse, the same or better than a year ago. Worse Same Better Total 58 64 67 189 We want to test the hypothesis H0 that the underlying distribution is uniform, i.e. p1 = p2 = p3 = 1/3. Let us take level of significance � = 0.05. Then the threshold c in the chi-squared test α= � H0 : T ≈ c H1 : T > c is found from the condition that �23−1=2 (c, →) = 0.05 which gives c = 5.9. We compute chi-squared statistic T = (58 − 189/3)2 (64 − 189/3)2 (67 − 189/3)2 + + = 0.666 < 5.9 189/3 189/3 189/3 which means that we accept H0 at the level of significance 0.05. Goodness-of-fit for continuous distribution. Let X1 , . . . , Xn be an i.i.d. sample from unknown distribution P and consider the fol­ lowing hypotheses: � H0 : P = P 0 H1 : P = ∞ P0 for some particular, possibly continuous distribution P0 . To apply the chi-squared test above we will group the values of Xs into a finite number of subsets. To do this, we will split a set of all possible outcomes X into a finite number of intervals I1 , . . . , Ir as shown in figure 10.4. 0.4 0.35 p.d.f. of P0 0.3 0.25 0.2 PSfrag replacements p02 0.15 0.1 0.05 0 p0r p01 I1 I2 ··· Ir Figure 10.4: Discretizing continuous distribution. x The null hypothesis H0 , of course, implies that for all intervals P(X ≥ Ij ) = P0 (X ≥ Ij ) = p0j . Therefore, we can do chi-squared test for H0� : P(X ≥ Ij ) = p0j for all j ≈ r H1� : otherwise. Asking whether H0� holds is, of course, a weaker question that asking if H0 holds, because H0 implies H0� but not the other way around. There are many distributions different from P that have the same probabilities of the intervals I1 , . . . , Ir as P. On the other hand, if we group into more and more intervals, our discrete approximation of P will get closer and closer to P, so in some sense H0� will get ’closer’ to H0 . However, we can not split into too many intervals either, because the �2r−1 -distribution approximation for statistic T in Pearson’s theorem is asymptotic. The rule of thumb is to group the data in such a way that the expected count in each interval np0i = nP0 (X ≥ Ii ) ∼ 5 is at least 5. (Matlab, for example, will give a warning if this expected number will be less than five in any interval.) One approach could be to split into intervals of equal probabilities p0i = 1/r and choose their number r so that n np0i = ∼ 5. r Example. Let us go back to the example from Lecture 2. Let us generate 100 observa­ tions from Beta distribution B(5, 2). X=betarnd(5,2,100,1); Let us fit normal distribution N(µ, ν 2 ) to this data. The MLE µ̂ and ν̂ are mean(X) = 0.7421, std(X,1)=0.1392. Note that ’std(X)’ in Matlab will produce the square root of unbiased estimator (n/n − 1)ν̂ 2 . Let us test the hypothesis that the sample has this fitted normal distribution. [H,P,STATS]= chi2gof(X,’cdf’,@(z)normcdf(z,0.7421,0.1392)) outputs H = 1, P = 0.0041, STATS = chi2stat: 20.7589 df: 7 edges: [1x9 double] O: [14 4 11 14 14 16 21 6] E: [1x8 double] Our hypothesis was rejected with p-value of 0.0041. Matlab split the real line into 8 intervals of equal probabilities. Notice ’df: 7’ - the degrees of freedom r − 1 = 8 − 1 = 7.