Goodness Of Fit Tests For Normality

The Pearson and Kolmogorov-Smirnov tests are used to test whether a random data set comes from a normally distributed population.

by David K. Neal

We often hear of measurements that are commonly assumed or found 
to be normally distributed.  Such measurements include IQ scores, heights 
or lengths of animals, and standardized test scores, among others.  However 
given a data set, it requires quite a bit of effort to verify that the data comes 
from a normal population.  In this article, we shall describe two goodness of fit tests for normality and provide Mathematica codes which can quickly implement the procedures.  We shall also use the procedures to test for uniformity.

Such tests are important in statistics since many theorems and concepts 
are based on the underlying assumption that one's sample data comes from 
a normally distributed population.  For example, suppose we wish to construct 
a (1-a)x100% confidence interval for the unknown mean µ of a population 
Ω.  We obtain a random sample of data x1, x2, . . . , xn, let x-  be the sample mean, and let S2 be the unbiased sample variance.   If we know that Ω is normally distributed, then for all sample sizes n ≥ 2, the expression (x-  - µ)/(S/n1/2) 
follows a T distribution with n-1 degrees of freedom  [Hogg and Tanis 1993, 
page 303].  We denote this by (x-  - µ)/(S/n1/2) ~ T(n-1).  If we let t be the 
T-score such that P(T(n-1) ≥ t) =  a/2, then the confidence interval is given by  
(x-  - t S/n1/2, x-  + t S/n1/2).  A normal population allows us to use small sample sizes and an exact distribution.  If  Ω were not normal, then we would need much larger sample sizes; moreover, the distribution of (x-  - µ)/(S/n1/2) could only be approximated by a standard normal distribution.

If Ω is normal, then we can also find a confidence interval for the variance
s2 since (n-1)S2/s2 follows a Chi Square distribution with n-1 degrees of 
freedom [Hogg and Tanis 1993, page 274].  If we have several independent 
normal populations with a common variance, then we can perform one-way 
or two-way analysis of variance.  Also, correlation analysis between two 
independent populations requires that one of the populations be normally distributed [Hogg and Tanis 1993, Section 9.4].  Therefore, it is nice to be able 
to assume that one has a normal population; but we should not overlook an 
actual test for whether or not such an assumption is warranted.

We shall first describe the Pearson test which is more common though maybe not as powerful as the Kolmogorov-Smirnov test.  However, it has its advantages for educational purposes since it is easily taught using basic distributions.  Moreover, it is easily adaptable to both simple and composite hypothesis tests and can be used for testing both discrete and continuous distributions  The author has used this method in elementary descriptive statistics courses as well as in undergraduate calculus-based mathematical statistics.  Students routinely use this program to test for normality before using their data sets for further statistical analysis.

THE HYPOTHESIS TESTS

Given a population Ω, we would like to perform a test in order to accept or 
reject the claim that Ω is normally distributed.  We shall consider two cases
for the null hypothesis Ho:  (i) Ω is normal with mean µ; and (ii) Ω is normal.  
The alternative hypotheses are respectively  (i) Either Ω is not normal or the 
mean does not equal µ; and (ii) Ω is not Normal.  In each case, we let a be the 
level of significance which is the probability of rejecting Ho given that Ho is true.

We begin the test with a random sample of data points x1, x2, . . . , xn  from Ω and again let x-  be the sample mean and S2 be the unbiased sample variance.  In each case, we assume the null hypothesis to be true and then compare what should happen with the data to what actually does happen.  If the comparison is too far off, then we must reject the null hypothesis.  More specifically, if that which does happen occurs with a probability less than a, then we say that it is too unlikely to occur if Ho were true.  We therefore reject Ho.

PEARSON GOODNESS OF FIT TEST

We consider the first null hypothesis and assume that Ω is normal with mean 
µ and unknown variance s2.  If X is a random measurement from Ω, then 
Z = (X - µ)/s follows a standard normal distribution.  Since W = (n-1)S2/s2 ~ C2(n-1), we obtain a T distribution with n-1 degrees of freedom by 
Z/Sqrt[W/(n-1)] = (X - µ)/S  [Hogg and Tanis 1993, page 299].  Hence, if  yi = 
(xi - µ)/S, for 1 ≤ i ≤ n, then the data {yi} is a random sample of measurements
from the distribution Y ~ T(n-1).  We note that the sample deviation S is the required variable in this distribution.  Thus, we do do not need to know the true value of the variance s2.

