-
-
Notifications
You must be signed in to change notification settings - Fork 302
/
goodness.tex
2302 lines (2123 loc) · 90.2 KB
/
goodness.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[draft]{article}
\usepackage{makeidx}
%\title{\(\rightarrow\)DRAFT\(\leftarrow\)\\
%Goodness of Fit Techniques}
\title{Goodness of Fit Tests\\
{\large Documentation on {\tt libcdhc.a}}\\
{\large and}\\
{\large A GRASS Tutorial on {\tt s.normal}}}
\author{James Darrell McCauley\thanks{USDA National Needs Fellow,
Department of Agricultural Enginering, Purdue University. Email:
{\tt mccauley@ecn.purdue.edu}}}
\makeindex
\addtolength{\oddsidemargin}{-.55in}
\addtolength{\evensidemargin}{-.55in}
\addtolength{\textwidth}{.1in}
\addtolength{\marginparwidth}{.45in}
\addtolength{\topmargin}{-.25in}
\addtolength{\textheight}{.5in}
\def\libname{{\tt cdhc}}
\def\returns#1{\sffamily\slshape Returns \(\mathsf{#1}\).}
\def\function#1#2{\centerline{%
\protect\index{#1}
\framebox[.9\marginparwidth][l]{\vbox{\noindent\textsf{#2}}}}
\vspace{.5\baselineskip}}
%\def\function#1#2{\marginpar{%
% \protect\index{#1}
% \framebox[.9\marginparwidth][l]{\vbox{\textsf{#2}}}\hfill}}
\newenvironment{example}{%
\vspace{\baselineskip}
\par\noindent\hrulefill\par
\noindent{\em Example:}}{%
\par\noindent\hrulefill\par
\vspace{\baselineskip}}
\begin{document}
\bibliographystyle{plain}
\maketitle
\begin{abstract}
The methods used by the GRASS program {\tt s.normal}
are presented. These are various goodness of fit statistics for testing
the null hypothesis of normality. Other additional tests found in
\libname\, a C programming library,
are also documented (this document serves two puposes:
a tutorial for the GRASS geographic information system
and documentation for the library).
\end{abstract}
\section{Introduction}
This document is a programmer's
manual for \libname, a C programming library useful
for testing whether a sample is normally, lognormally,
or exponentially distributed.
Prototypes for library functions\footnote{%
Each function in the library returns a pointer to static double.
The \libname\ library was inspired by Johnson's
STATLIB collection of FORTRAN routines for testing
distribution assumptions~\protect\cite{johnson94}.
Some functions in \libname\
are loosely based on Johnson's work (they have been completely
rewritten, reducing memory requirements and number
of computations and fixing a few bugs). Others are based on
algorithms found in \emph{Applied Statistics}, \emph{Technometrics},
and other related journals.}
are given in the margins near
corresponding mathematical explanations. Hence, it is also
a user's guide for programs using \libname.
Readers should be equipped with at least one graduate course
in probability and statistics. Much of the background
and derivation/justification of each test has been
omitted. A good text for more background information
is {\em Goodness-of-Fit Techniques\/} by
D'Agostino and Stephens~\cite{dagostino86b} (see also references in text).
\subsection{Hypothesis Testing}
Before beginning the description of the tests, a few definitions
should be given. The general framework for mosts tests is that
the {\em null\/} hypothesis \(H_0\) is that a random variable \(x\)
follows a particular distribution \(F\left(x\right)\).
Generally, the {\em alternative\/} hypothesis is that
\(x\) does not follow \(F\left(x\right)\) (with no additional
usuable information; the Kotz Separate Families test in \S\ref{sec:kotz}
is one exception).
This may differ from the way that some have learned hypothesis testing
in that some tests are set up to reject the null hypothesis in
favor of the alternative.
A {\em simple\/} hypothesis implies that \(F\left(x\right)\)
is completely specified, e.g., \(x\sim N\left(0,1\right)\).
A {\em composite\/} hypothesis means that
one (or more) of the parameters of \(F\left(x\right)\)
is not completely specified, e.g., \(x\sim N\left(\mu,\sigma\right)\).
That is, the composite hypothesis may be:
\begin{displaymath}
H_0 : F\left(x\right) = F_0\left(x; \theta\right)
\end{displaymath}
where \(\theta=\left[\theta_1, \ldots,\theta_p\right]'\)
is a \(p\) vector of \emph{nuisance} parameters whose values
are unknown and must be estimated from data.
% Less is known
% about the theory of this later
% case, which is the most commonly encountered in practice.
\subsection{Probability Plots}
In addition to these analytical techniques, graphical
methods are valuable supplements. The most important
graphical technique is probability plotting. A \emph{probability plot}
\label{pplot}
is a plot of the cumulative distribution function \(F\left(x\right)\)
on the vertical axis versus \(x\) on the horizontal axis.
The vertical axis is scaled such that, if the data fit
the assumed distribution, the resulting plot will lie on
a straight line. Special plotting paper may be purchased
to do these plots; however, most modern scientific
plotting programs have this capability (e.g., {\tt gnuplot}).
Each test presented below
should be used in conjunction with a probability plot.
\subsection{Shape of Distributions}
Through much of the literature are references to Johnson
curves: \(S_U\) or \(S_B\) (see \S\ref{sec:johnson-su},
page~\pageref{sec:johnson-su}).
These refer to a system of distributions introduced by
Johnson~\cite{johnson49} where a standard normal random
variable \(Z\) is translated to \(\left(Z-\gamma\right)/\delta\)
and transformed using \(T\):
\begin{equation}
Y=T\left(\frac{Z-\gamma}{\delta}\right).
\end{equation}
Three families in Johnson's~\cite{johnson49} system are:
\begin{enumerate}
\item a family of bounded distributions, denoted by \(S_B\), where:
\begin{equation}
Y=T\left( \frac{e^x}{1+e^x} \right);
\end{equation}
\item a family lognormal distributions where:
\begin{equation}
Y=T\left( e^x \right);
\end{equation}
\item and a family of unbounded distributions, denoted by \(S_U\), where:
\begin{equation}
Y=\sinh\left(x\right) = T\left( e^x-e^{-x} \right).
\end{equation}
\end{enumerate}
In the \(S_B\) and \(S_U\) families, \(\gamma\) and \(\delta\)
govern the shape of the distribution. In the lognormal families,
\(\delta\) governs the shape while \(\gamma\) is only a scaling
factor~\cite{hoaglin85c}. Other approaches to exploring the
shape of a distribution include \(g\)- and
\(h\)-distributions~\cite{hoaglin85c} and
Pearson curves (see Bowman~\cite{bowman86}).
\subsection{Miscellaneous}
Many tests are presented here without mention of their relative
merits. Users are advised to consult the cited literature to
determine which test is appropriate for their situation. Sometimes
a certain test will have more \emph{power} than another; that is,
a test may have a better ability to reject a model when
the model is incorrect.
\section{Moments: \(b_2\) and \(\protect\sqrt{b_1}\)}
\function{omnibus\_moments(x,n)}
{double* \\
\hbox{omnibus\_moments(x,n)}\\
double *x;\\
int n;\\
\returns{\left[\sqrt{b_1},b_2\right]'}}%
Let \(x_1, x_2, \ldots, x_n\) be the \(n\)
observations with mean:
\begin{equation}
m_1 = \frac{1}{n}\sum_{j=1}{n} x_j.
\end{equation}
The central moments are defined as:
\begin{equation}
\label{eqn:moments}
m_i = \frac{1}{n}\sum_{j=1}{n}\left( x_j - m_i\right)^i,\: i=2,3,4.
\end{equation}
The sample skewness \(\left(\sqrt{b_1}\right)\)
and kurtosis \(\left(b_2\right)\) are defined as:
\begin{equation}
\sqrt{b_1} = m_3/m_2^{3/2} = \sqrt{n}
\left(\sum_{j=1}^n\left(x_i-\bar{x}\right)^3\right)/
\left( \sum_{j=1}^n\left(x_i-\bar{x}\right)^2 \right)^{3/2}
\end{equation}
and
\begin{equation}
\label{eqn:4th-sample-moment}
b_2 = m_4/m_2^2.
\end{equation}
These are invariant under both origin and scale changes~\cite{bowman86}.
When a distribution is specified, these are denoted as
\(\sqrt{\beta_1}\) and \(\beta_2\).
For a standard normal, \(\sqrt{\beta_1}=0\) and \(\beta_2=3\).
To use either or both of these statistics to test for
departure from normality, these are sometimes transformed
to their standardized to their normal equivalent
deviates, \(X\left(\sqrt{b_1}\right)\) and \(X\left(b_1\right)\).
For \(X\left(\sqrt{b_1}\right)\), D'Agostino and
Pearson~\cite{dagostino73} gave coefficients \(\delta\)
and \(\lambda\) (\(n=8\) to 1000) for:
\begin{equation}
X\left(\sqrt{b_1}\right) = \delta \sinh^{-1}
\left(\sqrt{b_1}/\lambda\right)
\end{equation}
that transforms \(\sqrt{b_1}\) to a standard normal
using a Johnson \(S_U\) approximation (Table~\ref{tbl:johnson}).
\label{sec:johnson-su}
An equivalent approximation~\cite{dagostino86}
that avoids the use of tables is given by:
\begin{enumerate}
\item Compute \(\sqrt{b_1}\) from the sample data.
\item Compute:
\begin{eqnarray}
Y &=& \sqrt{b_1} \left[\frac{\left(n+1\right)\left(n+3\right)}
{6\left(n-2\right)}\right]^{\frac{1}{2}}, \\
\beta_2 &=& \frac{3\left(n^2+27n-70\right)\left(n+1\right)\left(n+3\right)}
{\left(n-2\right)\left(n+5\right)\left(n+7\right)\left(n+9\right)},\\
W^2 &=& \sqrt{2\left(\beta_2-1\right)}-1, \\
\delta &=& 1/\sqrt{\log W}, \\
\mbox{and}\\
\alpha &=& \sqrt{2/\left(W^2-1\right)}.
\end{eqnarray}
\item Compute the standard normal variable:
\begin{equation}
Z = \delta \log\left[Y/\alpha + \sqrt{\left(Y/\alpha\right)^2+1}\,\right].
\end{equation}
\end{enumerate}
This procedure is applicable for \(n\ge8\).
%D'Agostino~\cite{dagostino86} also notes
%that the normal approximation given by
%\begin{equation}
%\sqrt{\beta_1}\left[\frac{\left(n+1\right)\left(n+3\right)}
%{6\left(n-2\right)}\right]^{\frac{1}{2}}
%\end{equation}
%is valid for \(n\ge150\)~\cite{dagostino86}.
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(\sqrt{b_1} = 0.2373\). Suppose that we wish to test the
hypothesis of normality:
\(H_0\): \(\sqrt{\beta_1}=0\) (normality)
\noindent versus the two-sided alternative
\(H_1\): \(\sqrt{\beta_1}\ne0\) (non-normality)
\noindent at a level of significance of 0.05.
Following the procedure given above,
\(Y =2.3454\),
\(\beta_2 = 3.0592\),
\(W^2 = 1.0294\),
\(\delta=12.6132\),
\(\alpha=8.2522\), and
\(Z=1.5367\).
At a 0.05 significance level for a two-sided test, we reject
the null hypothesis of normality if \(\left|Z\right|\ge1.96\). In
this instance, we cannot reject \(H_0\).
\end{example}
The fourth standardized moment \(b_2\) may be used to
test the normality hypothesis by the following
procedure~\cite{anscombe63}:
\begin{enumerate}
\item Compute \(b_2\) from the sample data.
\item Compute the mean and variance of \(b_2\):
\begin{equation}
E\left(b_2\right) = \frac{3\left(n-1\right)}{n+1}
\end{equation}
and
\begin{equation}
Var\left(b_2\right) = \frac{24n\left(n-2\right)\left(n-3\right)}
{\left(n+1\right)^2\left(n+3\right)\left(n+5\right)}.
\end{equation}
\item Compute the standardized value of \(b_2\):
\begin{equation}
y = \frac{b_2-E\left(b_2\right)}{Var\left(b_2\right)}.
\end{equation}
\item Compute the third standardized moment of \(b_2\):
\begin{equation}
\sqrt{\beta_1\left(b_2\right)} =
\frac{6\left(n^2-5n+2\right)}{\left(n+7\right)\left(n+9\right)}
\sqrt{\frac{6\left(n+3\right)\left(n+5\right)}
{n\left(n-2\right)\left(n-3\right)}}.
\end{equation}
\item Compute:
\begin{equation}
A=6+\frac{8}{\sqrt{\beta_1\left(b_2\right)}}\left[
\frac{2}{\sqrt{\beta_1\left(b_2\right)}} +
\sqrt{1+\frac{4}{\sqrt{\beta_1\left(b_2\right)}}}\,\right].
\end{equation}
\item Compute:
\begin{equation}
\label{eqn:z-b2}
Z = \left(\left(1-\frac{2}{9A}\right)-
\left[\frac{1-2/A}{1+y\sqrt{2/\left(A-4\right)}}\right]^{\frac{1}{3}}\right)/
\sqrt{2/\left(9A\right)}
\end{equation}
where \(Z\) is a standard normal variable with
zero mean and variance of one.
\end{enumerate}
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(b_2 =1.9148\). Suppose that we wish to test the
hypothesis of normality:
\(H_0\): \(\beta_2=3\) (normality)
\noindent versus the one-sided alternative
\(H_1\): \(\beta_2>3\) (non-normality)
\noindent at a level of significance of 0.05. We would
reject \(H_0\) if \(Z\) (eqn.~\ref{eqn:z-b2}) is larger
than 1.645 (Table~\ref{tbl:normal}). Following the procedure given above,
\(E\left(b_2\right)=2.9897\),
\(Var\left(b_2\right)=0.0401\),
\(y=-26.8366\),
\(\sqrt{\beta_1\left(b_2\right)}=0.0989\),
\(A=2163\), and
\(Z=-131.7\).
Therefore, we cannot reject \(H_0\).
\end{example}
\subsection{Omnibus Tests for Normality}
\section{Geary's Test of Normality}
\label{sec:geary}
\function{geary\_test(x,n)}
{double*\\
\hbox{geary\_test(x,n)}\\
double *x;\\
int n;\\
\returns{\left[\sqrt{a},y\right]'}}
Let \(x_1, x_2, \ldots, x_n\) be the \(n\)
observations. The ratio of the mean deviation
to the standard deviation is given as:
\begin{equation}\label{eqn:geary}
a = \frac{1}{n\sqrt{m_2}}\sum_{j=1}^n \left|x_i-\bar{x}\right|
\end{equation}
where \(\bar{x}=\sum_{i=1}^n x_i\) and \(m_2\) is defined
by eqn.~\ref{eqn:moments}.
This ratio can be transformed
a standard normal~\cite{dagostino86} via
\begin{equation}\label{eqn:geary-normal}
y = \frac{\sqrt{n}\left(a-0.7979\right)}{0.2123}.
\end{equation}
This test is valid for \(n\ge41\).
More generally, Geary~\cite{geary47} considered tests of the
form
\begin{equation}
a\left(c\right) =
\frac{1}{nm_2^{c/2}}
\sum_{j=1}^n \left|x_i-\bar{x} \right|^c \: \mbox{for}\: c\ge1
\end{equation}
where \(a\left(1\right)=a\) of eqn.~\ref{eqn:geary}, and
\(a\left(4\right)=b_2\) of eqn.~\ref{eqn:4th-sample-moment}.
D'Agostino and Rosman~\cite{dagostino74} conclude that
Geary's \(a\) test has good power for symmetric alternatives
and skewed alternatives with \(\beta_2 < 3\) when compared to
other tests, though for symmetric alternatives, \(b_2\)
(eqn.~\ref{eqn:4th-sample-moment}) can sometimes be more powerful and
for skewed alternatives, \(W\) (eqn~\ref{eqn:w-test})
or \(W'\) (eqn~\ref{eqn:w-prime-test})
usually dominate \(a\).
The Geary test (eqns.~\ref{eqn:geary}-\ref{eqn:geary-normal})
is seldom used today---D'Agostino~\cite{dagostino86} include it
in his summary work because it is of ``historical interest.''
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(a = 0.8823\). Suppose that we wish to test the
hypothesis of normality:
\(H_0\): normality
\noindent versus the two-sided alternative
\(H_1\): non-normality
\noindent at a level of significance of 0.05.
From eqn.~\ref{eqn:geary-normal}, \(y=9.9607\).
\end{example}
\section{Extreme Normal Deviates}
\label{sec:extreme}
\function{extremes(x,n)}
{double* \\
\hbox{extremes(x,n)}\\
double *x;\\
int n;\\
\returns{\left[x_n-\bar{x}, x_1-\bar{x}\right]'}}
Let \(x_1 \le x_2 \le \cdots \le x_n\) be the \(n\)
observations. Given a known normal deviation \(\sigma\),
the largest and smallest deviation from a normal population
may be computed:
\begin{equation}
u_n = \frac{x_n-\bar{x}}{\sigma}
\end{equation}
and
\begin{equation}
u_1 = -\frac{x_1-\bar{x}}{\sigma},
\end{equation}
respectively. These statistics are potentially
useful for detecting outliers for populations
with a known \(\sigma\) but an unknown mean.
Table 25 in Pearson and Hartley~\cite{pearson76}
gives percentage points for this statistic.
Pearson and Hartley~\cite{pearson76} also give examples
of the use of extreme deviates when an estimator of
\(\sigma\) (independent of the sample) is
known and when a combined ``internal''
and ``external'' estimate is used.
\section{EDF Statistics for Testing Normality}
[Note: This section follows closely the presentation of
Stephens~\cite{stephens86}.]
Let \(x_1 \le x_2 \le \cdots \le x_n\) be the \(n\)
observations. Suppose that the continuous distribution of \(x\)
is \(F\left(x\right)\). The empirical distribution function (EDF)
is \(F_n\left(x\right)\) defined by:
\begin{equation}
F_n\left(x\right) = \frac{1}{n}\left(\mbox{number of observations}
\le x\right); \: -\infty < x < \infty
\end{equation}
or
\begin{displaymath}
\begin{array}{rclll}
F_n\left(x\right)& = &0, & x<x_1\\
F_n\left(x\right)& = &\frac{1}{n}, & x_i\le x<x_{i+1}, & i=1,\ldots,n-1\\
F_n\left(x\right)& = &1, & x_n\le x.
\end{array}
\end{displaymath}
Thus \(F_n\left(x\right)\) is a step function calculated from
the data. As \(n\rightarrow\infty\), \\
\(\left|F_n\left(x\right)- F\left(x\right)\right|\)
decreases to zero with probability one~\cite{stephens86}.
EDF statistics that measure the difference between
\(F_n\left(x\right)\) and \(F\left(x\right)\) are divided
into two classes: supremum and quadratic.
On a graph of
\(F_n\left(x\right)\) and \(F\left(x\right)\) versus \(x_i\),
denote the largest vertical distance when
\(F_n\left(x\right)>F\left(x\right)\) as \(D^+\).
Also, let \(D^-\) denote the largest vertical distance when
\(F_n\left(x\right)<F\left(x\right)\). These two
measures are supremum statistics.
Quadratic statistics are given by the Cram\'er--von Mises family
\begin{equation}
\label{eqn:cramer-family}
Q = n\int_{-\infty}^{\infty}
\left(F_n\left(x\right) - F\left(x\right)\right)^2
\psi\left(x\right) d F\left(x\right)
\end{equation}
where \(\psi\left(x\right)\) is a weighting function~\cite{stephens86}.
To compute these statistics, the Probability Integral Transformation
is used: \(z=F\left(x\right)\) where \(F\left(x\right)\) is
the Gaussian distribution. The new variable, \(z\), is uniformly
distributed between 0 and 1. Then \(z\) has distribution
function \(F^*\left(z\right)=z\), \(0\le z\le1\).
A sample \(x_1, x_2, \ldots, x_n\) gives values \(z_i=F\left(x_i\right)\),
\(i=1, \ldots, n\), and \(F^*_n\left(z\right)\) is the EDF of
values \(z_i\). For testing normality,
\begin{equation}
z_{\left(i\right)} = \Phi\left(
\left(x_{\left(i\right)}-\hat{\mu}\right)/\hat{\sigma}
\right)
\end{equation}
where \(\hat{\mu}\) and \(\hat{\sigma}\) are estimated from
the data and \(\Phi\left(\cdot\right)\) denotes the cumulative
probability of a standard normal. For testing if the data
follows an exponential distribution \(\mbox{Exp}\left(\alpha,\beta\right)\),
where \(\alpha\) is known to be zero, \(\hat{\beta}\)
is estimated by \(\bar{x}\) (the sample mean) and
\begin{equation}
z_{\left(i\right)} = 1-\exp\left(-x_{\left(i\right)}/\bar{x}\right).
\end{equation}
Now, EDF statistics can be computed by comparing \(F^*_n\left(z\right)\)
and a uniform distribution for \(z\). These take the same values
as comparisons between \(F_n\left(x\right)\) and \(F\left(x\right)\):
\begin{equation}
F_n\left(x\right) - F\left(x\right) =
F^*_n\left(z\right) - F^*\left(z\right) =
F^*_n\left(z\right) - z.
\end{equation}
After ordering \(z\)-values,
\(z_{\left(1\right)}\le
z_{\left(2\right)} \le \cdots
\le z_{\left(n\right)}\) and computing \(\bar{z}=\sum_{i=1}^n z_i/n\),
the supremum statistics are
\begin{equation}
\label{eqn:dplus}
D^+=\max_{i=1,\ldots,n}\left(i/n-z_{\left(i\right)}\right)
\end{equation}
and
\begin{equation}
\label{eqn:dminus}
D^-=\max_{i=1,\ldots,n}\left(z_{\left(i\right)}-\left(i-1\right)/n\right).
\end{equation}
\subsection{Kolmogorov \(D\)}
\function{kolmogorov\_smirnov(x,n)}
{double* \\
\hbox{kolmogorov\_smirnov(x,n)}\\
double *x;\\
int n;\\
\returns{\left[D^n,D\right]'}}
\function{kolmogorov\_smirnov\_exp(x,n)}
{double* \\
\hbox{kolmogorov\_smirnov\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[D^e,D\right]'}}
The most well-known EDF statistic is Kolmogorov's \(D\), computed
from supremum statistics:
\begin{equation}
D = \sup_x\left|F_n\left(x\right) - F\left(x\right)\right| =
\max\left(D^+,D^-\right).
\end{equation}
The modified form for testing a completely specified
distribution~\cite{stephens86}:
\begin{equation}
D^*=D\left(\sqrt{n}+0.12+0.11/\sqrt{n}\right).
\end{equation}
For testing a normal distribution with \(\mu\) and \(\sigma\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
D^n=D\left(\sqrt{n}-0.01+0.85/\sqrt{n}\right).
\end{equation}
For testing an exponential distribution with \(\alpha\) and \(\beta\)
% origin and scale
unknown, \(D\) does not need modified~\cite{stephens86}.
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(D^n = 4.0314\) and \(y=\). Suppose that we wish to test the
hypothesis of normality:
\(H_0\): normality
\noindent versus the two-sided alternative
\(H_1\): non-normality
\noindent at a level of significance of 0.05.
\end{example}
\subsection{Kuiper's \(V\)}
\label{sec:kuiper}
\function{kuipers\_v(x,n)}
{double* \\
\hbox{kuipers\_v(x,n)}\\
double *x;\\
int n;\\
\returns{\left[V^n,V\right]'}}
\function{kuipers\_v\_exp(x,n)}
{double* \\
\hbox{kuipers\_v\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[V^e,V\right]'}}
Kuiper's~\cite{kuiper60} \(V\) is another statistic computed
from supremum statistics:
\begin{equation}
\label{eqn:kuipers-v}
V = D^+ + D^-.
\end{equation}
The modified form for testing a completely specified
distribution~\cite{stephens86}:
\begin{equation}
V^*=V\left(\sqrt{n}+0.155 +0.24\sqrt{n}\right).
\end{equation}
For testing a normal distribution with \(\mu\) and \(\sigma\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
V^n=V\left(\sqrt{n}+0.05+0.82/\sqrt{n}\right).
\end{equation}
For testing an exponential distribution with \(\alpha\) and \(\beta\)
unknown, \(V\) the modified equation is~\cite{stephens86}:
\begin{equation}
V^e=\left(V-0.2/\sqrt{n}\right)
\left(\sqrt{n}+0.24+0.35/\sqrt{n}\right).
\end{equation}
\subsection{Pyke's Statistics}
\label{sec:pyke}
For some purposes, eqns.~\ref{eqn:dplus} and~\ref{eqn:dminus}
may be modified to~\cite{pyke59}:
\begin{equation}
\label{eqn:cplus}
C^+=\max_{0\le i\le n}\left(\frac{i}{n+1}-z_{\left(i\right)}\right),\:
z_{\left(0\right)}=0,
\end{equation}
and
\begin{equation}
\label{eqn:cminus}
C^-=\max_{0\le i\le n}\left(z_{\left(i\right)}-\frac{i}{n+1}\right)
\end{equation}
(following the modification of notation by Durbin~\cite{durbin73}). Then,
\begin{equation}
C = \max\left(C^+,C^-\right).
\end{equation}
Durbin~\cite{durbin73} notes that these modifications to
eqns.~\ref{eqn:dplus} and~\ref{eqn:dminus} are related to
the fact that \(E\left(z_{\left(i\right)}\right)=i/\left(n+1\right)\).
Percentage points were given by Durbin~\cite{durbin69}.
\subsection{Brunk's \(B\)}
\label{sec:brunk}
As an alternative to Kuiper's \(V\) (eqn.~\ref{eqn:kuipers-v}),
Brunk~\cite{brunk62} suggests:
\begin{equation}
\label{eqn:brunks-b}
B = C^+ + C^-
\end{equation}
where \(C^+\) and \(C^-\) are given by eqns.~\ref{eqn:cplus}
and \ref{eqn:cminus}.
\subsection{Cram\'er--von Mises \(W^2\)}
\label{sec:cramer-von-mises}
\function{cramer\_von\_mises(x,n)}
{double* \\
\hbox{cramer\_von\_mises(x,n)}\\
double *x;\\
int n;\\
\returns{\left[W^{2,n},W^2\right]'}}
\function{cramer\_von\_mises(x,n)}
{double* \\
\hbox{cramer\_von\_mises\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[W^{2,e},W^2\right]'}}
Quadratic statistics are computed from
the Cram\'er--von Mises family given in eqn~\ref{eqn:cramer-family}.
When \(\psi\left(x\right)=1\) in eqn~\ref{eqn:cramer-family}, the statistic is
the Cram\'er--von Mises statistic \(W^2\):
\begin{equation}
W^2=\sum_{j=1}^n\left(Z_i - \left(2j-1\right)/\left(2n\right)\right)^2
+\frac{1}{12n}
\end{equation}
(When \(\psi\left(x\right)=
\left(F\left(x\right)\left(1 - F\left(x\right)\right)\right)^{-1}\),
this yields the Anderson--Darling statistic given below
in \S\ref{sec:anderson-darling}~\cite{stephens86}.)
The modified form for testing a completely specified
distribution~\cite{stephens86}:
\begin{equation}
W^{2,*} = \left(W^2-0.4/n +0.6/n^2\right)/\left(1 + 1/n\right).
\end{equation}
For testing a normal distribution with \(\mu\) and \(\sigma\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
W^{2,n}=W^2\left(1.0 + 0.5/n\right).
\end{equation}
For testing an exponential distribution with \(\alpha\) and \(\beta\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
W^{2,e}=W^2\left(1.0 + 2.8/n -3/n^2\right).
\end{equation}
\subsection{Watson \(U^2\)}
\label{sec:watson}
\function{watson\_u2(x,n)}
{double* \\
\hbox{watson\_u2(x,n)}\\
double *x;\\
int n;\\
\returns{\left[U^{2,n}, U^{2}\right]'}}
\function{watson\_u2\_exp(x,n)}
{double* \\
\hbox{watson\_u2\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[U^{2,e}, U^{2}\right]'}}
\begin{equation}
U^2=W^2-n\left(\bar{Z}-0.5\right)^2
\end{equation}
where \(W^2\) is the Cram\'er--von Mises statistic
(\S\ref{sec:cramer-von-mises}).
The modified form for testing a completely specified
distribution~\cite{stephens86}:
\begin{equation}
U^{2,*} = \left(U^2-0.1/n +0.1/n^2\right)/\left(1 + 0.8/n\right).
\end{equation}
For testing a normal distribution with \(\mu\) and \(\sigma\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
U^{2,n}=U^2\left(1.0 + 0.5/n\right).
\end{equation}
For testing an exponential distribution with \(\alpha\) and \(\beta\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
U^{2,e}=U^2\left(1.0 + 2.3/n -3/n^2\right).
\end{equation}
\subsection{Anderson--Darling \(A^2\)}
\label{sec:anderson-darling}
\function{anderson\_darling(x,n)}
{double* \\
\hbox{anderson\_darling(x,n)}\\
double *x;\\
int n;\\
\returns{\left[A^{2,n}, A^{2}\right]'}}
\function{anderson\_darling\_exp(x,n)}
{double* \\
\hbox{anderson\_darling\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[A^{2,e}, A^{2}\right]'}}
Anderson and Darling~\cite{anderson54} present
another EDF test statistic which is sensitive at the
tails of the distribution (rather than near
the median).
When \(\psi\left(x\right)=
\left(F\left(x\right)\left(1 - F\left(x\right)\right)\right)^{-1}\)
in eqn.(\ref{eqn:cramer-family}),
this yields the Anderson--Darling statistic~\cite{anderson54,stephens86}:
\begin{equation}
A^2 = -n - \frac{1}{n} \sum_{j=1}^n \left(2j-1\right)
\left[ \ln z_j + \ln\left(1-z_{n-j+1}\right)\right].
\end{equation}
Equivalently~\cite{stephens86},
\begin{equation}
A^2 = -n - \frac{1}{n} \sum_{j=1}^n\left[ \left(2j-1\right)
\ln z_j + \left(2n+1-2j\right) \ln\left(1-z_{j}\right)\right].
\end{equation}
Anderson and Darling~\cite{anderson54} give
the following asymptotic significance values of \(A^2\):
\begin{center}
\begin{tabular}{cc}\hline
Significance & Significance \\
Level & Point \\ \hline \hline
0.10 & 1.933\\
0.05 & 2.492\\
0.01 & 3.857\\ \hline
\end{tabular}
\end{center}
Anderson and Darling~\cite{anderson54} state that
sample size should be at least 40; however, Stephens~\cite{stephens86}
give the same asymptotic values (for more significance levels)
for a sample size \(\ge5\).
For testing a completely specified distribution, \(A^2\)
is used unmodified.
For testing a normal distribution with \(\mu\) and \(\sigma\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
A^{2,n}=A^2\left(1.0 + 0.75/n+2.25/n^2\right).
\end{equation}
For testing an exponential distribution with \(\alpha\) and \(\beta\)
unknown, the modified equation is~\cite{stephens86}:
\begin{equation}
A^{2,e}=A^2\left(1.0 + 5.4/n -11/n^2\right).
\end{equation}
\subsection{Durbin's Exact Test}
\label{sec:durbin}
\function{durbins\_exact(x,n)}
{double* \\
\hbox{durbins\_exact(x,n)}\\
double *x;\\
int n;\\
\returns{\left[K_m,\sqrt{n}K_m\right]'}}
Durbin~\cite{durbin61} presented a modified
Kolmogorov test. The discussion that follows
has been adapted from Durbin's work~\cite{durbin61}.
Let \(x_1, x_2, \ldots, x_n\) be the \(n\)
i.i.d.\ observations and suppose that
it is desired to test the hypothesis that
they come from the continuous distribution \(F\left(x\right)\).
If the null hypothesis is true, then \(u_j=F\left(x_j\right)\)
(\(j=1,\ldots,n)\) are independent \(U\left(0,1\right)\)
variables and are randomly scattered on the (0,1) interval.
Clustering may indicated a departure from the null hypothesis.
Denoting the ordered \(u\)'s by
\(0 \le u_{\left(1\right)} \le \cdots \le u_{\left(n\right)} \le 1\),
let \(c_1=u_{\left(1\right)}\),
\(c_2=u_{\left(j\right)}-u_{\left(j-1\right)}\)
(\(j=2,\ldots,n\)), and \(c_{n+1}=1-u_{\left(n\right)}\).
Since the interest is in relative magnitudes of \(c\)'s, these
are ordered:
\(c_{\left(1\right)} \le c_{\left(2\right)} \cdots
\le c_{\left(n\right)}\). Then, the following transformation
is applied:
\begin{equation}
\label{eqn:durbin:g}
g_j=\left(n+2-j\right)\left(c_{\left(j\right)}
-c_{\left(j-1\right)}\right)\:
\left(c_{\left(0\right)}=0;\: j=1,\ldots,n+1\right).
\end{equation}
Durbin~\cite{durbin61} shows that the \(g\)'s, which depend
on the {\em ordered\/} intervals, have the same
distribution as the {\em unordered\/} \(c\)'s.
Letting
\begin{equation}
\label{eqn:durbin:w}
w_r = \sum_{j=1}^r g_j
\end{equation}
it follows that \(w_1, \ldots, w_n\) have the same distribution
as the ordered \(U\left(0,1\right)\) variables
\(u_{\left(1\right)}, \ldots, u_{\left(n\right)}\).
From eqns.~\ref{eqn:durbin:g} and \ref{eqn:durbin:w}, \(w_j\)
can be expressed as:
\begin{equation}
w_j=c_{\left(1\right)} + \cdots
+ c_{\left(j-1\right)} + \left(n+2-j\right)c_{\left(j\right)},\:
\left(j=1,\ldots,n\right),
\end{equation}
where
\(c_{\left(1\right)} \le \cdots \le c_{\left(n\right)}\) is
the ordered set of intervals.
In addition to two other test, Durbin~\cite{durbin61} introduces
the {\em modified Kolmogorov test}. The test statistic is:
\begin{equation}
K_m = \max_{r=1,\ldots,n}\left(\frac{r}{n}-w_r\right).
\end{equation}
The test procedure is to reject when \(K_m\) is greater
than the value tabulated for a one-sided Kolmogorov test.
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(K_m = 0.4127\). To test the
hypothesis of normality:
\(H_0\): normality
\noindent versus the one-sided alternative
\(H_1\): non-normality
\noindent at a level of significance of 0.05, we would
reject \(H_0\) if \(K_m\) is larger
than 0.895 (critical value of \(D\) for \(\alpha=0.05\).
Therefore, we cannot reject \(H_0\).
\end{example}
\section{Chi-Square Test}
\function{chi\_square(x,n)}
{double* \\
\hbox{chi\_square(x,n)}\\
double *x;\\
int n;\\
\returns{\left[x^2,k-3\right]'}}
\function{chi\_square\_exp(x,n)}
{double* \\
\hbox{chi\_square\_exp(x,n)}\\
double *x;\\
int n;\\
\returns{\left[x^2,k-2\right]'}}
According to Shapiro~\cite{shapiro90},
the chi-square goodness of fit test is the oldest
procedure for testing distributional assumptions.
It is useful for testing normality and exponentiality
when the number of observations is large (because its power
is poor for small samples when compared to other tests).
It is also useful when data are discrete~\cite{shapiro90}.
The basic idea is to divide the \(n\) data into \(k\) cells
and compare the observed number in each cell with the
expected number in each cell. The resulting statistic
is distributed as a chi-square random variable with
\(k-1-t\) degrees of freedom, where \(t\) is the number
of parameters estimated. The number of cells is taken
as
\begin{equation}
k=\mbox{(int)} 4\left[0.75\left(n-1\right)^2\right]^{1/5}.
\end{equation}
\marginpar{what should the notation be for rounding? For ceil,
we use \(\lceil x\rceil\). For floor, we use \(\lfloor x\rfloor\).}
The ratio \(n/k\) should be at least 5; otherwise another
test should be used~\cite{shapiro90}. In this implementation,
\(k\) is decremented by one until \(n/k\ge5\).
Let \(x_{\left(1\right)},
x_{\left(2\right)},\ldots, x_{\left(k\right)}\)
be the upper boundaries of cells. Choose \(x_{\left(i\right)}\)
so that the probability of being in any cell
is the same:
\begin{equation}
P\left(x\le x_{\left(i\right)}\right) = \frac{i}{k},\:
i=1,2,\ldots,k
\end{equation}
In thse implmentation, only the case of raw data, as opposed
to pre-tabulated data, is considered (i.e., equal probability cells).
For testing the normality hypothesis,
let \(x_{\left(0\right)}=-\infty\) and
\(x_{\left(k\right)}=\infty\).
The values of \(x_{\left(i\right)}\) are:
\begin{equation}
x_{\left(i\right)} = \bar{x} + s\,Z_{i/k}
\end{equation}
where \(\bar{x}\) and \(s\) are estimated
mean and variance parameters and \(Z_{i/k}\)
are percentiles of the standard normal distribution.
The test statistic is
\begin{equation}
\label{eqn:chi-square}
x^2 = \frac{k}{n}\sum_{i=1}^k f_i^2-n
\end{equation}
where \(f_i\) is the number of observations in cell \(i\).
The hypothesis of normality is rejected at an \(\alpha\)
level if \(x^2\) is greater \(x^2_{\alpha}\), a
\(\chi^2\) random variable with \(k-3\) degrees of freedom.
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(x^2 = 952.7\) with \(\nu=45\) degrees of freedom.
Since \(\chi^2_{45,0.05}\approx30.33\) (Table~\ref{tbl:chisq}),
we reject \(H_0\) at an \(\alpha=0.05\) level.
\end{example}
For testing the exponentiality hypothesis,
let \(x_{\left(0\right)}=0\) and
\(x_{\left(k\right)}=\infty\).
The values of \(x_{\left(i\right)}\) are:
\begin{equation}
x_{\left(i\right)} = -\frac{1}{\lambda}\ln\left(1-\frac{i}{k}\right),
i=1,2,\ldots,k-1.
\end{equation}
The parameter \(\lambda\) is estimated from
\begin{equation}
\hat{\lambda} = n \left(\sum_{i=1}^n x_i\right)^{-1}
\end{equation}
where \(x_i\) is the \(i\)th observation in
the sample. Equation~(\ref{eqn:chi-square})
is the statistic used for testing exponentiality. The hypothesis
of exponentiality is rejected at an \(\alpha\) level if
\(x^2\), a \(\chi^2\) random variable with \(k-2\)
degrees of freedom.
\begin{example}
For the sample data given in Table~\ref{tbl:pine} (\(n=584\)),
\(x^2 = 308.11\) with \(\nu=46\) degrees of freedom.
Since \(\chi^2_{46,0.05}\approx31.16\) (Table~\ref{tbl:chisq}),
we reject \(H_0:\) exponentiality, at an \(\alpha=0.05\) level.
\end{example}
\section{Analysis of Variance Tests}