Phil Lucht Math & Physics Archive
Home / Math and Physics Files / Math / Scheid and numerical / Numerical Recipes in Fortran

f14-2

PDF · 6 pages · 62.9 KB
Open PDF file

Sample pages (about pp. 609-614) from Chapter 14 of Numerical Recipes in Fortran 77, a published textbook by others, filed in Phil's numerical-methods folder. It covers Student's t-test for equal variances, the unequal-variance t-test and the paired-sample t-test, with Fortran routines ttest, tutest, tptest and avevar. It also begins the F-test for variances and ends the preceding section on median and mode.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)
14.2DoTwoDistributionsHavetheSameMeansorVariances? 609Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If a distribution has a strong central tendency, so that most of its area is under a single peak, then the median is an estimator of the central value. It is a morerobust estimator than the mean is: The median fails as an estimator only if the area in the tails is large, while the mean fails if the first moment of the tails is large; it is easy to construct examples where the first moment of the tails is large eventhough their area is negligible. To find the median of a set of values, one can proceed by sorting the set and thenapplying(14.1.14). Thisisaprocessoforder NlogN. Youmightrightlythink that this is wasteful, since it yields much more information than just the median (e.g., the upper and lower quartile points, the deciles, etc.). In fact, we saw in§8.5 that the element x (N+1) /2can be located in of order Noperations. Consult that section for routines. Themodeof a probability distribution function p(x)is the value of xwhere it takesonamaximumvalue. Themodeisusefulprimarilywhenthereisasingle,sharp maximum, in which case it estimates the central value. Occasionally, a distribution will bebimodal, with two relative maxima; then one may wish to know the two modes individually. Note that, in such cases, the mean and median are not very useful, since they will give only a “compromise”value between the two peaks. CITED REFERENCES AND FURTHER READING: Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Chapter 2. Stuart, A., andOrd, J.K. 1987, Kendall’sAdvanced Theory of Statistics , 5th ed. (London:Griffin and Co.) [previous eds. published as Kendall, M., and Stuart, A., The Advanced Theory of Statistics ], vol. 1, §10.15 Norusis,M.J.1982, SPSSIntroductoryGuide:BasicStatisticsandOperations ;and1985, SPSS- X Advanced Statistics Guide (New York: McGraw-Hill). Chan,T.F.,Golub,G.H.,andLeVeque,R.J.1983, AmericanStatistician ,vol.37,pp.242–247.[1] Cram´er, H. 1946, Mathematical Methods of Statistics (Princeton: Princeton University Press), §15.10. [2] 14.2 Do Two Distributions Have the Same Means or Variances? Not uncommonly we want to know whether two distributions have the same mean. For example, a first set of measured values may have been gathered before some event,a secondset after it. We want to knowwhetherthe event,a “treatment” or a “change in a control parameter,” made a difference. Ourfirst thoughtis to ask“howmanystandarddeviations”onesamplemeanis from the other. That numbermay in fact be a useful thing to know. It does relate to the strength or “importance” of a difference of means if that difference is genuine . However, by itself, it says nothing about whether the difference isgenuine, that is, statistically significant. A difference of means can be very small compared to the standard deviation, and yet very significant, if the number of data points is large. Conversely, a difference may be moderately large but not significant, if the data 610 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).are sparse. We will be meeting these distinct concepts of strengthandsignificance several times in the next few sections. A quantity that measures the significance of a difference of means is not the number of standard deviations that they are apart, but the number of so-called standard errors that they are apart. The standard error of a set of values measures theaccuracywithwhichthesamplemeanestimatesthepopulation(or“true”)mean. Typically the standard error is equal to the sample’s standard deviation divided by the square root of the number of points in the sample. Student’s t-test for Significantly Different Means Applyingtheconceptofstandarderror,theconventionalstatistic formeasuring the significance of a difference of means is termed Student’s t . When the two distributions are thought to have the same variance, but possibly different means, then Student’s tis computed as follows: First, estimate the standard error of the difference of the means, sD, from the “pooled variance” by the formula sD=/radicalBigg/summationtext i∈A(xi−xA)2+/summationtext i∈B(xi−xB)2 NA+NB−2/parenleftbigg1 NA+1 NB/parenrightbigg (14.2.1 ) where each sum is over the points in one sample, the first or second, each mean likewisereferstoonesampleortheother,and NAandNBarethenumbersofpoints in the first and second samples, respectively. Second, compute tby t=xA−xB sD(14.2.2 ) Third, evaluate the significance of this value of tfor Student’s distribution with NA+NB−2degrees of freedom, by equations (6.4.7) and (6.4.9), and by the routine betai(incomplete beta function) of §6.4. The significance is a number between zero and one, and is the probability that |t|could be this large or larger just by chance, for distributions with equal means. Therefore,a small numerical value of the significance (0.05 or 0.01)means that the observed difference is “very significant.” The function A(t|ν)in equation (6.4.7) is one minus the significance. As a routine, we have SUBROUTINE ttest(data1,n1,data2,n2,t,prob) INTEGER n1,n2REAL prob,t,data1(n1),data2(n2) C USES avevar,betai Given the arrays data1(1:n1) anddata2(1:n2) , this routine returns Student’s tast, anditssignificanceas prob,smallvaluesof probindicatingthatthearrayshavesignificantly different means. The data arrays are assumed to be drawn from populations with the same true variance. REAL ave1,ave2,df,var,var1,var2,betaicall avevar(data1,n1,ave1,var1) call avevar(data2,n2,ave2,var2) df=n1+n2-2 Degrees of freedom. var=((n1-1)*var1+(n2-1)*var2)/df Pooled variance. t=(ave1-ave2)/sqrt(var*(1./n1+1./n2)) prob=betai(0.5*df,0.5,df/(df+t**2)) See equation (6.4.9). returnEND 14.2DoTwoDistributionsHavetheSameMeansorVariances? 611Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).which makes use of the following routine for computing the mean and variance of a set of numbers, SUBROUTINE avevar(data,n,ave,var) INTEGER nREAL ave,var,data(n) Given array data(1:n) , returns its mean as aveand its variance as var. INTEGER j REAL s,epave=0.0 do 11j=1,n ave=ave+data(j) enddo 11 ave=ave/n var=0.0 ep=0.0do 12j=1,n s=data(j)-ave ep=ep+svar=var+s*s enddo 12 var=(var-ep**2/n)/(n-1) Corrected two-pass formula (14.1.8). returnEND The next case to consider is where the two distributions have significantly differentvariances,but we neverthelesswant to knowif theirmeans arethe same or different. (A treatment for baldness has caused some patients to loseall their hair andturnedothersintowerewolves,but we want toknowif it helpscurebaldness on the average !) Be suspiciousoftheunequal-variance t-test: Iftwodistributionshave very different variances, then they may also be substantially different in shape; in thatcase,thedifferenceofthemeansmaynotbeaparticularlyusefulthingtoknow. To find out whether the two data sets have variances that are significantly different, you use the F-test, described later on in this section. The relevant statistic for the unequal variance t-test is t=xA−xB [Var (xA)/N A+Var (xB)/N B]1/2(14.2.3 ) This statistic is distributed approximately as Student’s twith a number of degrees of freedom equal to /bracketleftbigg Var (xA) NA+Var (xB) NB/bracketrightbigg2 [Var (xA)/N A]2 NA−1+[Var (xB)/N B]2 NB−1(14.2.4 ) Expression(14.2.4)is in generalnot aninteger,but equation(6.4.7)doesn’tcare. The routine is SUBROUTINE tutest(data1,n1,data2,n2,t,prob) INTEGER n1,n2 REAL prob,t,data1(n1),data2(n2) C USES avevar,betai Given the arrays data1(1:n1) anddata2(1:n2) , this routine returns Student’s tast, anditssignificanceas prob,smallvaluesof probindicatingthatthearrayshavesignificantly 612 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).different means. The data arrays are allowed to be drawn from populations with unequal variances. REAL ave1,ave2,df,var1,var2,betai call avevar(data1,n1,ave1,var1)call avevar(data2,n2,ave2,var2) t=(ave1-ave2)/sqrt(var1/n1+var2/n2) df=(var1/n1+var2/n2)**2/((var1/n1)**2/(n1-1)+(var2/n2)**2/(n2-1))prob=betai(0.5*df,0.5,df/(df+t**2))return END Our final example of a Student’s ttest is the case of paired samples . Here we imagine that much of the variance in bothsamples is due to effects that are point-by-point identical in the two samples. For example, we might have two jobcandidateswhohaveeachbeenratedbythesametenmembersofahiringcommittee. We want to know if the means of the ten scores differ significantly. We first try ttestabove, and obtain a value of probthat is not especially significant (e.g., >0.05). But perhaps the significance is being washed out by the tendency of some committee members always to give high scores, others always to give low scores, which increases the apparent variance and thus decreases the significance of any difference in the means. We thus try the paired-sample formulas, Cov (x A,x B)≡1 N−1N/summationdisplay i=1(xAi−xA)(xBi−xB)( 14.2.5 ) sD=/bracketleftbiggVar (xA)+Var (xB)−2Cov (xA,x B) N/bracketrightbigg1/2 (14.2.6 ) t=xA−xB sD(14.2.7 ) where Nis thenumberineachsample(numberofpairs). Noticethatit is important that a particular value of ilabel the corresponding points in each sample, that is, the ones that are paired. The significance of the tstatistic in (14.2.7) is evaluated forN−1degrees of freedom. The routine is SUBROUTINE tptest(data1,data2,n,t,prob) INTEGER n REAL prob,t,data1(n),data2(n) C USES avevar,betai Given the paired arrays data1(1:n) anddata2(1:n) , this routine returns Student’s tfor paired data as t, and its significance as prob, small values of probindicating a significant difference of means. INTEGER j REAL ave1,ave2,cov,df,sd,var1,var2,betai call avevar(data1,n,ave1,var1)call avevar(data2,n,ave2,var2) cov=0. do 11j=1,n cov=cov+(data1(j)-ave1)*(data2(j)-ave2) enddo 11 df=n-1cov=cov/dfsd=sqrt((var1+var2-2.*cov)/n) t=(ave1-ave2)/sd 14.2DoTwoDistributionsHavetheSameMeansorVariances? 613Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).prob=betai(0.5*df,0.5,df/(df+t**2)) return END F-Test for SignificantlyDifferent Variances TheF-testtests the hypothesis that two samples have different variances by trying to reject the null hypothesis that their variances are actually consistent. The statistic Fis the ratio of one variance to the other, so values either /greatermuch 1or/lessmuch 1 will indicate very significant differences. The distribution of Fin the null case is given in equation (6.4.11),which is evaluated using the routine betai. In the most commoncase, wearewillingto disprovethenullhypothesis(ofequalvariances)by either very large or verysmall values of F, so the correct significance is two-tailed , the sum of two incomplete beta functions. It turns out, by equation (6.4.3),that the two tails are always equal; we need compute only one, and doubleit. Occasionally,when the null hypothesisis stronglyviable, the identity of the two tails can become confused,givinganindicatedprobabilitygreaterthanone. Changingtheprobability to two minus itself correctlyexchangesthe tails. These considerationsand equation(6.4.3) give the routine SUBROUTINE ftest(data1,n1,data2,n2,f,prob) INTEGER n1,n2 REAL f,prob,data1(n1),data2(n2) C USES avevar,betai Given the arrays data1(1:n1) anddata2(1:n2) , this routine returns the value of f,a n d its significance as prob. Smallvalues of probindicate that the two arrays have significantly different variances. REAL ave1,ave2,df1,df2,var1,var2,betaicall avevar(data1,n1,ave1,var1) call avevar(data2,n2,ave2,var2) if(var1.gt.var2)then Make Fthe ratio of the larger variance to the smaller one. f=var1/var2df1=n1-1 df2=n2-1 else f=var2/var1 df1=n2-1 df2=n1-1 endifprob=2.*betai(0.5*df2,0.5*df1,df2/(df2+df1*f)) if(prob.gt.1.)prob=2.-prob returnEND CITED REFERENCES AND FURTHER READING: von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: Academic Press), Chapter IX(B). Norusis,M.J.1982, SPSSIntroductoryGuide:BasicStatisticsandOperations ;and1985, SPSS- X Advanced Statistics Guide (New York: McGraw-Hill). 614 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).14.3 Are Two Distributions Different? Given two sets of data, we can generalize the questions asked in the previous sectionandaskthesinglequestion: Arethetwosetsdrawnfromthesamedistribution function, or from different distribution functions? Equivalently,in proper statistical language, “Can we disprove, to a certain required level of significance, the null hypothesis that two data sets are drawn from the same population distribution function?” Disprovingthenullhypothesisineffectprovesthatthedatasetsarefromdifferent distributions. Failing to disprove the null hypothesis, on the other hand, only shows that the data sets can be consistent with a single distribution function. One can never provethat two data sets come from a single distribution, since (e.g.) no practical amount of data can distinguish between two distributions which differ only by one part in 10 10. Provingthattwodistributionsaredifferent,orshowingthattheyareconsistent, is a task that comes up all the time in many areas of research: Are the visible stars distributed uniformly in the sky? (That is, is the distribution of stars as a functionof declination — position in the sky — the same as the distribution of sky area as a function of declination?) Are educational patterns the same in Brooklyn as in the Bronx? (That is, are the distributions of people as a function of last-grade-attendedthe same?) Do two brands of fluorescent lights have the same distribution of burn-outtimes? Istheincidenceofchickenpoxthesameforfirst-born,second-born, third-born children, etc.? Thesefourexamplesillustratethefourcombinationsarisingfromtwodifferent dichotomies: (1) The data are either continuous or binned. (2) Either we wish tocompare one data set to a known distribution, or we wish to compare two equally unknown data sets. The data sets on fluorescent lights and on stars are continuous, since we can be given lists of individual burnout times or of stellar positions. Thedata sets on chicken pox and educational level are binned, since we are given tables of numbers of events in discrete categories: first-born, second-born, etc.; or 6th Grade, 7th Grade, etc. Stars and chicken pox, on the other hand, share the property that the null hypothesis is a known distribution (distribution of area in the sky, or incidence of chicken pox in the general population). Fluorescent lights andeducationallevel involvethe comparisonoftwo equallyunknowndatasets (the two brands, or Brooklyn and the Bronx). One can always turn continuous data into binned data, by grouping the events into specified ranges of the continuous variable(s): declinations between 0 and 10 degrees,10and20,20and30,etc. Binninginvolvesa loss ofinformation,however. Also, there is often considerable arbitrariness as to how the bins should be chosen. Alongwithmanyotherinvestigators,weprefertoavoidunnecessarybinningofdata. Theacceptedtestfordifferencesbetweenbinneddistributionsisthe chi-square test. For continuous data as a function of a single variable, the most generally accepted test is the Kolmogorov-Smirnovtest . We consider each in turn. Chi-Square Test Suppose that Niis the number of events observed in the ith bin, and that niis the number expected according to some known distribution. Note that the Ni’s are