We next divide the real line (-∞, ∞) into m equally probable subdivisions or “bins”.  As a rule of thumb, we will require that each bin have an expected frequency of around 5 data points.  Thus, the number of bins will be the greatest integer less than or equal to n/5, where n is the sample size.  The chance of a data point lying in any one bin will then be 1/m.  To divide the real line in such a manner, we simply choose uk so that P(T(n-1) ≤ uk) = k/m for k = 1, . . . , m-1.  We then have the required partition {-∞, u1, u2, . . . , um-1, ∞ }.  For larger sample sizes n, we would necessarily have more bins.  We can control the bin size though by increasing our choice of expected frequency to a more suitable amount.

We then count how many of the measurements yi lie in each bin and denote this frequency by fk, for 1 ≤ k ≤ m.   By construction, the expected number of measurements that should fall in each bin is ek = n/m.  We wish to compare the actual frequencies fk with these expected frequencies ek.  If the differences are too large, then we must reject the null hypothesis.

The test statistic that we use for the test is defined by  
                                                          m
                                                   Q = ∑ (fk - ek)2/ek.
                                                         k=1
For large samples, Q approximately follows a Chi Square distribution with 
m-1 degrees of freedom, C2(m-1) [Bain and Engelhardt 1992 page 454].

To complete the test, we simply compute the right tail value P(C2(m-1) ≥ Q) created by the test statistic.  It is preferable to compare the tail value with a, rather than comparing the test statistic with a chi-square score, since the tail value gives the probability of the data occurring as it did given that Ho were true.  If the probability is too small, that is if P(C2(m-1) ≥ Q) < a, then we must reject Ho.

THE COMPOSITE TEST

We use the second null hypothesis to test for normality without a specified 
mean.  If we reject Ho in the first case, then we still do not know if Ω is normal. 
It might be that Ω is not normal or that the mean is just not what we thought.  
If we subsequently accept Ho in the second case, then we are saying that Ω is normal but that the mean is just not known.  We could then find a confidence interval for the mean.

For the second null hypothesis, we only initially assume that Ω is normal 
and then estimate µ with the maximum likelihood estimator which is x-.  We 
now re-define yi by yi = (xi -x-)/S, for 1 ≤ i ≤ n.  We define the frequencies fk, the expected values ek, and the test statistic Q in the same manner as above.  However, since we have used one maximum likelihood estimate, we lose this one degree of freedom and Q is now approximately C2(m-2) [Bain and Engelhardt 1992 page 454].  We note that since S is the required variable in this distribution, we do not lose another degree of freedom for using S.  We now reject the null hypothesis if P(C2(m-2) ≥ Q) < a.

THE PROGRAM

Suppose we wish to test the claim, at the .05 level of significance, that combined SAT scores are normally distributed with a mean score of 1000.  We must base our decision solely on a random sample of data.  Consider then the following sample of scores which is entered in Mathematica code:

data={1270,1220,1160,950,1070,1110,1220,1150,1440,850,
     1280,1020,880,1080,1200,1120,970,1210,1030,1020,
     1400,720,1170,1070,860,920,1350,910,980,800,1070,
     1320,900,1200,1260,830,1330,1070,810,1240,820,
     1140,1170,870,1010,760,810,830,1240,910,1230,1100,
     970,1090,1300,1070,1040,1320,1230,890,960,920,
     1260,1060,1040,990,1040,1300,950,1090,1020,1270,
     1160,1050,1130,880,750,1200,1010,1060,1010,1080,
     1090,940,1170,920,990,1100,1270,930,790,1080,740,
     960,1170,1050,1250,1090,1010,1200,1080,1200,1070,
     970,1140,1360,1260,1170,1160,1310,1090,1040,1070,
     890,1280,1270,1190,1070,1410,1250,1070,910,1230};

