c14-3 (Numerical Recipes in C), страница 2

In such cases, the list of data points canbe easily converted to an unbiased estimator SN (x) of the cumulative distributionfunction of the probability distribution from which it was drawn: If the N events arelocated at values xi , i = 1, . . . , N , then SN (x) is the function giving the fractionof data points to the left of a given value x. This function is obviously constantbetween consecutive (i.e., sorted into ascending order) xi ’s, and jumps by the sameconstant 1/N at each xi . (See Figure 14.3.1.)Different distribution functions, or sets of data, give different cumulativedistribution function estimates by the above procedure. However, all cumulativedistribution functions agree at the smallest allowable value of x (where they arezero), and at the largest allowable value of x (where they are unity).

(The smallestand largest values might of course be ±∞.) So it is the behavior between the largestand smallest values that distinguishes distributions.One can think of any number of statistics to measure the overall differencebetween two cumulative distribution functions: the absolute value of the area betweenthem, for example.

Or their integrated mean square difference. The KolmogorovSmirnov D is a particularly simple measure: It is defined as the maximum valueof the absolute difference between two cumulative distribution functions. Thus,for comparing one data set’s SN (x) to a known cumulative distribution functionP (x), the K–S statistic isD=max−∞<x<∞|SN (x) − P (x)|(14.3.5)while for comparing two different cumulative distribution functions SN1 (x) andSN2 (x), the K–S statistic isD=max−∞<x<∞|SN1 (x) − SN2 (x)|(14.3.6)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Ri + Sii624Statistical Description of DataSN (x)DP(x)xxFigure 14.3.1.Kolmogorov-Smirnov statistic D. A measured distribution of values in x (shownas N dots on the lower abscissa) is to be compared with a theoretical distribution whose cumulativeprobability distribution is plotted as P (x). A step-function cumulative probability distribution SN (x) isconstructed, one that rises an equal amount at each measured point. D is the greatest distance betweenthe two cumulative distributions.What makes the K–S statistic useful is that its distribution in the case of the nullhypothesis (data sets drawn from the same distribution) can be calculated, at least touseful approximation, thus giving the significance of any observed nonzero value ofD.

A central feature of the K–S test is that it is invariant under reparametrizationof x; in other words, you can locally slide or stretch the x axis in Figure 14.3.1,and the maximum distance D remains unchanged. For example, you will get thesame significance using x as using log x.The function that enters into the calculation of the significance can be writtenas the following sum:QKS (λ) = 2∞X(−1)j−1 e−2j2λ2(14.3.7)j=1which is a monotonic function with the limiting valuesQKS (0) = 1QKS (∞) = 0(14.3.8)In terms of this function, the significance level of an observed value of D (asa disproof of the null hypothesis that the distributions are the same) is givenapproximately [1] by the formulahpp i Ne + 0.12 + 0.11/ Ne DProbability (D > observed ) = QKS(14.3.9)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).cumulative probability distributionChapter 14.14.3 Are Two Distributions Different?625where Ne is the effective number of data points, Ne = N for the case (14.3.5)of one distribution, andNe =N1 N2N1 + N2(14.3.10)#include <math.h>#include "nrutil.h"void ksone(float data[], unsigned long n, float (*func)(float), float *d,float *prob)Given an array data[1..n], and given a user-supplied function of a single variable func whichis a cumulative distribution function ranging from 0 (for smallest values of its argument) to 1(for largest values of its argument), this routine returns the K–S statistic d, and the significancelevel prob.

Small values of prob show that the cumulative distribution function of data issignificantly different from func. The array data is modified by being sorted into ascendingorder.{float probks(float alam);void sort(unsigned long n, float arr[]);unsigned long j;float dt,en,ff,fn,fo=0.0;sort(n,data);en=n;*d=0.0;for (j=1;j<=n;j++) {fn=j/en;ff=(*func)(data[j]);dt=FMAX(fabs(fo-ff),fabs(fn-ff));if (dt > *d) *d=dt;fo=fn;}en=sqrt(en);*prob=probks((en+0.12+0.11/en)*(*d));If the data are already sorted into ascending order, then this call can beomitted.Loop over the sorted data points.Data’s c.d.f.

after this step.Compare to the user-supplied function.Maximum distance.Compute significance.}#include <math.h>void kstwo(float data1[], unsigned long n1, float data2[], unsigned long n2,float *d, float *prob)Given an array data1[1..n1], and an array data2[1..n2], this routine returns the K–S statistic d, and the significance level prob for the null hypothesis that the data sets aredrawn from the same distribution. Small values of prob show that the cumulative distributionfunction of data1 is significantly different from that of data2.

The arrays data1 and data2are modified by being sorted into ascending order.{float probks(float alam);void sort(unsigned long n, float arr[]);unsigned long j1=1,j2=1;float d1,d2,dt,en1,en2,en,fn1=0.0,fn2=0.0;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).for the case (14.3.6) of two distributions, where N1 is the number of data points inthe first distribution, N2 the number in the second.The nature of the approximation involved in (14.3.9) is that it becomesasymptotically accurate as the Ne becomes large, but is already quite good forNe ≥ 4, as small a number as one might ever actually use.

(See [1].)So, we have the following routines for the cases of one and two distributions:626Chapter 14.Statistical Description of DataIf we are not done...Next step is in data1.Next step is in data2.Compute significance.}Both of the above routines use the following routine for calculating the functionQKS :#include <math.h>#define EPS1 0.001#define EPS2 1.0e-8float probks(float alam)Kolmogorov-Smirnov probability function.{int j;float a2,fac=2.0,sum=0.0,term,termbf=0.0;a2 = -2.0*alam*alam;for (j=1;j<=100;j++) {term=fac*exp(a2*j*j);sum += term;if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum) return sum;fac = -fac;Alternating signs in sum.termbf=fabs(term);}return 1.0;Get here only by failing to converge.}Variants on the K–S TestThe sensitivity of the K–S test to deviations from a cumulative distribution functionP (x) is not independent of x.

In fact, the K–S test tends to be most sensitive around themedian value, where P (x) = 0.5, and less sensitive at the extreme ends of the distribution,where P (x) is near 0 or 1. The reason is that the difference |SN (x) − P (x)| does not, in thenull hypothesis, have a probability distribution that is independent of x. Rather, its variance isproportional to P (x)[1 − P (x)], which is largest at P = 0.5. Since the K–S statistic (14.3.5)is the maximum difference over all x of two cumulative distribution functions, a deviation thatmight be statistically significant at its own value of x gets compared to the expected chancedeviation at P = 0.5, and is thus discounted.

A result is that, while the K–S test is good atfinding shifts in a probability distribution, especially changes in the median value, it is notalways so good at finding spreads, which more affect the tails of the probability distribution,and which may leave the median unchanged.One way of increasing the power of the K–S statistic out on the tails is to replaceD (equation 14.3.5) by a so-called stabilized or weighted statistic [2-4], for example theAnderson-Darling statistic,D* =max−∞<x<∞|SN (x) − P (x)|pP (x)[1 − P (x)](14.3.11)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.