This data was obtained from the Chance Database on the World Wide Web 
(http://www. geom.umn.edu/docs/snell/chance/teaching_aids/ETS_validity.
html).  The scores, as part of a larger data set, were studied by the Educational
Testing Service to determine the validity of SAT scores.

To execute the program, we first need to enter the following special packages :

<<Statistics`ContinuousDistributions`

<<Graphics`Graphics`

Then we compute the sample size, sample mean, unbiased sample deviation, the number of bins based on a desired expected frequency of around 5, and the actual
expected frequency per bin.

SampleSize=n=Length[data]
Xbar=N[Mean[data]]
S=N[StandardDeviation[data]]
NumberOfBins=m=Floor[n/5]
ExpectedFrequencyPerBin=N[n/m]

Next we create the partition of equally probable bins.  We note that the partition
points are not equally spaced but they are symmetric about 0.

Y=StudentTDistribution[n-1];u[0]=-Infinity;
u[m]=Infinity;Do[u[i]=Quantile[Y,i/m],{i,1,m-1}];
Table[u[i],{i,0,m}]

Now we must enter the value of the mean to be tested during the test for normality.  Here, we are testing whether the distribution is normal with a mean 
of 1000; so we enter the following command :

TestMean=mu=1000;

(If we were just testing for normality, then we would enter the command
TestMean=mu=Xbar;)

Next, we define the adjusted data yi = (xi - µ)/S.  We then use a short algorithm to count how many of the yi lie in each bin and display these frequencies :

x=Sort[data];Clear[y];
Do[y[i]=N[(x[[i]]-mu)/S],{i,1,n}];

i=1;
k=0;While[y[k+1]<=u[1],k++];
f[1]=k;

z=f[1];
Do[k=0;While[u[i-1]<y[z+k+1]<=u[i],k++];
   f[i]=k;z=z+k,{i,2,m-1}];
f[m]=n-z;

frequencies=Table[f[k],{k,1,m}]

If desired, we can visually compare the fit with a bar chart.  We note that, by construction, the expected number of occurrences in each bin is n/m, which in this case is 5.125.  From the above frequency list, we can see that the fit begins to fail towards the rightmost bins.  We can observe this further with the following display.

BarChart[frequencies]

Next we enter our desired level of significance and compute the test statistic.

LevelOfSignificance=a=.05;

TestStat=Q=N[Sum[(f[k]-n/m)^2/(n/m),{k,1,m}]]

To complete the test, we compute the right tail value and compare it with the level of significance.  Recall that the test statistic Q follows a Chi Square distribution with either m-1 or m-2 degrees of freedom.  Here since we have specified the mean to be tested, we use m-1 degrees.  If we were only testing for normality and estimating the mean with Xbar, then we would use m-2.

X=ChiSquareDistribution[m-1];
RightTailValue=P=1-CDF[X,Q]

Result=If[P<a,RejectHypothesis,AcceptHypothesis]

So we reject the claim that SAT scores are normally distributed with a mean 
of 1000; for if this hypothesis were true, then the probability of the data occurring as it did would only be .00000483589.  The tail value is also the highest level of significance for which we could accept the claim.  Since this value is so low, we would reject the claim for all reasonable levels of significance.

The SAT scores might still come from a normally distributed population though.  The population will just not have mean 1000.  We can now perform 
the composite test for normality only.  First we return to the code TestMean=mu=1000; and change it to the code below :

TestMean=mu=Xbar;

We re-enter the codes starting from this point; but we must change the degrees of freedom at the end to m-2.  In doing so, our last outputs become :

TestStat=Q=N[Sum[(f[k]-n/m)^2/(n/m),{k,1,m}]]

X=ChiSquareDistribution[m-2];
RightTailValue=P=1-CDF[X,Q]

Result=If[P<a,RejectHypothesis,AcceptHypothesis]

So we no longer have evidence to reject the claim of normality.  If the entire population of combined SAT scores were normally distributed, then there would be a 31.1275% chance of the data being distributed as it was.  Since the right tail value is now high, we would accept the claim of normality at all reasonable levels of significance.

Since we can now assume that the population of SAT scores is normally distributed, we can also use the T distribution Y~T(n-1) to find a confidence 
interval for the mean.  Suppose we want a 98% confidence interval.  We first 
enter this level :

ConfLevel=r=.98;

Then we can find the T-score and find the confidence interval directly with the following codes:

tscore=t=Quantile[Y,(r+1)/2];
ConfInterval=
   {N[Xbar-t*S/Sqrt[n]],N[Xbar+t*S/Sqrt[n]]}

Or we could use the built-in command (which assumes a normally distributed population).

<<Statistics`ConfidenceIntervals`

N[MeanCI[data,ConfidenceLevel->r]]

So, we ultimately conclude that this sample of SAT scores did come from a normally distributed population of scores with a mean somewhere between
1044.7 and 1113.51.

TEST FOR UNIFORMITY

With a few alterations to the above codes, we can also use a collection of data points to test for uniformity.  We now test the null hypothesis that the data 
comes from a population which is uniformly distributed over an interval [c, d], where either none, one, or both of the parameters c and d are specified.  If we do not specify one of these endpoints, then we must use the maximum likelihood estimate of the endpoint.  The endpoint c is estimated by the minimum of the 
data while d is estimated by the maximum.

As an example, we shall test the hypothesis that the Mathematica random number generator creates a uniform distribution.  To create a sample of say 200 data points between 0 and 30, we use the following commands :

Clear[x];
Do[x[i]=Random[Real,{0,30}],{i,1,200}];
data=Table[x[i],{i,1,200}];

Before we test for uniformity, we can quickly perform the composite test for 
normality.  We shall create a partition with an expected frequency of 8 per bin.  If we go back to the beginning of the above program, use n/8 for NumberOfBins and use Xbar for TestMean, then we obtain a distribution that does not appear normal :

frequencies=Table[f[k],{k,1,m}]

BarChart[frequencies]

The following is the result of the test for normality with the above generated random data :

TestStat=Q=N[Sum[(f[k]-n/m)^2/(n/m),{k,1,m}]]

X=ChiSquareDistribution[m-2];
RightTailValue=P=1-CDF[X,Q]

Since the tail value is so low, we of course reject normality.

After our data has been generated or entered as above, we can begin the test 
for uniformity.  For our purposes here, it is not necessary to see the actual data points; but if one wishes to do so, merely enter the command data.  It is helpful though to determine the minimum and maximum values :

SampleSize=n=Length[data]
Min[data]
Max[data]

Now for our first case, we shall test the null hypothesis that the code
Random[Real,{0,30}] is uniformly distributed on [0, 30].  We first enter 
these parameters.  We note that we can specify any values for c and d as long as
c ≤ Min[data] and d ≥ Max[data].  We then choose an appropriate number 
of bins.  Here, we shall divide the interval [0, 30] into equally probable bins that each have an expected frequency of 8.  Since we are testing for uniformity, the resulting partition will consist of equidistant points as well as equally probable bins.

c=0;
d=30;

NumberOfBins=m=Floor[n/8]

u[0]=c;Do[u[i]=c+(d-c)i/m,{i,1,m}];
Table[N[u[i]],{i,0,m}]

We next use a short algorithm to determine the frequencies of the data lying in each bin :

y=Sort[data];z=0;
Do[k=0;While[y[[z+k+1]]<=u[i],k++];
f[i]=k;z=z+k,{i,1,m-1}];f[m]=n-z;
frequencies=Table[f[k],{k,1,m}]

BarChart[frequencies]

If desired, we can define a level of significance before computing the test statistic; however, our decision could be reached solely from the tail value.

LevelOfSignificance=a=.05;

TestStat=Q=N[Sum[(f[k]-n/m)^2/(n/m),{k,1,m}]]

Since we have specified both parameters, we have not had to estimate either.  
Therefore, the test statistic follows a Chi Square distribution with m-1 degrees
of freedom.  Using this distribution, we find our tail value :

X=ChiSquareDistribution[m-1];
RightTailValue=P=1-CDF[X,Q]

Since the tail value is relatively high, we would accept the hypothesis at all reasonable levels of significance.  We can though compare the tail value with 
the pre-assigned level :

Result=If[P<a,RejectHypothesis,AcceptHypothesis]

Hence, based on the generated data, we can accept the claim that the 
Random[Real,{0,30}] code is uniformly distributed over [0, 30].

THE COMPOSITE UNIFORM TEST

If we have a raw data set and we simply wish to test for uniformity of distribution without specifying parameters c and d, then we must first estimate these endpoints by maximum likelihood.  Instead of specifying the values of c and d above as 0 and 30, we change the entries as follows :

c=Min[data];
d=Max[data];

However, for each maximum likelihood estimate, the test statistic loses a degree of freedom.  In this case, we would change X at the end as follows :

X=ChiSquareDistribution[m-3];
righttailvalue=P=1-CDF[X,Q]

If we had specified one parameter and estimated the other, then X would have m-2 degrees of freedom.

THE KOLMOGOROV-SMIRNOV TEST
     
The Pearson goodness of fit test requires grouping the sample points into a finite number of bins.  This grouping can be constructed in different ways, such as the
equally probable bins that we have used or even equidistant partition points.  Often, bins that accrue fewer than 5 data points are merged together to form a single bin.  In any case, the resulting test statistic and decision on acceptance or
rejection can be altered by various choices of grouping.

We shall now describe the Kolmogorov-Smirnov test which does not require any grouping.  This test can be used to test whether a data set comes from any specified continuous distribution.  We shall test the null hypothesis that the data set comes from a distribution with cumulative distribution function F[x].

The sample points must first be sorted into increasing order a1, . . . , an.  A sample distribution function G[x] is defined piecewise by G[x] = 0 for x < a1;
G[x] = i/n for ai  ≤ x  < ai+1 and i = 1, . . . , n-1; and G[x] = 1 for x ≥ an.  A test statistic is defined by 
         
                                                     D = n1/2 sup|F[x] - G[x]|,
                                                   
where the supremum is taken over all x in (-∞, ∞).

As the sample size n increases, then the test statistic converges to a distribution Y such that
                                                                     ∞
                                               P(Y > t) = 2  ∑ (-1)i+1 Exp[-2 i2 t2] 
                                                                    i=1
                                                                   
 [Degroot, 1986, p. 556].  Thus, the right tail value is defined by P(Y > D).  We reject
 the null hypothesis if this right tail value is less than the level of significance a.

As an example, we shall test our original data set of combined SAT scores for normality and for uniformity.  In each case though, the population parameters must be specified.  We first enter the data and obtain the sample statistics.

data={1270,1220,1160,950,1070,1110,1220,1150,1440,850,
     1280,1020,880,1080,1200,1120,970,1210,1030,1020,
     1400,720,1170,1070,860,920,1350,910,980,800,1070,
     1320,900,1200,1260,830,1330,1070,810,1240,820,
     1140,1170,870,1010,760,810,830,1240,910,1230,1100,
     970,1090,1300,1070,1040,1320,1230,890,960,920,
     1260,1060,1040,990,1040,1300,950,1090,1020,1270,
     1160,1050,1130,880,750,1200,1010,1060,1010,1080,
     1090,940,1170,920,990,1100,1270,930,790,1080,740,
     960,1170,1050,1250,1090,1010,1200,1080,1200,1070,
     970,1140,1360,1260,1170,1160,1310,1090,1040,1070,
     890,1280,1270,1190,1070,1410,1250,1070,910,1230};     
SampleSize=n=Length[data]
Xbar=N[Mean[data]]
S=N[StandardDeviation[data]]

We shall then test the null hypothesis that this data set comes a normally distributed population with mean 1079.11 and standard deviation 161.874.  To do so, we first sort the data into a list a, and define the piecewise continuous sampling distribution function, labeled here as h[x].  We then plot h[x] from the minimum a[[1]] of the data set  to the maximum a[[n]].

a=Sort[data];Clear[x];g[0,x_]=0;
Do[g[i_,x_]=If[a[[i]] <= x < a[[i+1]],i/n,0],{i,1,n-1}];
g[n,x_]=If[x>=a[[n]],1,0];
h[x_]=Sum[g[i,x],{i,0,n}];

b=Plot[h[x],{x,a[[1]],a[[n]]}]

We can then visually compare the sample distribution with the cumulative distribution function of the desired normal distribution.

X=NormalDistribution[Xbar,S];F[x_]=CDF[X,x];
c=Plot[F[x],{x,a[[1]],a[[n]]}]

Show[b,c]

Next, we must compute the maximum absolute difference between these two disributions.  We shall do so by graphing and inspection.

Plot[Abs[F[x]-h[x]],{x,a[[1]],a[[n]]}]

We can observe that the maximum is around .055.  If we desire, we can adjust the limits and re-plot to zoom in and get a more accurate value.  Using .055, we define the test statistic and compute the right tail value.

TestStatistic=d=N[.055*Sqrt[n]]

RightTailValue=
2 NSum[(-1)^(i+1) E^-(2 i^2 d^2),{i,1,Infinity}]

If the null hypothesis were true, then there would be an 85.0799% chance of the data being distributed as it is.  Obviously we accept normality since the tail value is so high.

Finally, we shall test for uniformity.  For the parameters we shall specify the minimum and maximum of the data set.  Thus, we are testing the null hypothesis that the data set comes from a population that is normally distributed over [a[[1]],a[[n]]]. We need only make this one change and re-run the program.

X=UniformDistribution[a[[1]],a[[n]]];F[x_]=CDF[X,x];
c=Plot[F[x],{x,a[[1]],a[[n]]}]

Show[b,c]

Plot[Abs[F[x]-h[x]],{x,a[[1]],a[[n]]}]

Here we are using a maximum absolute difference of .132.

TestStatistic=d=N[.132*Sqrt[n]]

RightTailValue=
2 NSum[(-1)^(i+1) E^-(2 i^2 d^2),{i,1,Infinity}]

If the population were uniform, then there would only be a 2.75% chance of the data being distributed as it is.  With this low tail value, we can reject the hypothesis that combined SAT scores are uniformly distributed.

CONCLUSION

Goodness of fit tests are an important part of data analysis; however, they are 
just too unwieldy to do by hand.  In particular, with standard textbook charts 
one usually cannot find the quantiles of the of the T distribution which are required to find the partition points in the Pearson test for normality.  Neither is the Kolmogorov-Smirnov test appealing to do by hand.  However with the above short programs, before one assumes a required condition of normality, one can use a data set to verify (or reject) such a condition.
    The Pearson test can easily be adjusted for simple hypothesis tests where the parameters are specified, or for composite tests where the parameters are not specified.  This Kolmogorov-Smirnov test though can only be used to test a specified distribution.  However it can be used for any continuous distribution with a known cumulative distribution function.  The Pearson test must be adjusted for each distribution since the ranges vary and thus the partitioning into bins must be adjusted.  The Pearson test is also be used to test discrete distributions, such as
Poisson, and for contingency table tests.  Though again each test is unique due to the required partitioning of the range.
     As a project, one can design a Pearson goodness of fit test for an exponential distribution with average value Ø and range [0, ∞).  If Ø is not specified, then 1/x- will be its maximum likelihood estimate.  The program can then be modified for the many continuous distributions having range [c, ∞) for some c.  One can then compare the results with a Kolmogorov-Smirnov test.

ACKNOWLEDGEMENTS 

The author would like to thank the referee for the many helpful suggestions for improving this article.

REFERENCES

L.J. Bain and M. Engelhardt. Introduction to Probability and Mathematical Statistics  Second Edition.  PWS-Kent, Boston (1992).

M. H. DeGroot.  Probability and Statistics  Second Edition. Addison-Wesley (1986).

R.V. Hogg and E.A. Tanis. Probability and Statistical Inference  Fourth Edition.  Macmillan Publishing Co., New York (1993).

ABOUT THE AUTHOR

David K. Neal is an associate professor of mathematics at Western Kentucky
University where he has taught for seven years.  His mathematical interests are
in pure and applied probability.  He enjoys enriching his courses with computer projects using Mathematica, which he now considers to be an indispensible part of both his teaching and research.  His other hobbies include music, sports, reading, and felines.

David K. Neal
Department of Mathematics
Western Kentucky University
Bowling Green, KY  42101
nealdk @wkuvx1.wku.edu