PROGRAM ROSTAT c***************************************************************************** c c VERSION 1.2 February.1991 c c--- this routine contains the current versions of statistical c routines being tested by T. Beers, K. Flynn, and K. Gebhardt for robust c estimation of simple statistics c c***************************************************************************** c c-- CHANGES SINCE V1.0 (see update document for more information) c c-- 10/10/90 -- include the original estimate of the biweight when c error corrections are desired c c-- 10/10/90 -- change dip test so it doesn't loop forever c c-- 10/10/90 -- add check so KSREJ does not clip to N=5 c c-- 10/11/90 -- changed calculation of acceleration constant c c-- 1/ 1/91 -- changed b1 & b2 test and added probability calculation c c-- CHANGES SINCE V1.1 (see update document for more information) c c-- 1/28/91 -- included the term of 1/1.349 in the calculation of the c confidence interval for the median and fourth spread c c-- 1/28/91 -- the calculation for the confidence interval for the 3-sigma c mean and stand. dev. was changed in confidence and PRINTER c was updated c c-- 1/28/91 -- PRINTER was changed so it now writes out the proper T-VALUES c which are calculated in the program c c-- 1/28/91 -- Changed calculation of B1-test so it gives the left side c probability when the skewness is less than 0 c c-- 1/30/91 -- Fixed default number of decimal places for output c c-- 1/31/91 -- Added check so that SIGMA does not clip to N = 5 c c-- 2/ 1/91 -- Added check so that KURTSKEW is skipped when N < 9 c c-- 2/ 1/91 -- Changed calculation of B2-test so it gives the left side c probability when the kurtosis is less than 3 c c-- 2/ 1/91 -- The subroutine KSREJ now calls EDFGOF directly, rather than c going through the entire NORMALITY subroutine c c-- 2/ 1/91 -- Updated PRINTER so that exit diagnostics ( values = -999 ) c are properly output c c-- CHANGES SINCE V1.2 (see update document for more information) c c-- 4/18/91 -- Updated main program so that it will not call TAIL,KSREJ,TRANS c if n is less than 7, but so that it will still call SIGMA to c n = 5 (check is inside SIGMA) c c-- 4/18/91 -- Updated CONFIDENCE so that it will return intervals = 0 if c the estimate of XSDEV = -999, indicating that it was not c calculated earlier in the program c c-- 4/24/91 -- Corrected typo in data statement of WEXT. Apparently W test c is unaffected. c c-- 4/26/91 -- Added a check in the W probability calculation of WEXT to c guard against underflows in call to ALNORM c c-- 4/29/91 -- Changed KURTSKEW so it can handle larger datasets c without blowing up. c c-- 8/27/91 -- Updated subroutine GAPPER to calculate approximate "per gap" c probabilities as well as distributional probabilities for gaps c c***************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) integer iunit2,iunit3,iunit4,ounit3 real*4 mcave,mcsdev,mletter(8),mcconf(32,4),conf(4) dimension xdata(iii),odata(iii),xletter(15,2),depths(15) dimension xerr(iii),oerr(iii) dimension xmsc(4),xmdc(4),xtsc(4),xjkc(4),xjlc(4) dimension xxscl(4),xxscu(4),xjbc(4),xjblc(4),xjgc(4),xjglc(4) dimension clvalues(4),cuvalues(4),tvalues(12),test(26) dimension sletter(15),rletter(8),cletter(8) dimension r(iii*2) dimension zstar(iii),zbig(100),jbig(100),zp(100),tp(100) character*32 filename character*70 outputfile1 character DATADIR*(*) parameter (DATADIR='/export/home2/posts/rostat/') data zero,cee /0.0,2.99792458E+05/ data iunit2,iunit3,iunit4,ounit3 /2,3,4,9/ common iunit2,iunit3,iunit4,ounit3 common /print/ jtnor,jtgap,jtcon,jtchi,jtasy,jtcorr,jtot c--- open the data tables (added 17 april 92 by tina bird) c c substitute your local path here!!! c c********************************************************************** open(unit=2,file=DATADIR // 'ttable.dat',status='old') open(unit=3,file=DATADIR // 'cmin.dat',status='old') open(unit=4,file=DATADIR // 'cplu.dat',status='old') open(unit=8,file=DATADIR // 'corr.dat',status='old') c********************************************************************** c--- get the data print *, 'INPUT DATA FILE' read(5,*) filename open(unit=20,file=filename,status='old') c--- open the output files write(6,*) 'Outputfile for short printout?' read(5,*) outputfile1 open(unit=9,file=outputfile1,status='new') c--- typical inputs print *, 'Including Velocity Errors? 1-yes, 0-no' read *,ierr nin=0 print *,'What is number of simulations in Boot?' read *,nsims conf(1)=68 ! confidence interval percentages in monte conf(2)=90 conf(3)=95 conf(4)=99 do 331 i=1,iii if (ierr) 50,50,51 51 read(20,*,end=15) odata(i),oerr(i) goto 52 50 read(20,*,end=15) odata(i) 52 nin=nin+1 331 continue 15 xlbiwt = zero ! initially c----- If cosmological errors are requested we assume by default we c are using ROSTAT for clusters of galaxies, and approximate errors c are calculated taking the biweight estimates of location and scale, c then reanalyzing the data set after scaling by CEE if (ierr.eq.1) then ! adjust velocities j=0 do 392 i=1,nin if (odata(i).ne.-32368.) then j=j+1 xerr(j)=oerr(i) xdata(j)=odata(i) endif 392 continue CALL XBIWT (XDATA,J,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) vt2=xlbiwt c-- make the cosmological correction here do 393 i=1,j 393 xdata(i)=(xdata(i)-xlbiwt)/(1.+xlbiwt/cee) else do 394 i=1,nin if (odata(i).ne.-32368.) then j=j+1 xdata(j)=odata(i) endif 394 continue endif n=j 395 continue c--- here we initialize the estimators xave = zero ! straight mean xsdev = zero ! standard deviation xadev = zero ! average deviation (from mean) xskew = zero ! skewness xcurt = zero ! curtosis xvar = zero ! variance xdf = zero ! fourth spread xdfs = zero ! standardized fourth spread xsave = zero ! 3-sigma mean xssdev = zero ! 3-sigma standard deviation xkave = zero ! K-S mean xkdev = zero ! K-S standard deviation xtrim5 = zero ! 5% trimmed mean xtrim10 = zero ! 10% trimmed mean xtrim20 = zero ! 20% trimmed mean xtrim30 = zero ! 30% trimmed mean xtrim40 = zero ! 40% trimmed mean xmed = zero ! median xbmed = zero ! broadened median xtrimn = zero ! tri-mean xmid = zero ! mid-mean xmadm = zero ! median absolute deviation (fr median) xadm = zero ! mean absolute deviation (fr median) xlbiwt = zero ! full bi-weight location estimate xsbiwt = zero ! full bi-weight scale estimate xlbiwt1 = zero ! one-step biweight location estimate xsbiwt1 = zero ! one-step biweight scale estimate xljack1 = zero ! jacknife location of PARAM (1) (mean) xsjack1 = zero ! jacknife scale of PARAM (1) xlljack1 = zero ! jacknife (log) location of PARAM (1) xlsjack1 = zero ! jacknife (log) scale of PARAM (1) xljack2 = zero ! jacknife location of PARAM (2) (sd) xsjack2 = zero ! jacknife scale of PARAM (2) xlljack2 = zero ! jacknife (log) location of PARAM (2) xlsjack2 = zero ! jacknife (log) scale of PARAM (2) xljack6a = zero ! jacknife location of biweight location xsjack6a = zero ! jacknife scale of biwt location xlljack6a = zero ! jacknife (log) location of biwt loc xlsjack6a = zero ! jacknife (log) scale of biwt location xljack6 = zero ! jacknife location of biwt spread xsjack6 = zero ! jacknife scale of biwt spread xlljack6 = zero ! jacknife (log) location of biwt spread xlsjack6 = zero ! jacknife (log) scale of biwt spread xljack7 = zero ! jacknife location of PARAM (7) (gapper) xsjack7 = zero ! jacknife scale of PARAM (7) xlljack7 = zero ! jacknife (log) location of PARAM (7) xlsjack7 = zero ! jacknife (log) scale of PARAM (7) xsgap = zero ! gap estigate of scale tindex1 = zero ! estimate of tail weight tindex2 = zero ! tindex1 relative to a gaussian tindex pow = zero ! power transformation for symmetry tee = zero ! Finch's symmetry statistic c--- obtain some initial information about the data set c--- sort the data and find the set of letter values CALL LETTERS (XDATA,N,DEPTHS,XLETTER,MLETTER, + SLETTER,RLETTER,CLETTER,NLETTER) xmed = xletter(1,1) ! median xfl = xletter(2,1) ! lower f value xfu = xletter(2,2) ! upper f value xdf = abs(xfu-xfl) ! fourth spread xdfs = xdf/1.349 ! standardized fourth spread c--- perform calculations and report back CALL MOMENT (XDATA,N,XAVE,XADEV,XSDEV,XVAR,XSKEW,XCURT) CALL TRIM (XDATA,N,.05,XTRIM5) CALL TRIM (XDATA,N,.1,XTRIM10) CALL TRIM (XDATA,N,.2,XTRIM20) CALL XMIDMEAN (XDATA,N,XMID) CALL TRIM (XDATA,N,.3,XTRIM30) CALL TRIM (XDATA,N,.4,XTRIM40) CALL BMED (XDATA,N,XBMED) CALL TRIMEAN (XDATA,N,XMED,XFL,XFU,XTRIMN) CALL XMAD (XDATA,N,XMED,XMADM) CALL XAD (XDATA,N,XMED,XADM) CALL XBIWT (XDATA,N,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) CALL GAPPER (XDATA,N,XSGAP,ZSTAR,NBIG,JBIG,ZBIG,ZP,TP) if (ierr.eq.1) then ! apply corrections d1sum = 0 do 55 I=1,n 55 d1sum = d1sum+1./(xerr(i)*xerr(i)) davg=(n/d1sum)**0.5 C----- obtain corrections for clusters of galaxies cicorr=davg*davg ! add in quadrature to location interval scorr=davg*davg/(1.+vt2/cee)**2 ! subtract from scale estimate sicorr = scorr*(1.+scorr/(2.*xsbiwt*xsbiwt))/n ! add to scale endif CALL JACKNIFE (XDATA,N,1,XLJACK1,XSJACK1, ! mean and standard + XLLJACK1,XLSJACK1,XLJACK2, ! deviation + XSJACK2,XLLJACK2,XLSJACK2) CALL JACKNIFE (XDATA,N,6,XLJACK6a,XSJACK6a, ! biwt loc and spread + XLLJACK6a,XLSJACK6a,XLJACK6, ! 6a is location + XSJACK6,XLLJACK6,XLSJACK6) CALL JACKNIFE (XDATA,N,7,XLJACK7,XSJACK7, ! gapper + XLLJACK7,XLSJACK7,XJ71,XJ72, + XJ73,XJ74) if (n .ge. 7) then CALL TAIL (XDATA,N,TINDEX1,TINDEX2) CALL KSREJ (XDATA,N,XKAVE,XKDEV,NKTOT) CALL TRANS (N,NLETTER,XLETTER,POW) endif CALL SIGMA (XDATA,N,XSAVE,XSSDEV,NSTOT) ! will work for n >= 5 c-- If isin is set equal to 1 then the canonical confidence interval c used will be calculated using the 3-sigma mean and standard c deviation. If isin = 0, the mean and standard deviation will be used . isin=1 if (isin.eq.0) then nsin=n xsin=xsdev else if (isin.eq.1) then nsin=nstot xsin=xssdev endif 769 CALL CONFIDENCE (N,NSIN,XSIN,XDFS,XSBIWT,XSJACK2, + XLSJACK2,XSJACK6,XLSJACK6,XSJACK7,XLSJACK7, + XMSC,XMDC,XTSC,XJKC,XJLC,XJBC,XJBLC, + XJGC,XJGLC,XXSCL,XXSCU, + CLVALUES,CUVALUES,TVALUES) if(nsims.eq.0) goto 779 ! skip bootstrap CALL MCCON (XDATA,N,CONF,XAVE,XSDEV,XLBIWT,XSBIWT, + MCCONF,NSIMS) c--- for mcconf(i,j) : j=1,2,3,4 --> interval percentages c i= 1-8 are mean limits c i= 9-16 are biwt loc limits c i= 17-24 are the stand dev limits c i= 25-32 are biwt spread limits 779 nt = 26 ! number of normality test results CALL NORMALITY (XDATA,N,XAVE,XSDEV,XSKEW,XCURT, + XSBIWT,XMED,TEST) CALL SYMMETRY (XDATA,N,TEE) c--- fill up a array of results to be saved and printed r(1) = n ! total number of points j=1 do 100 i=1,n ! sorted data (2 -> n+1) j=j+1 r(j) = xdata(i) 100 continue j=j+1 r(j) = nletter ! number of letter values (n+2) do 101 i=1,nletter ! letter values j=j+1 ! (n+3 -> n+2+7*nl) r(j) = depths(i) ! depths j=j+1 r(j)=xletter(i,1) ! lower j=j+1 r(j)=xletter(i,2) ! upper j=j+1 r(j) =mletter(i) ! midspread j=j+1 r(j)=sletter(i) ! letter spread j=j+1 r(j) = rletter(i) ! gaussian ratio spread j=j+1 r(j) = cletter(i) ! corrected pseudosigmas 101 continue j=j+1 ! classical estimators (n+7*nl+...) r(j) = xave ! straight mean (3) j=j+1 r(j) = xadev ! average deviation from mean (4) j=j+1 r(j) = xsdev ! standard deviation from mean (5) j=j+1 r(j) = xvar ! variance (6) j=j+1 r(j) = xskew ! skewness (7) j=j+1 r(j) = xcurt ! curtosis (8) j=j+1 ! location estimates (n+7*nl+...) r(j) = xave ! straight mean (9) j=j+1 r(j) = nstot ! data points for xsave (10) j=j+1 r(j) = xsave ! three sigma mean (11) j=j+1 r(j) = nktot ! data points for K-S (12) j=j+1 r(j) = xkave ! K-S mean (13) j=j+1 r(j) = xtrim5 ! 5% trimmed mean (14) j=j+1 r(j) = xtrim10 ! 10% trimmed mean (15) j=j+1 r(j) = xtrim20 ! 20% trimmed mean (16) j=j+1 r(j) = xmid ! midmean (17) j=j+1 r(j) = xtrim30 ! 30% trimmed mean (18) j=j+1 r(j) = xtrim40 ! 40% trimmed mean (19) j=j+1 r(j) = xmed ! median (20) j=j+1 r(j) = xbmed ! broadened median (21) j=j+1 r(j) = xtrimn ! trimean (22) j=j+1 r(j) = xlbiwt ! full biweight location (23) j=j+1 r(j) = XLBIWT1 ! one step biweight location (24) j=j+1 r(j) = xljack1 ! jacknife location of PARAM (1) (25) j=j+1 r(j) = xsjack1 ! jacknife scale of PARAM (1) (26) j=j+1 r(j) = xlljack1 ! jacknife log locat of PARAM (1) (27) j=j+1 r(j) = xlsjack1 ! jacknife log scale of PARAM (1) (28) j=j+1 r(j) = xlboot ! bootstrap location of PARAM (29) j=j+1 ! scale estimators (n+7*nl+...) r(j) = xsdev ! standard deviation (30) j=j+1 r(j) = nstot ! number for xsdev (31) j=j+1 r(j) = xssdev ! three sigma standard deviation (32) j=j+1 r(j) = nktot ! number for ksdev (33) j=j+1 r(j) = xkdev ! K-S standard deviation (34) j=j+1 r(j) = xadev ! average deviation from mean (35) j=j+1 r(j) = xadm ! mean abs deviation (fr median) (36) j=j+1 r(j) = xdf ! fourth spread (37) j=j+1 r(j) = xdf/1.349 ! standardized fourth spread (38) j=j+1 r(j) = xmadm ! median abs deviation (fr median) (39) j=j+1 r(j) = xmadm/0.6745 ! standardized MAD (40) j=j+1 r(j) = xsbiwt ! full biweight scale (41) j=j+1 r(j) = xsbiwt1 ! one-step biweight scale (42) j=j+1 r(j) = xljack2 ! jacknife locat of PARAM (2) (43) j=j+1 r(j) = xsjack2 ! jacknife scale of PARAM (2) (44) j=j+1 r(j) = xlljack2 ! jacknife log locat of PARAM (2) (45) j=j+1 r(j) = xlsjack2 ! jacknife log scale of PARAM (2) (46) j = j+1 r(j) = xljack6 ! jacknife locat of param (6) (47) j=j+1 r(j) = xsjack6 ! jacknife scale of param (6) (48) j=j+1 r(j) = xlljack6 ! jacknife log locat of param (6) (49) j = j+1 r(j) = xlsjack6 ! jacknife log scale of param (6) (50) j = j+1 r(j) = xljack7 ! jacknife locat of param (7) (51) j=j+1 r(j) = xsjack7 ! jacknife scale of param (7) (52) j=j+1 r(j) = xlljack7 ! jacknife log locat of param (7) (53) j = j+1 r(j) = xlsjack7 ! jacknife log scale of param (7) (54) j = j+1 r(j) = xsboot ! bootstrap scale of PARAM (55) j=j+1 ! tail index 1 (raw) (56) r(j) = tindex1 j=j+1 r(j) = tindex2 ! tail index 2 (scaled) (57) j=j+1 ! gap estimate of scale (58) r(j) = xsgap j=j+1 r(j) = vt2 ! original est. of Biweight loc.(59) j=j+1 r(j) = ierr ! including errors(1) or not(0) (60) jtnor = n+7*nletter+60 do 103 i=1,nt ! normality tests (jtnor+1->jtnor+nt) j=j+1 ! r(j) = test(i) 103 continue ! a-test statistic (jtnor+1) ! u-test statistic (2) ! w-test statistic (3) ! w-test probability (4) ! b1-test statistic (5) ! b2-test statistic (6) ! i-test statistic (7) ! i-test 90% prob (8) ! i-test 95% prob (9) ! dip-test statistic (10) ! dip-test lower modal (11) ! dip-test upper modal (12) ! KS-test statistic (13) ! KS-test probability (14) ! Kuiper V-test statistic (15) ! Kuiper V-test probability (16) ! Cramer Von-Mises W2-test stat (17) ! Cramer Von-Mises W2-test prob (18) ! Watson U2-test statistic (19) ! Watson U2-test statistic (20) ! Anderson-Darling A2-test stat (21) ! Anderson-Darling A2-test prob (22) ! b1 test probability (23) ! b2 test probability (24) ! ominbus test for b1 & b2 (25) ! chi-square prob for omni test (26) jtgap = jtnor+nt ngap = n-1 ! gap test results j=j+1 r(j) = ngap ! number of gaps (jtgap+1) do 112 i=1,ngap ! xsgap,zstar (jtgap+2->jtgap+1+ngap) j=j+1 r(j)= zstar(i) 112 continue j=j+1 r(j) = nbig ! number of big gaps (2+jtgap+ngap) do 113 i=1,nbig ! index and big gaps (3+jtgap+ngap-> j=j+1 ! 2+jtgap+ngap+4*nbig) r(j) = jbig(i) j=j+1 r(j) = zbig(i) j=j+1 r(j) = zp(i) ! gap probability j=j+1 r(j) = tp(i) ! distribution probability 113 continue jtcon = 2+jtgap+ngap+4*nbig do 104 i= 1,4 ! confidence intervals j=j+1 r(j) = xmsc(i) ! mean and standard deviation 104 continue ! (jtcon+1 -> jtcon+4) do 105,i=1,4 j=j+1 r(j) = xmdc(i) ! median and fourth spread 105 continue ! (jtcon+5 -> jtcon+8) do 106 i=1,4 j=j+1 r(j) = xtsc(i) ! biweight location and spread 106 continue ! (jtcon+9 -> jtcon+12) do 107 i=1,4 j=j+1 r(j) = xjkc(i) ! jacknife estimates 107 continue ! (jtcon+13 -> jtcon+16) do 108 i=1,4 j=j+1 r(j) = xjlc(i) ! log jacknife estimates 108 continue ! (jtcon+17 -> jtcon+20) do 121 i = 1,4 j = j+1 ! jacknife biwt r(j) = xjbc(i) ! (jtcon+21 -> jtcon+24) 121 continue do 122 i = 1,4 j=j+1 ! log jacknife biwt r(j) = xjblc(i) ! (jtcon+25 -> jtcon+28) 122 continue do 123 i = 1,4 j=j+1 ! jacknife gapped r(j) = xjgc(i) ! (jtcon+29 -> jtcon+32) 123 continue do 124 i = 1,4 j=j+1 ! log jacknife gapped r(j) = xjglc(i) ! (jtcon+33 -> jtcon+36) 124 continue do 125 i = 1,4 j=j+1 ! boot version 1 of the mean r(j) = mcconf(1,i) ! (jtcon+37 -> jtcon+44) j=j+1 r(j) = mcconf(2,i) 125 continue do 126 i = 1,4 j=j+1 ! boot version 2 of the mean r(j) = mcconf(3,i) ! (jtcon+45 -> jtcon+52) j=j+1 r(j) = mcconf(4,i) 126 continue do 127 i = 1,4 j=j+1 ! boot version 3 of the mean r(j) = mcconf(5,i) ! (jtcon+53 -> jtcon+60) j=j+1 r(j) = mcconf(6,i) 127 continue do 128 i = 1,4 j=j+1 ! boot version 1 of the biweight r(j) = mcconf(9,i) ! (jtcon+61 -> jtcon+68) j=j+1 r(j) = mcconf(10,i) 128 continue do 129 i = 1,4 j=j+1 ! boot version 2 of the biweight r(j) = mcconf(11,i) ! (jtcon+69 -> jtcon+76) j=j+1 r(j) = mcconf(12,i) 129 continue do 130 i = 1,4 j=j+1 ! boot version 3 of the biweight r(j) = mcconf(13,i) ! (jtcon+77 -> jtcon+84) j=j+1 r(j) = mcconf(14,i) 130 continue do 131 i = 1,4 j=j+1 ! boot version 1 of the stand dev r(j) = mcconf(17,i) ! (jtcon+85 -> jtcon+92) j=j+1 r(j) = mcconf(18,i) 131 continue do 132 i = 1,4 j=j+1 ! boot version 2 of the stand dev r(j) = mcconf(19,i) ! (jtcon+93 -> jtcon+100) j=j+1 r(j) = mcconf(20,i) 132 continue do 133 i = 1,4 j=j+1 ! boot version 3 of the stand dev r(j) = mcconf(21,i) ! (jtcon+101 -> jtcon+108) j=j+1 r(j) = mcconf(22,i) 133 continue do 134 i = 1,4 j=j+1 ! boot version 1 of the biwt scale r(j) = mcconf(25,i) ! (jtcon+109 -> jtcon+116) j=j+1 r(j) = mcconf(26,i) 134 continue do 135 i = 1,4 j=j+1 ! boot version 2 of the biwt scale r(j) = mcconf(27,i) ! (jtcon+117 -> jtcon+124) j=j+1 r(j) = mcconf(28,i) 135 continue do 136 i = 1,4 j=j+1 ! boot version 3 of the biwt scale r(j) = mcconf(29,i) ! (jtcon+125 -> jtcon+132) j=j+1 r(j) = mcconf(30,i) 136 continue do 137 i = 1,4 j=j+1 ! boot version 4 of the mean r(j) = mcconf(7,i) ! (jtcon+133 -> jtcon+140) j=j+1 r(j) = mcconf(8,i) 137 continue do 138 i = 1,4 j=j+1 ! boot version 4 of the biwt r(j) = mcconf(15,i) ! (jtcon+141 -> jtcon+148) j=j+1 r(j) = mcconf(16,i) 138 continue do 139 i = 1,4 j=j+1 ! boot version 4 of the stand dev r(j) = mcconf(23,i) ! (jtcon+149 -> jtcon+156) j=j+1 r(j) = mcconf(24,i) 139 continue do 140 i = 1,4 j=j+1 ! boot version 4 of the biwt scale r(j) = mcconf(31,i) ! (jtcon+157 -> jtcon+164) j=j+1 r(j) = mcconf(32,i) 140 continue j = j+1 r(j)=isin ! switch for 3-sigma clipped or jtchi = jtcon+165 ! usual mean,deviation input into ! xxsc confidence intervals ! (jtcon+165) do 109 i=1,4 j=j+1 r(j) = xxscl(i) ! stand dev chi square (low) j=j+1 ! (jtchi+1 -> jtchi+8) r(j) = xxscu(i) ! stand dev chi square (hi) 109 continue do 110 i=1,4 ! chi-square values j=j+1 r(j) = clvalues(i) ! lower (68,90,95,99) j=j+1 ! (jtchi+9 -> jtchi+16) r(j) = cuvalues(i) ! upper (68,90,95,99) 110 continue do 111 i=1,12 ! t values j=j+1 ! t(n-1) (jtchi+17 -> jtchi+20) r(j) = tvalues(i) ! t(n-1)/1.075 (jtchi+21 ->jtchi+24) 111 continue ! t(0.7(n-1)) (jtchi+25 -> jtchi+28) ! (68,90,95,99) jtasy = jtchi + 28 j=j+1 r(j) = tee ! asymmetry index (jtasy+1) j=j+1 r(j) = pow ! power law transformation (jtasy+2) jtcorr = jtasy + 2 j=j+1 r(j) = cicorr ! correction to location confidence j=j+1 ! (jtcorr+1) r(j)= scorr ! correction to scale estimate j=j+1 ! (jtcorr+2) r(j) = sicorr ! correction to scale confidence ! (jtcorr+3) jtot = j c--- print out the result array in a readible format CALL PRINTER (OUNIT3,FILENAME,R) 999 stop end c------------------------------------------------------------------------------ SUBROUTINE BIASAD (EST,NSIMS,ACC,SIMVAL,Z,LIMIT1,LIMIT2, + LIMITA1,LIMITA2,IFAULT) c------------------------------------------------------------------------------ c--- Finds adjustment required to implement bias-corrected percentile method c c uses functions ALNORM and PPND Algorithms AS 66 and AS 111 c dimension simval(nsims) data half,two/0.5E0,2.0E0/ j = 0 k = 0 c c--- m = (number of values less than est) + (half the c number equal to est), rounded up if not integral c do 30 m = 1,nsims if (simval(m) - est) 10,20,30 10 j = j+1 20 k = k+1 30 continue m = (j + k + 1)/2 if (m .eq. 0) ifault = 3 if (m .eq. nsims) ifault = 4 if (ifault .gt. 0) return zed = two*ppnd(float(m)/float(nsims),ifault) c c--- ifail cannot exceed 0 simce m .ge. 1 and m .le. nsims - 1 c fn1 = float(nsims + 1) limit1 = fn1*alnorm(zed - z,.false.) + half limit2 = fn1*alnorm(zed + z,.false.) + half c c-- calculate limits from the acceleration constant c xzed1 = zed/two + ((zed/two - z)/(1-acc*(zed/two-z))) xzed2 = zed/two + ((zed/two + z)/(1-acc*(zed/two+z))) limita1 = fn1*alnorm(xzed1,.false.) + half limita2 = fn1*alnorm(xzed2,.false.) + half return end c------------------------------------------------------------------------------ SUBROUTINE BMED (XDATA,N,XBMED) c------------------------------------------------------------------------------ c--- This subroutine calculates the broadened median for a set of N c ORDERED statistics in XDATA. The value of the broadened median c is returned as XBMED. XBMED is defined differently for different c values of N, and is not defined for N<5. XBMED is calculated by c taking a weighted average of the inner depths as specified c page 313 of UREDA. The weights are as follows: c c 5 <= n <= 12 n odd --> middle 3 ordered stats are averaged c n even --> middle 4 ordered stats are averaged c with weights 1/6 1/3 1/3 1/6 c c n > 12 n odd --> middle 5 ordered stats are averaged c n even --> middle 6 ordered stats are averaged c with weights 1/10,1/5,1/5,1/5,1/5,1/10 c c The broadened median can also be easily found by taking an alpha c trimmed mean with alpha defined as: c c 5 <= n <= 12 alpha = (.5 - 1.5/n) c c n > 12 alpha = (.5 - 2.5/n) c**************************************************************************** implicit real*4 (a-h,o-z) dimension xdata(n) data n2,d6,d5,d3,d10,d1 /2,6.0,5.0,3.0,10.0,1.0/ data n1,n3 /1,3/ do3 = d1/d3 do5 = d1/d5 do6 = d1/d6 do10 = d1/d10 c--- bmed is not defined for n<5 if (n .le. 12) then if (float(n)/float(n2) - int(n/n2) .eq. 0) then i1 = n/n2 - n1 i2 = n/n2 i3 = n/n2 + n1 i4 = n/n2 + n2 xbmed = do6*(xdata(i1)+xdata(i4))+do3*(xdata(i2) + +xdata(i3)) else i1 = int(n/n2) i2 = int(n/n2) + n1 i3 = int(n/n2) + n2 xbmed = (xdata(i1) + xdata(i2) + xdata(i3))*do3 endif else if (float(n)/float(n2) - int(n/n2) .eq. 0) then i1 = n/n2 - n2 i2 = n/n2 - n1 i3 = n/n2 i4 = n/n2 + n1 i5 = n/n2 + n2 i6 = n/n2 + n3 xbmed = do10*(xdata(i1)+xdata(i6)) + do5*(xdata(i2) + +xdata(i3)+xdata(i4)+xdata(i5)) else i1 = int(n/n2) - n1 i2 = int(n/n2) i3 = int(n/n2) + n1 i4 = int(n/n2) + n2 i5 = int(n/n2) + n3 xbmed = do5*(xdata(i1)+xdata(i2)+xdata(i3)+xdata(i4) + +xdata(i5)) endif endif return end c------------------------------------------------------------------------------ SUBROUTINE BOOTSTRAP (X,N,IP,NSIMS,X1SIM,X2SIM) c------------------------------------------------------------------------------ c--- This routine returns a BOOTSTRAP array (X1SIM,X2SIM) of the c simulated values. The initial seed is set by a call to the system clock. c c****************************************************************************** parameter (iii = 10000) real*4 x(n),xdata(iii),x1sim(nsims),x2sim(nsims) integer iran1(iii) c--- fill up an array of random numbers of length n (between 1 and n) c xx=SECNDS(0.0) ! call to system clock c k1=2*jifix(10.*xx)+1 c k2=2*jifix(100.*xx)+1 c do 6 i=1,100 c xx=RAN(K1) c6 xx=RAN(K2) idum = -1 do 20 j=1,nsims do 5 i=1,n iran1(i) = 1 + int(n*RAN1(IDUM)) 5 continue do 10 i=1,n xdata(i) = x(iran1(i)) 10 continue if (ip.eq.1 .or. ip.eq.2)then x1sim(j) = PARAM(XDATA,N,1) x2sim(j) = PARAM(XDATA,N,2) else if (ip.eq.6) then CALL XBIWT(XDATA,N,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) x1sim(j) = xlbiwt x2sim(j) = xsbiwt endif 20 continue return end c------------------------------------------------------------------------------ SUBROUTINE CONFIDENCE (N,NSTOT,XSDEV,XDF,XSBIWT,XSJACK2, + XLSJACK2,XSJACK6,XLSJACK6,XSJACK7,XLSJACK7, + XMSC,XMDC,XTSC,XJKC,XJLC,XJBC,XJBLC, + XJGC,XJGLC,XXSCL,XXSCU, + CLVALUES,CUVALUES,TVALUES) c------------------------------------------------------------------------------ c--- CONFIDENCE calculates the confidence intervals for a set of location c parameters, in this case the straight mean XAVE, the median XMED, c and the bi-weight estimator XLBIWT. The confidence intervals are c returned in the variables: c LOCATIONS and their CONFIDENCE INTERVALS c XMSC: uses mean and standard deviation c XMDC: uses median and fourth spread c XTSC: uses bi-weight location and spread c c c SCALES and their CONFIDENCE INTERVALS c c XJKC: uses jacknifed estimates c XJLC: uses log - jacknifed estimates c XJBC: uses jacknifed biwt estimates c XJBLC: uses log - jacknifed estimates c XJGC: uses jacknifed gapped estimates c XJGLC: uses log -jacknifed gapped estimates c XXSC: uses standard deviation (chi-square) c (returns lower and upper values (XXSCL,XXSCU)) c c NOTE: the above variables are arrays with return the c 68% (1), 90% (2), 95% (3), and 99% (4) c confidence values c c TVALUES: tvalues of the confidence estimator c c (1) 0.68t for XMSC, XJKC, and XJLC c (2) 0.90t for XMSC, XJKC, and XJLC c (3) 0.95t for XMSC, XJKC, and XJLC c (4) 0.99t for XMSC, XJKC, and XJLC c c (5) 0.68t for XMDC c (6) 0.90t for XMDC c (7) 0.95t for XMDC c (8) 0.99t for XMDC c c (9) 0.68t for XTSC c (10) 0.90t for XTSC c (11) 0.95t for XTSC c (12) 0.99t for XTSC c c c CVALUES: cvalues of the confidence estimator c c (1) 0.68c for XXSC (u and l) c (2) 0.90c fpr XXSC c (3) 0.95c for XXSC c (4) 0.99c for XXSC c c c--- NOTE that the appropriate t-values are also calculated and returned c these values are also useful for standard hypothesis testing c c c**************************************************************************** implicit real*4 (a-h,o-z) integer n dimension tvalues(12),clvalues(4),cuvalues(4) dimension xmsc(4),xmdc(4),xtsc(4),xjkc(4),xjlc(4), + xxscl(4),xxscu(4),xjbc(4),xjblc(4),xjgc(4),xjglc(4) c--- obtain the appropriate t-values and c-values for each test ndof = n-1 nsdof = nstot-1 tprob = .68 t68 = TVALUE (NDOF,TPROB) ts68 = TVALUE (NSDOF,TPROB) c68l = CVALUE (0,NSDOF,TPROB) c68u = CVALUE (1,NSDOF,TPROB) tprob = .90 t90 = TVALUE (NDOF,TPROB) ts90 = TVALUE (NSDOF,TPROB) c90l = CVALUE (0,NSDOF,TPROB) c90u = CVALUE (1,NSDOF,TPROB) tprob = .95 t95 = TVALUE (NDOF,TPROB) ts95 = TVALUE (NSDOF,TPROB) c95l = CVALUE (0,NSDOF,TPROB) c95u = CVALUE (1,NSDOF,TPROB) tprob = .99 t99 = TVALUE (NDOF,TPROB) ts99 = TVALUE (NSDOF,TPROB) c99l = CVALUE (0,NSDOF,TPROB) c99u = CVALUE (1,NSDOF,TPROB) c--- XMSC, XJKC, and XJLC use the straight percentage points of the c t-distribution with N-1 degrees of freedom tvalues(1) = t68 tvalues(2) = t90 tvalues(3) = t95 tvalues(4) = t99 c--- XMDC uses the approximation t* = t(n-1 DOF)/1.075 tvalues(5) = t68/1.075 tvalues(6) = t90/1.075 tvalues(7) = t95/1.075 tvalues(8) = t99/1.075 c--- XTSC uses the approximation t* = t with int(0.7(n-1)) DOF idof = int (0.7*ndof) tprob = .68 tvalues(9) = TVALUE (IDOF,TPROB) tprob = .90 tvalues(10) = TVALUE (IDOF,TPROB) tprob = .95 tvalues(11) = TVALUE (IDOF,TPROB) tprob = .99 tvalues(12) = TVALUE (IDOF,TPROB) c--- XXSC uses the percentage points of the chi-square distribution clvalues(1) = c68l clvalues(2) = c90l clvalues(3) = c95l clvalues(4) = c99l cuvalues(1) = c68u cuvalues(2) = c90u cuvalues(3) = c95u cuvalues(4) = c99u c--- obtain interval estimators for central locations if(xsdev .ne. -999) then xmsc(1) = ts68*xsdev/nstot**0.5 xmsc(2) = ts90*xsdev/nstot**0.5 xmsc(3) = ts95*xsdev/nstot**0.5 xmsc(4) = ts99*xsdev/nstot**0.5 endif xmdc(1) = tvalues(5)*xdf/n**0.5 xmdc(2) = tvalues(6)*xdf/n**0.5 xmdc(3) = tvalues(7)*xdf/n**0.5 xmdc(4) = tvalues(8)*xdf/n**0.5 xtsc(1) = tvalues(9)*xsbiwt/n**0.5 xtsc(2) = tvalues(10)*xsbiwt/n**0.5 xtsc(3) = tvalues(11)*xsbiwt/n**0.5 xtsc(4) = tvalues(12)*xsbiwt/n**0.5 c--- obtain interval estimators for scale xjkc(1) = tvalues(1)*xsjack2 xjkc(2) = tvalues(2)*xsjack2 xjkc(3) = tvalues(3)*xsjack2 xjkc(4) = tvalues(4)*xsjack2 xjlc(1) = tvalues(1)*xlsjack2 xjlc(2) = tvalues(2)*xlsjack2 xjlc(3) = tvalues(3)*xlsjack2 xjlc(4) = tvalues(4)*xlsjack2 xjbc(1) = tvalues(1)*xsjack6 xjbc(2) = tvalues(2)*xsjack6 xjbc(3) = tvalues(3)*xsjack6 xjbc(4) = tvalues(4)*xsjack6 xjblc(1) = tvalues(1)*xlsjack6 xjblc(2) = tvalues(2)*xlsjack6 xjblc(3) = tvalues(3)*xlsjack6 xjblc(4) = tvalues(4)*xlsjack6 xjgc(1) = tvalues(1)*xsjack7 xjgc(2) = tvalues(2)*xsjack7 xjgc(3) = tvalues(3)*xsjack7 xjgc(4) = tvalues(4)*xsjack7 xjglc(1) = tvalues(1)*xlsjack7 xjglc(2) = tvalues(2)*xlsjack7 xjglc(3) = tvalues(3)*xlsjack7 xjglc(4) = tvalues(4)*xlsjack7 if(xsdev .ne. -999) then xxscu(1) = sqrt(xsdev**2/(clvalues(1)/float(nsdof))) xxscu(2) = sqrt(xsdev**2/(clvalues(2)/float(nsdof))) xxscu(3) = sqrt(xsdev**2/(clvalues(3)/float(nsdof))) xxscu(4) = sqrt(xsdev**2/(clvalues(4)/float(nsdof))) xxscl(1) = sqrt(xsdev**2/(cuvalues(1)/float(nsdof))) xxscl(2) = sqrt(xsdev**2/(cuvalues(2)/float(nsdof))) xxscl(3) = sqrt(xsdev**2/(cuvalues(3)/float(nsdof))) xxscl(4) = sqrt(xsdev**2/(cuvalues(4)/float(nsdof))) endif return end c------------------------------------------------------------------------------ SUBROUTINE DIPTST (X,N,DIP,XL,XU,IFAULT,GCM,LCM,MN,MJ) c------------------------------------------------------------------------------ c--- Algorithm AS 217 Appl. Stat. (1985) Vol.34, No.3 c Hartigan, P.M. c c--- This routine provides an estimator which guages the probability that c the data is drawn from a parent population having a single mode. As c such, it can be used as a conservative test for multi-modality of c a distribution. c--- Does the dip calculation for an ordered vector X c using the greatest convex minorant and the least concave c majorant, skepping through the data using the change c points of these distributions. It returns the dip c statistic 'DIP' and the modal interval (XL,XU) c c--- NOTE that the array X is expected to already be sorted low c to high on input c c****************************************************************************** dimension x(n), mn(n),mj(n),lcm(n) integer gcm(n), high c c--- check that n is positive c ifault = 1 if (n .le. 0) return ifault = 0 c c--- check if n is one c if (n .eq. 1) goto 4 c c--- check that x is sorted c ifault = 2 do 3 k = 2,n if (x(k) .lt. x(k-1)) return 3 continue ifault = 0 c c--- check for all values of x identical and for 1 .lt. n .lt. 4 c if (x(n) .gt. x(1) .and. n .ge. 4) goto 5 4 xl = x(1) xu = x(n) dip = 0.0 return c c--- low contains the index of the current estimate of the c lower end of the modal interval, high contains the c index for the current upper end. c 5 fn = float(n) low = 1 high = n dip = 1.0/fn xl = x(low) xu = x(high) c c--- establish the indices over which combination is c necessary for the convex minorant fit c mn(1) = 1 do 28 j=2,n mn(j) = j-1 25 mnj = mn(j) mnmnj = mn(mnj) a = float(mnj - mnmnj) b = float(j - mnj) if (mnj .eq. 1 .or. (x(j) - x(mnj))*a .lt. (x(mnj)-x(mnmnj)) + *b) goto 28 mn(j) = mnmnj goto 25 28 continue c c--- establish the indices over which combination is c necessary for the concave majorant fit. c mj(n) = n na = n-1 do 34 jk = 1,na k = n - jk mj(k) = k + 1 32 mjk = mj(k) mjmjk = mj(mjk) a = float(mjk - mjmjk) b = float(k - mjk) if (mjk .eq. n .or. (x(k) - x(mjk))*a .lt. (x(mjk)-x(mjmjk)) + *b) goto 34 mj(k) = mjmjk goto 32 34 continue c c--- start the cycling c collect the change points for the gcm from high to low c icount=0 40 ic = 1 gcm(1) = high 42 igcm1 = gcm(ic) ic = ic + 1 gcm(ic) = mn(igcm1) if (gcm(ic) .gt. low) goto 42 icx = ic c c--- collect the change points for the lcm from low to high c ic = 1 lcm(1) = low 44 lcm1 = lcm(ic) ic = ic + 1 lcm(ic) = mj(lcm1) if(lcm(ic) .lt. high) goto 44 icv = ic c c--- icx,ix,ig are counters for the convex minorant c icv,iv,ih are counters for the concave majorant c ig = icx ih = icv c c--- find the largest distance greater the 'DIP' c between the gcm and the lcm from low th high c ix = icx - 1 iv = 2 d = 0.0 if(icx .ne. 2 .or. icv .ne. 2) goto 50 d=1.0/fn goto 65 50 igcmx = gcm(ix) lcmiv = lcm(iv) if (igcmx .gt. lcmiv) goto 55 c c--- if the next point of either the gcm or lcm is c from the lcm then calculate distance here c lcmiv1 = lcm(iv - 1) a = float(lcmiv - lcmiv1) b = float(igcmx - lcmiv1 - 1) dx = (x(igcmx) - x(lcmiv1)*a)/(fn*(x(lcmiv) - x(lcmiv1))) + - b/fn ix = ix - 1 if (dx .lt. d) goto 60 d = dx ig = ix +1 ih = iv goto 60 c c--- if the next point of either the gcm or lcm is c from the gcm then calculate distance here c 55 lcmiv = lcm(iv) igcm = gcm(ix) igcm1 = gcm(ix + 1) a = float(lcmiv - igcm1 + 1) b = float(igcm - igcm1) dx = a/fn - ((x(lcmiv) - x(igcm1))*b)/(fn*(x(igcm) + - x(igcm1))) iv = iv + 1 if (dx .lt. d) goto 60 d = dx ig = ix + 1 ih = iv - 1 60 if (ix .lt. 1) ix = 1 if (iv .gt. icv) iv = icv if (gcm(ix) .ne. lcm(iv)) goto 50 65 if (d .lt. dip) goto 100 c c--- caclulate the dips for the current low and high c c the dip for the convex minorant c dl = 0.0 if (ig .eq. icx) goto 80 icxa = icx - 1 do 76 j = ig,icxa temp = 1.0/ fn jb = gcm(j+1) je = gcm(j) if (je - jb .le. 1) goto 74 if (x(je) .eq. x(jb)) goto 74 a = float(je - jb) const = a/(fn*(x(je) - x(jb))) do 72 jr = jb,je b = float(jr - jb +1) t = b/ fn - (x(jr) - x(jb)) * const if (t .gt. temp) temp = t 72 continue 74 if (dl .lt. temp) dl = temp 76 continue c c--- the dip for the concave majorant c 80 du = 0.0 if (ih .eq. icv) goto 90 icva = icv - 1 do 88 k = ih,icva temp = 1.0/fn kb = lcm(k) ke = lcm(k + 1) if (ke - kb .le. 1) goto 86 if (x(ke) .eq. x(kb)) goto 86 a = float(ke - kb) const = a/ (fn * (x(ke) - x(kb))) do 84 kr = kb,ke b = float(kr - kb - 1) t = (x(kr) - x(kb)) * const - b/fn if (t .gt. temp) temp = t 84 continue 86 if (du .lt. temp) du = temp 88 continue c c--- determine the current maximum c 90 dipnew = dl if (du .gt. dl) dipnew = du if (dip .lt. dipnew) dip = dipnew low = gcm(ig) high = lcm(ih) c c--- recycle c icount=icount+1 if(icount.gt.100) then dip=-999 xl=-999 xu=-999 return endif goto 40 c 100 dip = .5 * dip xl = x(low) xu = x(high) return end c------------------------------------------------------------------------------ SUBROUTINE EDFGOF (N, X, ITEST, XBAR, S2, D, V, W2, W2P, U2, + U2P, A2, A2P, IFAULT) c------------------------------------------------------------------------------ c--- Algorithm AS 248.1 Applied Statistics (91989), Vol. 38. No. 3 c c Tests a sample for uniformity, normality, or exponentiality c using goodness-of-fit statistics based on the empirical c distribution function parameter (iii = 10000) integer i, ifault, itest, j, n real x(n), z(III), xbar, s2, d, v, w2, w2p, u2, u2p, + a2, a2p, rn, + rootn, rn2, zero, pt01, pt05, pt1, pt11, pt12, pt155, + pt16, pt2, pt24, pt26, pt3, pt35, pt4, pt5, pt6, pt75, pt8, + pt82, pt85, one, twop25, wcut2(3), wcoef2(4,3), wcut3(3), + wcoef3(4,3), ucut2(3), ucoef2(4,3), ucut3(3), ucoef3(4,3), + acut2(3), acoef2(4,3), acut3(3), acoef3(4,3) data zero, pt01, pt1, pt05, pt11, pt12, pt155, pt16, pt2, + pt24 / 0.0, 0.01, 0.05, 0.1, 0.11, 0.12, 0.155, 0.16, 0.2, + 0.24/ data pt26, pt3, pt35, pt4, pt5, pt6, pt75, pt8, pt82, pt85, + one / 0.26, 0.3, 0.35, 0.4, 0.5, 0.6, 0.75, 0.8, 0.82, 0.85, + 1.0 / data twop25 / 2.25 / data wcut2, wcut3, ucut2, ucut3, acut2, acut3 / 0.0275, 0.051, + 0.092, 0.035, 0.074, 0.160, 0.0262, 0.048, 0.094, 0.029, + 0.062, 0.120, 0.200, 0.340, 0.600, 0.260, 0.510, 0.950 / data ((wcoef2(i,j), j=1,3), i=1,4) / -13.953, 775.5, -12542.61, + -5.903, 179.546, -1515.29, 0.886, -31.62, 10.897, 1.111, + -34.242, 12.832 / data ((wcoef3(i,j), j=1,3), i=1,4) / -11.334, 459.098, -5652.1, + -5.779, 132.89, -866.58, 0.586, -17.87, 7.417, 0.447, + -16.592, 4.849 / data ((ucoef2(i,j), j=1,3), i=1,4) / -13.642, 766.31, + -12432.74, -6.3328, 214.57, -2022.28, 0.8510, -32.006, + -3.45, 1.325, -38.918, 16.45 / data ((ucoef3(i,j), j=1,3), i=1,4) / -11.703, 542.5, -7574.59, + -6.3288, 178.1, -1399.49, 0.8071, -25.166, 8.44, 0.7663, + -24.359, 4.539 / data ((acoef2(i,j), j=1,3), i=1,4) / -13.436, 101.14, -223.73, + -8.318, 42.796, -59.938, 0.9177, -4.279, -1.38, 1.2937, + -5.709, 0.0186 / data ((acoef3(i,j), j=1,3), i=1,4) / -12.2204, 67.459, -110.3, + -6.1327, 20.218, -18.663, 0.9209, -3.353, 0.300, 0.731, + -3.009, 0.15 / ifault = 1 if (itest .lt. 1 .or. itest .gt. 3) return ifault = 2 if ( n .lt. 2 .or. (itest .eq. 2 .and. n .eq. 2)) return ifault = 3 do 10 i= 2, n if (x(i) .lt. x(i-1)) return 10 continue rn = n rootn = sqrt(rn) rn2 = rn * rn if(itest .eq. 1) then c C Test for uniformity (F(X) is completely specified) c ifault = 0 CALL STATS (N, X, D, V, W2, U2, A2, IFAULT) if (ifault .ne. 0) return C C Modifications when F(X) is completely specified C d = d * (rootn + pt12 + pt11 / rootn) v = v * (rootn + pt155 + pt24 / rootn) w2 = (w2 - pt4 / rn + pt6 / rn2) * (one + one / rn) u2 = (u2 - pt1 / rn + pt1 / rn2) * (one + pt8 / rn) return else C C Estimate the mean by XBAR C xbar = zero do 20 i = 1, n xbar = xbar + x(i) 20 continue xbar = xbar / rn if (itest .eq. 2) then C C Test for normality (MU and SIGMA**2 unspecified) C First estimate the variance by S**2 C ifault = 5 s2 = zero do 30 i = 1, n s2 = s2 + (x(i) - xbar) ** 2 30 continue s2 = s2 / (rn - one) if (s2 .le. zero) return C C Compute Z(I) = F((XI)-XBAR)/S) C s = sqrt(s2) do 40 i = 1, n z(i) = ALNORM((X(I) - XBAR) / S, .FALSE.) if (z(i).le.0) z(i)=1.E-07 if (z(i).ge.1) z(i)=1.-1.E-07 40 continue CALL STATS (N, Z, D, V, W2, U2, A2, IFAULT) if (ifaulT .NE. 0) return C C Modifications when F(X) is the normal distribution C d = d * (rootn - pt01 + pt85 / rootn) v = v * (rootn + pt05 + pt82 / rootn) w2 = w2 * (one + pt5 / rn) w2p=0 if (w2.gt.2) goto 41 !skip PVALUE calc. w2p = PVALUE (W2, WCUT2, WCOEF2) 41 u2 = u2 * (one + pt5 / rn) u2p=0 if (u2.gt.2) goto 42 !skip PVALUE calc. u2p = PVALUE (U2, UCUT2, UCOEF2) 42 a2 = a2 * (one + pt75 / rn + twop25 / rn2) a2p=0 if (a2.gt.2) return a2p = PVALUE (A2, ACUT2, ACOEF2) return else C C Test for exponentiality (scale parameter unspecified) C ifault = 6 if (xbar .le. zero) return do 50 i = 1, n z(i) = one - exp(-x(i) / xbar) 50 continue CALL STATS (N, Z, D, V, W2, U2, A2, IFAULT) if (ifault .ne. 0) return C C Modifications when F(X) is the exponential distribution C d = (d - pt2 / rn) * (rootn + pt26 + pt5 / rootn) v = (v - pt2 / rn) * (rootn + pt24 + pt35 / rootn) w2 = w2 * (one + pt16 / rn) w2p = PVALUE (W2, WCUT3, WCOEF3) u2 = u2 * (one + pt16 / rn) u2p = PVALUE (U2, UCUT3, UCOEF3) a2 = a2 * (one + pt3 / rn) a2p = PVALUE (A2, ACUT3, ACOEF3) return endif endif end c------------------------------------------------------------------------------ SUBROUTINE GAPPER (XDATA,N,XSGAP,ZSTAR,NBIG,JBIG,ZBIG,ZP,TP) c------------------------------------------------------------------------------ c--- This routine returns a list of gaps in a data set, weighted by c their location with respect to the middle of the data. It also c obtains values for these gaps relative to their average size. This c normalized gap provides a measure of how large a given gap is c compared to the rest of the gaps. Wainer and Schacht (1978; c Psychometrika 43, 203) assess the liklihood that a gap of given c size could occur in a normal distribution. c ZSTAR -- weighted and normalized gaps c NBIG -- number of "big" gaps (> 2.25) c JBIG -- indices of lower data value giving rise c to a "big" gap c ZBIG -- the size of the "big" gap c ZP -- the "per gap" probability of occurence c TP -- the "distribution" probability of finding this c many (or more) gaps at this significance level c Thr routine also reports an estimate of scale , XSGAP, based on c the sum of the weighted gaps (see Wainer and Thissen 1976; c Psychometrika 41, 9.) c*************************************************************************** parameter (iii = 10000) data pi /3.1415927/ real xdata(n),zstar(n),ygap(iii),zbig(100),zp(100),tp(100) real big(iii) integer jbig(iii) c--- determine gap, weight, weighted gap vectors, and sum of gaps ysum = 0. do 10 i=1,n-1 gap = xdata(i+1) - xdata(i) weight = i*(n-i) ygap(i) = (gap*weight)**0.5 big(i) = i ysum = ysum + gap*weight 10 continue c--- obtain scale estimator using sum of gaps xsgap = sqrt(pi)/(n*(n-1))*ysum c--- set flag so only enter this loop once c--- determine mid-mean of the weighted gaps if(iflag.eq.0) then iflag = 1 CALL SORT2 (N-1,YGAP,BIG) CALL XMIDMEAN (YGAP,N-1,GMID) c--- obtain standardized weighted gaps if(gmid.eq.0) gmid = 1 do 20 i=1,n-1 zstar(i) = ygap(i)/gmid 20 continue c--- save gaps larger than 2.25 as possibly significant j=0 do 30 i=1,n-1 if(zstar(i).ge.2.25) then j=j+1 zbig(j)=zstar(i) ! proposed significant gap jbig(j)=ifix(big(i)) ! original index of gap endif 30 continue nbig = j ! number of big gaps c--- obtain probability of gap occurence iz225 = 0 iz250 = 0 iz275 = 0 iz300 = 0 iz325 = 0 iz350 = 0 do 40 i=1,nbig if(zbig(i).ge.2.25.and.zbig(i).lt.2.50) then zp(i) = 0.03 iz225 = iz225+1 else if(zbig(i).ge.2.50.and.zbig(i).lt.2.75) then zp(i) = 0.014 iz250 = iz250+1 else if(zbig(i).ge.2.75.and.zbig(i).lt.3.00) then zp(i) = 0.006 iz275 = iz275+1 else if(zbig(i).ge.3.00.and.zbig(i).lt.3.25) then zp(i) = 0.002 iz300 = iz300+1 else if(zbig(i).ge.3.25.and.zbig(i).lt.3.50) then zp(i) = 0.001 iz325 = iz325+1 else if(zbig(i).ge.3.50) then zp(i) = 0.0005 iz350 = iz350+1 endif 40 continue c--- Obtain number of gaps which are significant at various levels iz225 = iz225+iz250+iz275+iz300+iz325+iz350 iz250 = iz250+iz275+iz300+iz325+iz350 iz275 = iz275+iz300+iz325+iz350 iz300 = iz300+iz325+iz350 iz325 = iz325+iz350 if(iz225.ne.0) then tp225 = BETAI(float(IZ225),float(N-IZ225),0.03) else tp225 = -9.99 endif if(iz250.ne.0) then tp250 = BETAI(float(IZ250),float(N-IZ250),0.014) else tp250 = -9.99 endif if(iz275.ne.0) then tp275 = BETAI(float(IZ275),float(N-IZ275),0.006) else tp275 = -9.99 endif if(iz300.ne.0) then tp300 = BETAI(float(IZ300),float(N-IZ300),0.002) else tp300 = -9.99 endif if(iz325.ne.0) then tp325 = BETAI(float(IZ325),float(N-IZ325),0.001) else tp325 = -9.99 endif if(iz350.ne.0) then tp350 = BETAI(float(IZ350),float(N-IZ350),0.0005) else tp350 = -9.99 endif do 50 i = 1,nbig if(zbig(i).ge.2.25.and.zbig(i).lt.2.50) then tp(i) = tp225 else if(zbig(i).ge.2.50.and.zbig(i).lt.2.75) then tp(i) = tp250 else if(zbig(i).ge.2.75.and.zbig(i).lt.3.00) then tp(i) = tp275 else if(zbig(i).ge.3.00.and.zbig(i).lt.3.25) then tp(i) = tp300 else if(zbig(i).ge.3.25.and.zbig(i).lt.3.50) then tp(i) = tp325 else if(zbig(i).ge.3.50) then tp(i) = tp350 endif 50 continue endif return end c------------------------------------------------------------------------------ SUBROUTINE GCF(GAMMCF,A,X,GLN) c------------------------------------------------------------------------------ c--- taken from Numerical Recipes c****************************************************************************** parameter (itmax=100,eps=3.e-7) GLN=GAMMLN(A) gold=0. a0=1. a1=x b0=0. b1=1. fac=1. do 11 n=1,itmax an=float(n) ana=an-a a0=(a1+a0*ana)*fac b0=(b1+b0*ana)*fac anf=an*fac a1=x*a0+anf*a1 b1=x*b0+anf*b1 if(a1.ne.0.)then fac=1./a1 g=b1*fac if(abs((g-gold)/g).lt.eps)go to 1 gold=g endif 11 continue PAUSE 'A too large, ITMAX too small' 1 exp1=-x+a*alog(x)-gln if (exp1.lt.-50.) exp1=-50. gammcf=exp(exp1)*g return end c------------------------------------------------------------------------------ SUBROUTINE GSER(GAMSER,A,X,GLN) c------------------------------------------------------------------------------ c--- taken from Numerical Recipes c****************************************************************************** parameter (itmax=100,eps=3.e-7) GLN=GAMMLN(A) if(x.le.0.) then if(x.lt.0.) pause gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,itmax ap=ap+1. del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*eps)go to 1 11 continue pause 'a TOO LARGE, itmax TOO SMALL' 1 gamser=sum*exp(-x+a*log(x)-gln) return end c------------------------------------------------------------------------------ SUBROUTINE JACKNIFE (XDATA,N,IP,X1LJACK,X1SJACK, + X1LLJACK,X1LSJACK,X2LJACK, + X2SJACK,X2LLJACK,X2LSJACK) c------------------------------------------------------------------------------ c--- given an array XDATA of length N, and external function PARAM, c this routine calculates the expected error in the estimate c of the supplied function via the statistical jacknife. In other c words, it obtains the spread in the parameter when calculated by c dropping one data point at a time. For reference, c see the article by Diaconis and Efron in Scientific American. c c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) integer n,nletter real xdata(n),xout(iii),x1jparam(iii),x1p(iii),x1lp(iii) real x2jparam(iii),x2p(iii),x2lp(iii) real depths(15) dimension sletter(15),rletter(8) c--- find the parameter considering all values if (ip.eq.1 .or. ip.eq.2) then ! mean and stand dev x1all = PARAM(XDATA,N,1) x2all = PARAM(XDATA,N,2) else if (ip.eq.6) then call XBIWT(XDATA,N,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) ! BIWEIGHT x1all = xlbiwt x2all = xsbiwt else if (ip.eq.7) then ! gapper x1all = PARAM(XDATA,N,IP) x2all = x1all endif if(x1all.gt.0) then x1lall = alog10(x1all) else x1lall = 0 endif if(x2all.gt.0) then x2lall = alog10(x2all) else x2lall = 0 endif c--- create the pseudovalues data array with the missing value identified c by setting it to -32368., then calling PARAM, and repeating process do 10 i=1,n temp=xdata(i) ! store original value in temp xdata(i)=-32368. ! replace with dropped value if (ip.eq.1 .or. ip.eq.2) then x1jparam(i)=PARAM(XDATA,N,1) x2jparam(i)=PARAM(XDATA,N,2) else if (ip.eq.6) then CALL UPDATE(N,XDATA,M,XOUT) CALL XBIWT(XOUT,M,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) x1jparam(i) = xlbiwt ! calculate parameter value x2jparam(i) = xsbiwt else if (ip.eq.7) then x1jparam(i)=PARAM(XDATA,N,IP) x2jparam(i)=x1jparam(i) endif if(x1jparam(i).gt.0) x1 = alog10(x1jparam(i)) if(x2jparam(i).gt.0) x2 = alog10(x2jparam(i)) x1p(i) = n*x1all - (n-1)*x1jparam(i) ! pseudovalue x2p(i) = n*x2all - (n-1)*x2jparam(i) x1lp(i) = n*x1lall - (n-1)*x1 ! log of pseudovalue x2lp(i) = n*x2lall - (n-1)*x2 xdata(i)=temp ! restore original value 10 continue c--- obtain estimate of central location and spread x1sum = 0 x2sum = 0 x1xsum1 = 0 x2xsum1 = 0 x1lsum = 0 x2lsum = 0 x1xlsum1 = 0 x2xlsum1 = 0 x1xsum2 = 0 x2xsum2 = 0 x1xlsum2 = 0 x2xlsum2 = 0 do 20 i=1,n x1sum = x1sum + x1p(i) x2sum = x2sum + x2p(i) x1xsum1 = x1xsum1 + x1p(i)*x1p(i) x2xsum1 = x2xsum1 + x2p(i)*x2p(i) x1lsum = x1lsum + x1lp(i) x2lsum = x2lsum + x2lp(i) x1xlsum1 = x1xlsum1 + x1lp(i)*x1lp(i) x2xlsum1 = x2xlsum1 + x2lp(i)*x2lp(i) 20 continue x1xsum2 = x1sum*x1sum/n x2xsum2 = x2sum*x2sum/n x1xlsum2 = x1lsum*x1lsum/n x2xlsum2 = x2lsum*x2lsum/n x1ljack = x1sum/n x2ljack = x2sum/n x1sjack = (abs(x1xsum1-x1xsum2)/(n*(n-1)))**0.5 x2sjack = (abs(x2xsum1-x2xsum2)/(n*(n-1)))**0.5 x1lljack = x1lsum/n x2lljack = x2lsum/n x1lsjack = (abs(x1xlsum1-x1xlsum2)/(n*(n-1)))**0.5 x2lsjack = (abs(x2xlsum1-x2xlsum2)/(n*(n-1)))**0.5 return end c------------------------------------------------------------------------------ SUBROUTINE KSREJ (DATA,N,KSAVE,KSDEV,NTOT) c------------------------------------------------------------------------------ c--- This routine performs a "K-S" clipping of the data, similar to the c cleaning technique described in Yahil and Vidal (1977). The procedure c first checks (with the K-S test) whether the sample is consistent with c a Gaussian distribution. If not the procedure then deletes c the data point most deviant from the sample mean, and redetermines the c sample mean and standard deviation. Each point is tested in turn c in order of their apparent deviation, until no further rejection occurs. c The list is now considered "clean," and the final estimates of mean and c standard deviation are accepted. c c--- This routine will exit back to the main program if the initial number of c data points N <= 5, or if the sample is clipped to NTOT=5 c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) dimension data(n),tdata(iii),sdata(iii),ssdata(iii) real ksave,ksdev if(n.le.5) then ! return to main routine if input n <= 5 ksave=-999 ksdev=-999 ntot = n return endif c--- copy the original data set do 5 i=1,n sdata(i) = data(i) 5 continue c--- find intial estimate of sample mean and deviation ave = PARAM(SDATA,N,1) ! find new mean sdev = PARAM(SDATA,N,2) ! find new sdev nout=0 c--- test once to see if all is OK ITEST=2 ! FOR NORMAL DISTRIBUTION CALL EDFGOF(N,XDATA,ITEST,XBAR,S2,D,V,W2,W2P,U2,U2P, + A2,A2P,IFAULT) if(d.le..775) test=.25 if(d.gt..775.and.d.le..819) test=.15 if(d.gt..819.and.d.le..895) test=.10 if(d.gt..895.and.d.le..995) test=.05 if(d.gt..995.and.d.le.1.035) test=.025 if(d.gt.1.035) test=.01 if (test.gt.0.20) then ! if consistent, return now ksave = ave ksdev = sdev ntot = n return endif c--- create a list of data ordered by their distance from the sample mean 20 continue do 10 i=1,n-nout tdata(i) = abs(sdata(i)-ave) 10 continue CALL SORT2 (N-NOUT,TDATA,SDATA) ! sorts TDATA hold SDATA c--- drop first deviant point temp = sdata(n-nout) sdata(n-nout) = -32368. ave = PARAM(SDATA,N-NOUT,1) ! find new mean sdev = PARAM(SDATA,N-NOUT,2) ! find new sdev c test if new mean and standard deviation give an acceptable Gaussian fit c if so, then delete offending point from consideration CALL UPDATE (N-NOUT,SDATA,M,SSDATA) CALL SORT (M,SSDATA) if(m.eq.5) then ! get out of here ksave=-999 ksdev=-999 ntot=5 return endif c--- test to see if consistent with a Gaussian ITEST=2 ! FOR NORMAL DISTRIBUTION CALL EDFGOF(M,SSDATA,ITEST,XBAR,S2,D,V,W2,W2P,U2,U2P, + A2,A2P,IFAULT) if(d.le..775) test=.25 if(d.gt..775.and.d.le..819) test=.15 if(d.gt..819.and.d.le..895) test=.10 if(d.gt..895.and.d.le..995) test=.05 if(d.gt..995.and.d.le.1.035) test=.025 if(d.gt.1.035) test=.01 if(test.le.0.20) then ! delete from list nout = nout + 1 else ksave=ave ksdev=sdev ntot=n-nout-1 return endif goto 20 end c------------------------------------------------------------------------------ SUBROUTINE KURTSKEW (N,SKEW,KURT,Z1,Z2,K2) c------------------------------------------------------------------------------ c--- This subroutine obtains the statistics z1 and z2, which are c normally distributed under the hypothesis of normality of c the underlying population. It also returns the statistic k2, c which is distributed like Chi^2 when the underlying population c is normally distributed. c c The user inputs the sample skewness, SKEW, and the sample kurtosis, c KURT. For details of the calculation see D'Agostino, Belanger, c and D'Agostino (1990), American Statistician 44, 316. c c For values of n<9 this should not be used c c***************************************************************************** implicit real*4 (a-h,o-z) real*4 k2,kurt,nn integer n nn = n if (nn.lt.9) return ! this test is not recommended so leave c--- This is the skewness test y = skew*((nn+1)*(nn+3)/(6*(nn-2)))**0.5 beta2 = 3*(nn*nn+27*nn-70)*(nn+1)*(nn+3)/((nn-2)*(nn+5) + *(nn+7)*(nn+9)) w = (-1.+(2.*(beta2-1.))**0.5)**0.5 d = 1./(log(w))**0.5 a = (2./(w**2-1.))**0.5 z1 = d*log(y/a+((y/a)**2+1.)**0.5) c--- This is the kurtosis test e = 3*(nn-1)/(nn+1) var =24*nn*(nn-2)*(nn-3)/((nn+1)**2*(nn+3)*(nn+5)) x = (kurt-e)/(var)**0.5 beta1 = (6*(nn**2-5*nn+2)/((nn+7)*(nn+9))*(6* + (nn+3)*(nn+5)/(nn*(nn-2)*(nn-3)))**0.5)**2 A = 6.+8./(beta1**0.5)*(2./(beta1**0.5)+(1.+4./beta1)**0.5) exp1=(1.-(2./A))/(1.+x*(2./(A-4.))**0.5) exp2=abs(exp1)**0.333 if (exp1.lt.0.) exp2=-exp2 z2 = (1.-(2./9./A)-exp2)/ + (2./(9.*A))**0.5 c--- This is the omnibus test k2 = (z1)**2+(z2)**2 return end c------------------------------------------------------------------------------ SUBROUTINE LETTERS (XDATA,N,DEPTHS,XLV,MLV,SLV,RLV,CLV,NLV) c------------------------------------------------------------------------------ c--- NOTE -- this routine also sorts the data low to high ! c--- For the batch of values in XDATA, find the selected quantiles known c as the letter values. Upon exit, XLV contains the letter values, c D contains the corresponding depths, and NLV is the number of pairs c of letter values. Specifically, XLV(1,1) and XLV(1,2) are both c set equal to the median, whose depth DEPTHS(1), is (N+1)/2. The rest of c the letter values come in pairs and are stored in XLV in order from c the hinges out to the extremes. Thus XLV(2,1) and XLV(2,2) are the c lower hinge and upper hinge, respectively, and XLV(NLV,1) and XLV(NLV,2) c are the lower extreme (minimum) and upper extremem (maximum), respectively. c c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) integer n,nlv,i,j,k,pt1,pt2 dimension xdata(n),depths(15),xlv(15,2),slv(15),glv(8),rlv(8) real mlv(8),clv(8),cclv(iii,8) data glv/1.,1.349,2.301,3.068,3.726,4.308,4.836,5.320/ c--- read in correction factors if (n.gt.100) then xn=float(n) cclv(n,1)=1.-(2./xn)*(1.-.86*((xn/2.)-(n/2))) cclv(n,2)=1.71-(3./xn)*(1.-.8*((xn/4.)-(n/4))) cclv(n,3)=2.27-(5./xn)*(1.-.77*((xn/8.)-(n/8))) cclv(n,4)=2.76-(10./xn)*(1.-.76*((xn/16.)-(n/16))) cclv(n,5)=3.19-(15./xn)*(1.-.76*((xn/32.)-(n/32))) cclv(n,6)=3.58-(30./xn)*(1.-.75*((xn/64.)-(n/64))) cclv(n,7)=3.94-(55./xn)*(1.-.75*((xn/128.)-(n/128))) do 10 i=1,7 cclv(n,i)=1.349*cclv(n,i) 10 continue cclv(n,8)=1.0 else do 11 i=1,100 read (8,*) id,(cclv(i,j),j=1,7) cclv(i,8) = 1.0 11 continue rewind (unit=8) endif c--- sort the data low to high CALL SORT (N,XDATA) c--- handle median separately because it is not a pair of lv depths(1) = float(n+1)/2.0 j = (n/2) + 1 pt2 = n + 1 - j xlv(1,1) = (xdata(j) + xdata(pt2))/2.0 xlv(1,2) = xlv(1,1) k=n i=2 20 k= (k+1)/2 j = (k/2) + 1 depths(i) = float(k+1)/2.0 pt2 = k + 1 - j xlv(i,1) = (xdata(j) + xdata(pt2))/2.0 pt1 = n-k+j pt2 = n+1-j xlv(i,2) = (xdata(pt1) + xdata(pt2))/2.0 i = i + 1 if(depths(i-1) .gt. 2.0) goto 20 nlv = i if(nlv.gt.8) nlv=8 depths(i) = 1.0 xlv(i,1) = xdata(1) xlv(i,2) = xdata(n) c--- obtain list of letter mids,spreads, and ratios to Gaussian values slv(1)=0 rlv(1)=0 do 30 i=2,min(8,nlv) slv(i)= abs(xlv(i,1)-xlv(i,2)) rlv(i) = slv(i)/glv(i) mlv(i) = (xlv(i,1)+xlv(i,2))/2. 30 continue c--- obtain corrections to spreads as a function of n do 40 i=2,min(8,nlv) clv(i) = slv(i)/cclv(n,i-1) 40 continue return end c------------------------------------------------------------------------------ SUBROUTINE MCCON (X,N,CONF,XAVE,XSDEV,XLBIWT,XSBIWT, + MCCONF,NSIMS) c------------------------------------------------------------------------------ c--- This routine obtains the BOOTSTRAP confidence c intervals for an input parameter c***************************************************************************** parameter (iii = 10000) real msim(iii),ssim(iii),bsim(iii),sbsim(iii) real mcave,mcsdev,conf(4) real x(n),xout(iii),mcconf(32,4),xnew(iii),freq(iii) real infav(iii),infsd(iii),infbt(iii),infsbt(iii) logical lsb,lbb c--- obtain sampling distribution ! sample mean and stand dev CALL BOOTSTRAP (X,N,1,NSIMS,MSIM,SSIM) ! biweight location and spread CALL BOOTSTRAP (X,N,6,NSIMS,BSIM,SBSIM) c c-- calculate the acceleration constants for the data set c using the Jacknife Estimation of Frangos and Schucany c Comp. Stat. & Data Anal. 9 (1990) 271-281 c do 15 i=1,n temp=x(i) ! store original value in temp x(i)=-32368. ! replace with dropped value CALL UPDATE(N,X,M,XOUT) xnave = PARAM(XOUT,M,1) xnsdev = PARAM(XOUT,M,2) CALL XBIWT(XOUT,M,XNLBIWT,XNSBIWT,XNLBIWT1,XNSBIWT1) infav(i) = (n-1)*(xave-xnave) infsd(i) = (n-1)*(xsdev-xnsdev) infbt(i) = (n-1)*(xlbiwt-xnlbiwt) infsbt(i) = (n-1)*(xsbiwt-xnsbiwt) x(i) = temp ! reset original data point 15 continue snum1=0 sden1=0 snum2=0 sden2=0 snum3=0 sden3=0 snum4=0 sden4=0 do 16 i=1,n snum1= snum1 + infav(i)**3 sden1= sden1 + infav(i)**2 snum2= snum2 + infsd(i)**3 sden2= sden2 + infsd(i)**2 snum3= snum3 + infbt(i)**3 sden3= sden3 + infbt(i)**2 snum4= snum4 + infsbt(i)**3 sden4= sden4 + infsbt(i)**2 16 continue aave=snum1/(6*(sden1**1.5)) asdev=snum2/(6*(sden2**1.5)) abwt=snum3/(6*(sden3**1.5)) asbwt=snum4/(6*(sden4**1.5)) c--- obtain monte-carlo estimators of intervals c--- get confidence iunterval of the mean from MSIM do 10 i=1,4 lsb=.true. sbound=0.0 lbb=.true. bbound=0.0 ifault = 0 CALL MONTE (XAVE,CONF(i),NSIMS,AAVE,MSIM,LSB, + SBOUND,LBB,BBOUND,MCCONF(1,i),MCCONF(2,i),MCCONF(3,i), + MCCONF(4,i),MCCONF(5,i),MCCONF(6,i),MCCONF(7,i), + MCCONF(8,i),IFAULT) C-- get biweight conf int on location from BSIM lsb=.true. sbound=0.0 lbb=.true. bbound=0.0 ifault = 0 CALL MONTE (XLBIWT,CONF(i),NSIMS,ABWT,BSIM,LSB, + SBOUND,LBB,BBOUND,MCCONF(9,i),MCCONF(10,i),MCCONF(11,i), + MCCONF(12,i),MCCONF(13,i),MCCONF(14,i),MCCONF(15,i), + MCCONF(16,i),IFAULT) c--- get confidence interval of the stand dev from SSIM lsb=.true. sbound=0.0 lbb=.true. bbound=0.0 ifault = 0 CALL MONTE (XSDEV,CONF(i),NSIMS,ASDEV,SSIM,LSB, + SBOUND,LBB,BBOUND,MCCONF(17,i),MCCONF(18,i),MCCONF(19,i), + MCCONF(20,i),MCCONF(21,i),MCCONF(22,i),MCCONF(23,i), + MCCONF(24,i),IFAULT) c-- get biweight spread conf limits from SBSIM lsb=.true. sbound=0.0 lbb=.true. bbound=0.0 ifault = 0 CALL MONTE (XSBIWT,CONF(i),NSIMS,ASBWT,SBSIM,LSB, + SBOUND,LBB,BBOUND,MCCONF(25,i),MCCONF(26,i),MCCONF(27,i), + MCCONF(28,i),MCCONF(29,i),MCCONF(30,i),MCCONF(31,i), + MCCONF(32,i),IFAULT) 10 continue return end c------------------------------------------------------------------------------ SUBROUTINE MDIAN1 (X,N,XMED) c------------------------------------------------------------------------------ c--- Taken from Numerical Recipes, page 460. c Given an array X of N numbers, returns their median value XMED, and c the array is sorted low to high c c******************************************************************************* real x(n) call sort(n,x) n2=n/2 if (2*n2.eq.n) then xmed=0.5*(x(n2)+x(n2+1)) else xmed=x(n2+1) endif return end c------------------------------------------------------------------------------ SUBROUTINE MOMENT(DATA,N,AVE,ADEV,SDEV,VAR,SKEW,CURT) c------------------------------------------------------------------------------ c--- given an array DATA of length N, this routine returns its c mean AVE, average deviation ADEV, standard deviation SDEV, c variance VAR, skewness SKEW, and kurtosis CURT c c--- stolen (unabashedly) from NUMERICAL RECIPES c c**************************************************************************** implicit real*4 (a-h,o-z) dimension data(n) if(n.le.1)pause 'n MUST BE AT LEAST 2' s=0. do 11 j=1,n s=s+data(j) 11 continue ave=s/n adev=0. var=0. skew=0. curt=0. do 12 j=1,n s=data(j)-ave adev=adev+abs(s) p=s*s var=var+p p=p*s skew=skew+p p=p*s curt=curt+p 12 continue adev=adev/n var=var/(n-1) sdev=sqrt(var) if(var.ne.0.)then skew=skew/float(n)/(var*(n-1)/float(n))**1.5 curt=curt/float(n)/(var*(n-1)/float(n))**2 else skew = 0.0 curt = 0.0 endif return end c------------------------------------------------------------------------------ SUBROUTINE MONTE(EST,CONFI,NSIMS,ACC,SIMVAL,LSB, + SBOUND,LBB,BBOUND,ALOW,AHI,BLOW,BHI, + CLOW,CHI,DLOW,DHI,IFAULT) c------------------------------------------------------------------------------ c c--- Algorithm AS 214 Appl. Stat. (1985) Vol. 34, No.3 c c--- Sets up Monte Carlo confidence intervals c uses function PPND - AS 111 c uses function ALNORN -- AS 66 c c--- The user supplies the estimate EST of a parameter, and MONTE c obtains the appropriate confidence intervals c c ALOW: lower conf. limit, symmetric c AHI: upper conf. limit, symmetric c BLOW: lower conf. limit, percentile c BHI: upper conf. limit, percentile c CLOW: lower conf. limit, bias corrected c CHI: upper conf. limit, bias corrected c DLOW: lower conf. limit, bias corrected accelerated c DHI: upper conf. limit, bias corrected accelerated c c*************************************************************************** logical lsb,lbb dimension simval(nsims) data least/5/ data zero,half,one,two,hun,big + /0.0E0,0.5E0,1.0E0,2.0E0,1.0E2,1.e20/ ifault = 0 sbnd = sbound bbnd = bbound if (confi .ge. hun) ifault = 6 if (confi .le. zero) ifault = 7 if (ifault .gt. 0) return c c--- symmetric mci c c--- find mean and variance c v1 = zero v2 = zero do 10 j = 1,nsims fj = float(j) v2 = v2 + (fj - one)*(simval(j) - v1)**2/fj v1 = (simval(j) + (fj - one)*v1)/fj 10 continue stderr = zero if (v2 .gt. zero) stderr = sqrt(v2/float(nsims - 1)) alpha = half*(hun - confi)/hun z = -ppnd(alpha,ifault) if (ifault .eq. 0) goto 20 ifault = 6 return 20 ahi = z*stderr alow = est - ahi ahi = est + ahi c c--- calculate bias adjustment, so that the number c of values that must be ordered is known c call biasad(est,nsims,acc,simval,z,limit1,limit2, + limita1,limita2,ifault) limitl = alpha*float(nsims + 1) + half limitu = (one - alpha)*float(nsims + 1) + half if (limitl .lt. least) ifault = 5 if (limitl .eq. limitu) ifault = 7 if (ifault .gt. 0) return l1 = max0(limit1, 2*limitl) l2 = min0(limit2, 2*limitu - nsims - 1) if (lsb) sbnd = -big if (lbb) bbnd = big c c--- select and order the l1 smallest and the nsims+1-l2 largest c values c call sort(nsims,simval) if (sbnd .le. simval(1)) goto 40 if (lsb) goto 30 ifault = 1 return 30 sbnd = simval(1)*two if (sbnd .gt. zero) sbnd = -sbnd 40 if (bbnd .ge. simval(nsims)) goto 60 if (lbb) goto 50 ifault = 2 return 50 bbnd = simval(nsims)*two if (bbnd .lt. zero) bbnd = -bbnd c c--- equal tails mci (percentile method) c 60 if (limitl .gt. 0) blow = simval(limitl) if (limitl .le. 0) blow = simval(1) if (limitu .le. nsims) bhi = simval(limitu) if (limitu .gt. nsims) bhi = simval(nsims) c c--- bias-corrected percentile method c clow = sbnd if (limit1 .gt. 0) clow = simval(limit1) if (limit1.le.0) clow = simval(1) chi = bbnd if (limit2 .le. nsims) chi = simval(limit2) if (limit2.gt.nsims) chi = simval(nsims) c-- accelerated bias corrected dlow = sbnd if (limita1 .gt. 0) dlow = simval(limita1) if (limita1.le.0) dlow = simval(1) dhi = bbnd if (limita2 .le. nsims) dhi = simval(limita2) if (limita2.gt.nsims) dhi = simval(nsims) return end c------------------------------------------------------------------------------ SUBROUTINE NORMALITY (XDATA,N,XAVE,XSDEV,XSKEW,XCURT, + XSBIWT,XMED,TEST) c------------------------------------------------------------------------------ c--- This routine obtains the tests for non-normality of a data sample. c A number of tests are performed. In order, these are: c c (1) The a- , u-, w- tests referred to by Yahil c and Vidal (1977) (Ap.J. 214, 347) with an improved W-test c calculation due to Royston (A.S. 181) -- critical values for c a and u in REFS, critical value for w calculated and returned c c (2) The b1 and b2 tests (essentially a renormalized set of third and c fourth moments) -- probability given in D'Agostino, Belanger, c D'Agostino (1990), The American Statistician, 44, 316-321 c c (3) The i-test described by Iglewicz (1983) (UREDA, p.426) -- critical c values are calculated and returned c c (4) The DIP statistic test for unimodality (Hartigan A.S. 217) -- c critical values given in Hartigan and Hartigan (1985) Ann. Stat. c 13, 70. c c (5) A modified K-S statistic -- see Goodness of Fit (D-Agistino and c Stephens, Chp. 4) c c (6) A modified Kuiper statistic -- see Goodness of Fit (D-Agistino and c Stephens, Chp. 4) c c (7) A modified Cramer Von-Mises -- see Goodness of Fit (D-Agistino and c Stephens, Chp. 4) c c (8) A modified Watson test -- see Goodness of Fit (D-Agistino and c Stephens, Chp. 4) c c (9) A modified Anderson-Darling test --see Goodness of Fit (D-Agistino and c Stephens, Chp. 4) c c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) dimension xdata(n),a(iii) real test(26),i90,i95,istat,isum,ntilda,max,k2 integer mn(iii),mj(iii),gcm(iii),lcm(iii) logical upper if (n.lt.5) return do 9 i=1,26 9 test(i) = 0.0 c---------------------- TEST 1 -- A, U, W --------------------- a1sum=0.0 do 10 i=1,n ! a-test a1sum = a1sum + abs(xdata(i)-xave) 10 continue a1 = a1sum/(n*xsdev) test(1) = a1 u1 = (xdata(n)-xdata (1))/xsdev ! u-test test(2) = u1 c--- a slightly different version of the W-test n2=n/2 ifault = 0 ssq = 0. c--- obtain SSQ do 19 i=1,n 19 ssq = ssq+(xdata(i)-xave)**2 CALL WCOEF (A,N,N2,EPS,IFAULT) CALL WEXT (XDATA,N,SSQ,A,N2,EPS,W,PW,IFAULT) test(3) = w test(4) = pw c---------------------- TEST 2 -- b1,b2 TESTS ----------------------- if (n.lt.9) then ! test should not be used for n<9 test(23)=-999 test(24)=-999 test(25)=-999 test(26)=-999 goto 323 endif CALL KURTSKEW (N,XSKEW,XCURT,ZB1,ZB2,K2) test(23)=ALNORM(ZB1,.true.) ! 2x for the 2-sided test if (xskew.lt.0.) test(23)=1.-test(23) ! gives P in the left tail test(24)=ALNORM(ZB2,.true.) ! 2x for the 2-sided test if (xcurt.lt.3.) test(24)=1.-test(24) ! gives P in the left tail test(25)=K2 test(26)=GAMMQ(1.,0.5*K2) ! 2 degrees of freedom 323 test(5) = xskew test(6) = xcurt c---------------------- TEST 3 -- I-TEST ---------------------------- isum = 0.0 do 100 i=1,n isum = isum + (xdata(i)-xmed)**2 100 continue istat = isum/((n-1)*xsbiwt*xsbiwt) test(7) = istat C--- use an approximate formula for the percentage points ntilda = alog10(float(n-1)) i90 = 0.6376-1.1535*ntilda+0.1266*ntilda*ntilda if (n.lt.50) then i95 = 1.9065-2.5465*ntilda+0.5652*ntilda*ntilda else i95 = 0.7824-1.1021*ntilda+0.1021*ntilda*ntilda endif test(8) = 10**i90+0.982 test(9) = 10**i95+0.982 c---------------------- TEST 4 -- DIP TEST ---------------------------- if(n.ge.7) CALL DIPTST (XDATA,N,DIP,XL,XU,IFAULT,GCM,LCM,MN,MJ) test(10) = dip test(11) = xl test(12) = xu c---------------------- TEST 5 -- MODIFIED KS-TEST-------------------- c c-- Use the algorithm of Davis and Stephens c ITEST=2 ! FOR NORMAL DISTRIBUTION CALL EDFGOF(N,XDATA,ITEST,XBAR,S2,D,V,W2,W2P,U2,U2P, + A2,A2P,IFAULT) test(13) = d if(d.le..775) test(14)=.25 if(d.gt..775.and.d.le..819) test(14)=.15 if(d.gt..819.and.d.le..895) test(14)=.10 if(d.gt..895.and.d.le..995) test(14)=.05 if(d.gt..995.and.d.le.1.035) test(14)=.025 if(d.gt.1.035) test(14)=.01 c---------------------- TEST 6 -- MODIFIED KUIPER STATISTIC --------- test(15) = v if(v.le.1.32) test(16)=.25 if(v.gt.1.32.and.v.le.1.386) test(16)=.15 if(v.gt.1.386.and.v.le.1.489) test(16)=.10 if(v.gt.1.489.and.v.le.1.585) test(16)=.05 if(v.gt.1.585.and.v.le.1.693) test(16)=.025 if(v.gt.1.693) test(16)=.01 c-------------- TEST 7 -- MODIFIED CRAMER-VON-MISES STATISTIC ------- test(17) = W2 test(18) = W2P ! probability c----------------------- TEST 8 -- MODIFIED WATSON STATISTIC --------- test(19) = U2 test(20) = U2P ! probability c--------------- TEST 9 -- MODIFIED ANDERSON-DARLING STATISTIC-------- test(21) = A2 test(22) = A2P ! probability return end c----------------------------------------------------------------------------- SUBROUTINE NSCOR2(S,N,N2,IFAULT) c----------------------------------------------------------------------------- c--- Algorithm AS 177.3 Appl. Stat. (1982) Vol. 31, No.2 c c--- Approximation for rankits. c real s(n2),eps(4),dl1(4),dl2(4),gam(4),lam(4),bb,d, + b1,an,ai,e1,e2,l1,correc,ppnd data eps(1),eps(2),eps(3),eps(4) + /.419885e0,.450536e0,.456936e0,.468488e0/, + dl1(1),dl1(2),dl1(3),dl1(4) + /.112063e0,.121770e0,.239299e0,.215159e0/, + dl2(1),dl2(2),dl2(3),dl2(4) + /.080122e0,.111348e0,-.211867e0,-.115049e0/, + gam(1),gam(2),gam(3),gam(4) + /.474798e0,.469051e0,.208597e0,.259784e0/, + lam(1),lam(2),lam(3),lam(4) + /.282765e0,.304856e0,.407708e0,.414093e0/, + bb/-.283833e0/,d/-.106136e0/,b1/.5641896e0/ ifault = 3 if (n2 .ne. n/2) return ifault = 1 if (n .le. 1) return ifault = 0 if (n .gt. 2000) ifault = 2 s(1) = b1 if (n .eq. 2) return c c--- calculate normal areas for 3 largest rankits c an = n k = 3 if (n2 .lt. k) k = n2 do 5 i = 1,k ai = i e1 = (ai - eps(i))/(an + gam(i)) e2 = e1**lam(i) s(i) = e1 + e2*(dl1(i) + e2*dl2(i))/ an - correc(i,n) 5 continue if (n2 .eq. k) goto 20 c c--- calculate normal areas for remaining rankits c do 10 i = 4,n2 ai = i l1 = lam(4) + bb/(ai + d) e1 = (ai - eps(4))/ (an + gam(4)) e2 = e1**l1 s(i) = e1 + e2*(dl1(4) + e2*dl2(4))/ an - correc(i,n) 10 continue c c--- convert normal tail areas to normal deviates c 20 do 30 i = 1,n2 30 s(i) = -ppnd(s(i),ifault) return end c----------------------------------------------------------------------------- SUBROUTINE PRINTER (OUNIT,FILENAME,R) c----------------------------------------------------------------------------- c--- This routine writes out the results of the statistical routine c ROSTAT in a format which is comprehensible by the average human c****************************************************************************** parameter (iii=10000) integer ounit real*4 r(iii) character*32 filename common/print/ jtnor,jtgap,jtcon,jtchi,jtasy,jtcorr,jtot 3 format(/) 2 format(//) 1 format('1') c--- Declare the number of decimals you want to output if(abs(r(2)).gt.50.) id=1 if(abs(r(2)).le.50.) id=3 c--- write out sorted data write (ounit,1) write (ounit,5) filename 5 format(t25,'FILE = ',a) write (ounit,2) n = ifix(r(1)) write (ounit,6) n 6 format(t25,'SORTED DATA, N = ',i3) write (ounit,7) 7 format (t25,'====================') write (ounit,3) write (ounit,10) (r(j),j=2,n+1) 10 format(5(1x,f10.)) id=3 c--- letter value matrix write (ounit,2) write (ounit,2) write (ounit,11) 11 format(t25,'LETTER VALUE MATRIX') write (ounit,111) 111 format(t25,'===================') write (ounit,2) write (ounit,12) 12 format(t10,'DEPTH',T20,'LOWER',T30,'UPPER',T40,'MID',T50, + 'SPREAD',T60,'RATIO',t70,'PSEUDO') write (ounit,3) nl = ifix(r(n+2)) do 14 i=1,nl k = n+3+(i-1)*7 write (ounit,13) i,r(k),r(k+1),r(k+2),r(k+3),r(k+4), + r(k+5),r(k+6) 14 continue 13 format(1x,i2,t5,f10.,t15,f10.,t25,f10.,t35,f10., + t45,f10.,t55,f10.,t65,f10.) c--- classical estimators write (ounit,1) write (ounit,5) filename write (ounit,2) write (ounit,15) 15 format(t25,'CLASSICAL ESTIMATORS') write (ounit,115) 115 format(t25,'====================') write (ounit,2) k=n+7*nl write (ounit,16) 16 format(t10,'MEAN',t20,'SDEV',T30,'ADEV',T40,'SKEW',T50, + 'KURT',T60,'VAR') WRITE (OUNIT,3) write (ounit,17) r(k+3),r(k+5),r(k+4),r(k+7),r(k+8),r(k+6) 17 format(t5,f10.,t15,f10.,t25,f10.,t35,f10.,t45,f10., + t57,E10.4) c--- location estimates write (ounit,2) write (ounit,20) 20 format(t25,'ESTIMATES OF CENTRAL LOCATION') write (ounit,201) 201 format(t25,'=============================') write (ounit,2) write (ounit,21) 21 format (t10,'MEAN',t20,'3S-MEAN',T30,'N-3SIG',T40,'KS-MEAN',T50, + 'N-KS ',T60,'MEDIAN') write (ounit,3) write(ounit,22) r(k+9),r(k+11),ifix(r(k+10)),r(k+13), + ifix(r(k+12)),r(k+20) 22 format (t5,f10.,t15,f10.,t25,i10,t35,f10., + t45,i10,t55,f10.) write (ounit,2) write (ounit,211) 211 format(T8,'TRIM-5',t18,'TRIM-10',T28,'TRIM-20',T38,'MID-MEAN', + T48,'TRIM-30',T58,'TRIM-40') write (ounit,3) write(ounit,221) r(k+14),r(k+15),r(k+16),r(k+17),r(k+18),r(k+19) 221 format (t5,f10.,t15,f10.,t25,f10.,t35,f10., + t45,f10.,t55,f10.) write (ounit,2) write (ounit,23) 23 format (t10,'BMED',T20,'TRIMEAN',T30,'BIWT',T40,'BIWT-1') write (ounit,3) write (ounit,24) r(k+21),r(k+22),r(k+23),r(k+24) 24 format (t5,f10.,t15,f10.,t25,f10.,t35,f10.) if (r(k+60).eq.0.) goto 2000 write (ounit,2) write (ounit,*) ' The above central location estimates are ' write (ounit,*) ' with respect to the biweight estimate = ', + r(k+59) 2000 continue c--- scale estimates write (ounit,1) write (ounit,5) filename write (ounit,2) write (ounit,25) 25 format(t25,'ESTIMATES OF DATA SCALE') write (ounit,251) 251 format(t25,'=======================') write (ounit,2) write (ounit,26) 26 format(t10,'SDEV',T20,'3S-DEV',T30, 'N-3SIG',T40,'KS-DEV', + T50,'N-KS') write (ounit,3) write (ounit,27) r(k+30),r(k+32),ifix(r(k+31)), + r(k+34),ifix(r(k+33)) 27 format(t5,f10.,t15,f10.,t25,i10,t35,f10.,t45,i10) write (ounit,2) write (ounit,28) 28 format(T10,'SIG-F',T20,'SIG-MAD',t30,'SIG-GAP',t40,'S-BIWT', + T50,'S-BIWT1') write (ounit,3) write (ounit,29) r(k+38),r(k+40),r(k+58),r(k+41),r(k+42) 29 format (t5,f10.,t15,f10.,t25,f10.,t35, + f10.,t45,f10.) c--- confidence intervals write (ounit,1) write (ounit,5) filename write (ounit,3) if (r(jtcon+165).eq.1) write (ounit,50) if (r(jtcon+165).eq.0) write (ounit,502) 50 format(t5,'ERRORS IN LOCATION: FORMULAE + (XMSC IS 3-SIGMA CLIPPED)') 502 format(t5,'ERRORS IN LOCATION: FORMULAE + ( XMSC IS NOT CLIPPED )') write (ounit,501) 501 format(t5,'========================================== + ==============') write (ounit,3) write (ounit,51) 51 format(t20,'LOCATION',T30,'PLU/MIN',T40,'LOWER',T50,'UPPER') write(ounit,52) 52 format(t5,'XMSC') if (r(jtcon+165).eq.1.) k=n+7*nl+11 ! 3-sigma mean if (r(jtcon+165).eq.0.) k=n+7*nl+9 ! mean do 55 i=1,4 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 55 continue 53 format(t15,f10.,t25,f10.,t35,f10.,t45,f10.) write(ounit,54) 54 format(t5,'XMDC') k=n+7*nl+20 do 56 i=5,8 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 56 continue write(ounit,57) 57 format(t5,'XTSC') k=n+7*nl+23 do 59 i=9,12 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 59 continue write(ounit,3) write(ounit,80) 80 format(t20,'ERRORS IN LOCATION: BOOTSTRAP') write(ounit,801) 801 format(t20,'==============================') write (ounit,3) write(ounit,81) 81 format(t20,'LOCATION',T30,'LENGTH',T40,'LOWER',T50,'UPPER') write(ounit,82) 82 format(t5,'MEAN-1') k=n+7*nl+3 do 85 i=37,44,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 85 continue write(ounit,86) 86 format(t5,'MEAN-2') k=n+7*nl+3 do 87 i=45,52,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 87 continue write(ounit,88) 88 format(t5,'MEAN-3') k=n+7*nl+3 do 89 i=53,60,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 89 continue write(ounit,90) 90 format(t5,'MEAN-4') k=n+7*nl+3 do 91 i=133,140,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 91 continue write(ounit,810) 810 format(t5,'BIWT-1') k=n+7*nl+23 do 811 i=61,68,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 811 continue write(ounit,812) 812 format(t5,'BIWT-2') k=n+7*nl+23 do 813 i=69,76,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 813 continue write(ounit,814) 814 format(t5,'BIWT-3') k=n+7*nl+23 do 815 i=77,84,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 815 continue write(ounit,92) 92 format(t5,'BIWT-4') k=n+7*nl+23 do 93 i=141,148,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 93 continue write (ounit,1) write (ounit,5) filename write (ounit,3) if(r(jtcon+165).eq.1 ) write (ounit,65) if(r(jtcon+165).eq.0 ) write (ounit,655) 65 format(t10,'ERRORS IN SCALE: CLASSICAL (3-SIGMA CLIPPED)') 655 format(t10,'ERRORS IN SCALE: CLASSICAL ( NO CLIPPING )') write (ounit,651) 651 format(t10,'=============================================') write (ounit,3) write (ounit,66) 66 format(t15,'SCALE',T25,'LO-ERROR',T35,'HI-ERROR',T45,'LOWER', + T55,'UPPER') write (ounit,3) write (ounit,642) 642 format(t5,'XXSC') if (r(jtcon+165).eq.1.) k=n+7*nl+32 ! 3-sigma standard deviation if (r(jtcon+165).eq.0.) k=n+7*nl+30 ! standard deviation do 67 i=1,4 l=jtchi+1+(i-1)*2 write (ounit,68) r(k),r(k)-r(l),r(l+1)-r(k),r(l),r(l+1) 67 continue 68 format(t10,f10.,t20,f10.,t30,f10.,t40, + f10.,t50,f10.) c write (ounit,641) c641 format(t5,'XJLC') c k=n+7*nl+45 c do 64 i=17,20 c l=jtcon+i c rv = 10**r(k) c rlv = 10**(r(k)-r(l)) c ruv = 10**(r(k)+r(l)) c rle = rv-rlv c rue = ruv-rv c write(ounit,68) rv,rle,rue,rlv,ruv c64 continue c write (ounit,643) c643 format(t5,'XJBLC') c k=n+7*nl+49 c do 644 i=25,28 c l=jtcon+i c rv = 10**r(k) c rlv = 10**(r(k)-r(l)) c ruv = 10**(r(k)+r(l)) c rle = rv-rlv c rue = ruv-rv c write(ounit,68) rv,rle,rue,rlv,ruv c644 continue c write (ounit,645) c645 format(t5,'XJGLC') c k=n+7*nl+53 c do 646 i=33,36 c l=jtcon+i c rv = 10**r(k) c rlv = 10**(r(k)-r(l)) c ruv = 10**(r(k)+r(l)) c rle = rv-rlv c rue = ruv-rv c write(ounit,68) rv,rle,rue,rlv,ruv c646 continue write (ounit,3) write (ounit,667) 667 format(t20,'ERRORS IN SCALE: JACKNIFE') write (ounit,668) 668 format(t20,'==========================') write(ounit,3) write (ounit,511) 511 format(t20,'SCALE',T30,'PLU/MIN',T40,'LOWER',T50,'UPPER') write (ounit,3) write(ounit,60) 60 format(t5,'XJKC') k=n+7*nl+30 do 62 i=13,16 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 62 continue write(ounit,601) 601 format(t5,'XJBC') k=n+7*nl+41 do 621 i=21,24 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 621 continue write(ounit,602) 602 format(t5,'XJGC') k=n+7*nl+58 do 622 i=29,32 l=jtcon+i write(ounit,53) r(k),r(l),r(k)-r(l),r(k)+r(l) 622 continue c--- chi-square and t-values and optional corrections write (ounit,1) write (ounit,5) filename write (ounit,2) write (ounit,70),n 70 format(t5,'CHI-SQUARE AND T-VALUES FOR N = ',i3) write (ounit,701) 701 format(t5,'===================================') write (ounit,2) write (ounit,71) 71 format(t10,'CHI-LOW',T20,'CHI-HIGH') write (ounit,3) do 72 i=1,4 l=jtchi+9+(i-1)*2 write (ounit,73) r(l),r(l+1) 72 continue 73 format(t5,f10.,t15,f10.) write (ounit,2) write (ounit,2) write (ounit,74) 74 format(t5,'T-VALUES') write (ounit,741) 741 format(t5,'========') write (ounit,3) write (ounit,75) r(jtchi+17),r(jtchi+18),r(jtchi+19),r(jtchi+20) 75 format (t5,'t(n-1)',t17,f10.,t27,f10.,t37,f10.,t47, + f10.) write (ounit,3) write (ounit,76) r(jtchi+21),r(jtchi+22),r(jtchi+23),r(jtchi+24) 76 format (t5,'t(n-1)/1.075',t17,f10.,t27,f10.,t37,f10., + t47,f10.) write (ounit,3) write (ounit,77) r(jtchi+25),r(jtchi+26),r(jtchi+27),r(jtchi+28) 77 format (t5,'t(0.7*(n-1))',t17,f10.,t27,f10.,t37,f10., + t47,f10.) write (ounit,2) write (ounit,2) write (ounit,774) 774 format(t5,'OPTIONAL COSMOLOGICAL CORRECTIONS') write (ounit,7741) 7741 format(t5,'=================================') write (ounit,3) write (ounit,78) r(jtcorr+1),r(jtcorr+2),r(jtcorr+3) 78 format(' LOCATION CONFIDENCE CORRECTION + // (add in quadrature): ',f10., + // ' SCALE CORRECTION (subtract in quadrature): ',f10., + // ' SCALE CONFIDENCE CORRECTION (add in quadrature): ', + // f10.) write(ounit,1) write(ounit,5) filename write(ounit,3) write(ounit,816) 816 format(t20,'ERRORS IN SCALE: BOOTSTRAP') write(ounit,817) 817 format(t20,'===========================') write(ounit,3) write(ounit,818) 818 format(t20,'SCALE',T30,'LENGTH',T40,'LOWER',T50,'UPPER') write(ounit,819) 819 format(t5,'SDEV-1') k=n+7*nl+30 do 820 i=85,92,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 820 continue write(ounit,821) 821 format(t5,'SDEV-2') k=n+7*nl+30 do 822 i=93,100,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 822 continue write(ounit,823) 823 format(t5,'SDEV-3') k=n+7*nl+30 do 824 i=101,108,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 824 continue write(ounit,831) 831 format(t5,'SDEV-4') k=n+7*nl+30 do 832 i=149,156,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 832 continue write(ounit,825) 825 format(t5,'S-BWT-1') k=n+7*nl+41 do 826 i=109,116,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 826 continue write(ounit,827) 827 format(t5,'S-BWT-2') do 828 i=117,124,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 828 continue write(ounit,829) 829 format(t5,'S-BWT-3') do 830 i=125,132,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 830 continue write(ounit,833) 833 format(t5,'S-BWT-4') do 834 i=157,164,2 l=jtcon+i write(ounit,53) r(k),r(l+1)-r(l),r(l),r(l+1) 834 continue c--- tails write (ounit,1) write (ounit,5) filename write (ounit,2) write (ounit,30) 30 format(t25,'TAILS OF THE DATA SET') write (ounit,301) 301 format(t25,'=====================') write (ounit,2) write (ounit,311) r(jtasy+1),r(jtasy+2) 311 format(t5,'ASYMMETRY INDEX FOR DATA SET: ',f6.3,// + t5,'TRANSFORMATION POWER: ',f6.3,//) k = n+7*nl write (ounit,31) r(k+56),r(k+57) 31 format (t5,'RAW TAIL INDEX FOR DATA SET: ',f5.3,// + t5,'SCALED TAIL INDEX FOR DATA SET: ',f5.3,//) write (ounit,32) 32 format(t5'COMPARISON INDICES FOR OTHER DISTRIBUTIONS',// + ,t5,'UNIFORM = 0.84'/ + ,t5,'TRIANGULAR = 0.91'/ + ,t5,'GAUSSIAN = 1.00'/ + ,T5,'CN(.05,3) = 1.02'/ + ,T5,'CN(.05,10) = 1.03'/ + ,T5,'LOGISTIC = 1.05'/ + ,T5,'CN(.20,3) = 1.08'/ + ,T5,'DBL EXP = 1.23'/ + ,T5,'CN(.20,10) = 1.26'/ + ,T5,'SLASH = 1.43'/ + ,T5,'CAUCHY = 1.62'////) c--- gaps write (ounit,33) 33 format(t25,'WEIGHTED GAPS IN THE DATA') write (ounit,331) 331 format(t25,'=========================') write (ounit,2) ngap = ifix(r(jtgap+1)) write (ounit,34) (r(k),k=jtgap+2,jtgap+1+ngap) 34 format(5(1x,f7.3)) nbig = ifix(r(jtgap+2+ngap)) write (ounit,2) write (ounit,35) nbig 35 format(t5,'THERE ARE ',i3,' BIG GAPS IN THE DATA SET') write (ounit,2) write (ounit,36) 36 format(t5,'IBIG',T18,'ZBIG',T30,'ZPROB',T40,'TPROB') write (ounit,3) do 37 i=1,nbig k=jtgap+3+ngap+(i-1)*4 write (ounit,38) ifix(r(k)),r(k+1),r(k+2),r(k+3) 37 continue 38 format (t5,i3,t15,f7.3,t25,f10.4,t35,f10.4) c--- normality test results write (ounit,1) write (ounit,5) filename write (ounit,2) write (ounit,40) 40 format (t25,'RESULTS OF THE NORMALITY TESTS') write (ounit,401) 401 format (t25,'==============================') write (ounit,2) write (ounit,41) r(jtnor+1) write (ounit,42) r(jtnor+2) write (ounit,43) r(jtnor+3),r(jtnor+4) write (ounit,44) r(jtnor+5),r(jtnor+23) write (ounit,45) r(jtnor+6),r(jtnor+24) write (ounit,1053) r(jtnor+25),r(jtnor+26) write (ounit,46) r(jtnor+7),r(jtnor+8),r(jtnor+9) write (ounit,47) r(jtnor+10),r(jtnor+11),r(jtnor+12) write (ounit,48) r(jtnor+13),r(jtnor+14) write (ounit,49) r(jtnor+15),r(jtnor+16) write (ounit,1050) r(jtnor+17),r(jtnor+18) write (ounit,1051) r(jtnor+19),r(jtnor+20) write (ounit,1052) r(jtnor+21),r(jtnor+22) 41 format(T3,'A-TEST STATISTIC: ',F8.3) 42 FORMAT(T3,'U-TEST STATISTIC: ',F8.3) 43 FORMAT(T3,'W-TEST STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 44 FORMAT(T3,'B1-TEST STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 45 FORMAT(T3,'B2-TEST STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 1053 FORMAT(T3,'B1 & B2 OMNI TEST : ',F8.3,T37, + 'PROBABILITY: ',F8.3) 46 FORMAT(T3,'I-TEST STATISTIC: ',F8.3,T37, + '90% POINT: ',F8.3, + T59,'95% POINT: ',F8.3) 47 FORMAT(T3,'DIP STATISTIC: ',F8.3,T37, + 'LOWER MODAL: ', + F8.,T59,'UPPER MODAL:',F8.) 48 FORMAT(T3,'KS STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 49 FORMAT(T3,'KUIPER V-STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 1050 FORMAT(T3,'CRAMER VON-MISES W2-STAT:',F8.3,T37, + 'PROBABILITY: ',F8.3) 1051 FORMAT(T3,'WATSON U2-STATISTIC: ',F8.3,T37, + 'PROBABILITY: ',F8.3) 1052 FORMAT(T3,'ANDERSON-DARLING A2-STAT:',F8.3,T37, + 'PROBABILITY: ',F8.3) return end c------------------------------------------------------------------------------ SUBROUTINE SIGMA (DATA,N,SXAVE,SXSDEV,NTOT) c------------------------------------------------------------------------------ c--- This routine performs a "3-sigma" clipping of the data using the c cleaning technique described in Yahil and Vidal (1977). The procedure c calls for deleting the data point most deviant from the sample mean, c redetermining the sample mean and standard deviation. If the removed c point deviates by more than 3 standard deviations from the new mean c it is rejected as a contamination. Each point in the list is tested c in order of their apparent deviation, until no further rejection occurs. c The list is now considered "clean," and the final estimates of mean and c standard deviation are accepted. c c--- Note that this routine will exit to the main program if the input number c of data point N <= 5, or if clipped to NTOT=5 c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) dimension data(n),tdata(iii),sdata(iii) if(n.le.5) then sxave=-999 sxsdev=-999 ntot=n return endif c--- copy the original data set do 5 i=1,n sdata(i) = data(i) 5 continue c--- find intial estimate of sample mean and deviation ave = PARAM(SDATA,N,1) ! find new mean sdev = PARAM(SDATA,N,2) ! find new sdev nout=0 c--- create a list of data ordered by their distance from the sample mean 20 continue if(n-nout.eq.5) then ! get out of here sxave=-999 sxsdev=-999 ntot=5 return endif do 10 i=1,n-nout tdata(i) = abs(sdata(i)-ave) 10 continue CALL SORT2 (N-NOUT,TDATA,SDATA) ! sorts TDATA hold SDATA c--- drop first deviant point temp = sdata(n-nout) sdata(n-nout) = -32368. ave = PARAM(SDATA,N-NOUT,1) ! find new mean sdev = PARAM(SDATA,N-NOUT,2) ! find new sdev c--- test if most deviant point is greater than 3-sigma from sample mean c calculated without including it in the estimate, c if so, then delete it from consideration if(abs(temp-ave) .gt. 3*sdev) then ! delete from list nout = nout + 1 else sdata(n-nout)=temp ave = PARAM(SDATA,N-NOUT,1) ! calculate mean including temp sdev = PARAM(SDATA,N-NOUT,2) ! calculate sdev including temp sxave=ave sxsdev=sdev ntot=n-nout return endif goto 20 end c------------------------------------------------------------------------------ SUBROUTINE SORT(N,RA) c------------------------------------------------------------------------------ c--- Routine to do a heapsort of a data array RA c Stolen (unabashedly) from NUMERICAL RECIPES c c**************************************************************************** implicit real*4 (a-h,o-z) dimension ra(n) l=n/2+1 ir=n 10 continue if(l.gt.1)then l=l-1 rra=ra(l) else rra=ra(ir) ra(ir)=ra(1) ir=ir-1 if(ir.eq.1)then ra(1)=rra return endif endif i=l j=l+l 20 if(j.le.ir)then if(j.lt.ir)then if(ra(j).lt.ra(j+1))j=j+1 endif if(rra.lt.ra(j))then ra(i)=ra(j) i=j j=j+j else j=ir+1 endif go to 20 endif ra(i)=rra go to 10 end c------------------------------------------------------------------------------ SUBROUTINE SORT2(N,RA,RB) c------------------------------------------------------------------------------ c--- Routine to do a heapsort of a data array RA and RB c Stolen (unabashedly) from NUMERICAL RECIPES c c**************************************************************************** implicit real*4 (a-h,o-z) dimension ra(n),rb(n) l=n/2+1 ir=n 10 continue if(l.gt.1)then l=l-1 rra=ra(l) rrb=rb(l) else rra=ra(ir) rrb=rb(ir) ra(ir)=ra(1) rb(ir)=rb(1) ir=ir-1 if(ir.eq.1)then ra(1)=rra rb(1)=rrb return endif endif i=l j=l+l 20 if(j.le.ir)then if(j.lt.ir)then if(ra(j).lt.ra(j+1))j=j+1 endif if(rra.lt.ra(j))then ra(i)=ra(j) rb(i)=rb(j) i=j j=j+j else j=ir+1 endif go to 20 endif ra(i)=rra rb(i)=rrb go to 10 end c------------------------------------------------------------------------------ SUBROUTINE STATS (N, Z, D, V, W2, U2, A2, IFAULT) c------------------------------------------------------------------------------ C ALGORITHM AS 248.2 APPL.STATIST. (1989), VOL. 38, NO. 3) C C Computes the goodness-of-fit statistics D, V, W**2, U**2 C and A**2 from the transformed Z values parameter (iii = 10000) integer i, ifault, n, ni real*4 z(n), d, v, w2, u2, a2, ri, rn, dplus, dminus, d1, d2, + sumz, twon, zm1, a2sum, zero, small, half, one, two, twelve C C Initialize constants C data zero / 0.0 / , small / 1.0e-37 / , half / 0.5 /, + one / 1.0 /, two / 2.0 /, twelve / 12.0 / rn = n twon = two*rn C C Calculating the Kologorov statistics DPLUS, DMINUS, D C and the Kuiper statistic V. C dplus = zero dminus = zero ifault = 4 do 10 i = 1, n if (z(i) .le. zero .or. z(i) .ge. one) return ri = i d1 = ri / rn - z(i) if (d1 .gt. dplus) dplus = d1 d2 = z(i) - (ri - one) / rn if (d2 .gt. dminus) dminus = d2 10 continue ifault = 0 d = dplus if (dminus .gt. dplus) d = dminus v = dplus + dminus C C Calculating the Cramer-Von Mises statistic W2 and the C Watson statistic U2. C w2 = one / (twelve * rn) sumz = zero do 20 i = 1,n w2 = w2 + (z(i) - (two * i - one) / twon) ** 2 sumz = sumz + z(i) 20 continue u2 = w2 - rn * (sumz / rn - half) ** 2 c C Calculating the Anderson-Darling statistic A2 C a2sum = zero ni = n do 30 i = 1, n if (z(i) .lt. small) z(i) = small zm1 = one - z(ni) if (zm1 .lt. small) zm1 = small a2sum = a2sum + (two * i - one) * (log(z(i)) + log(zm1)) ni = ni - 1 30 continue a2 = (-a2sum) / rn - rn return end c------------------------------------------------------------------------------ SUBROUTINE SYMMETRY (XDATA,N,TEE) c------------------------------------------------------------------------------ c-- for a reference see Finch,S."Robust Univariate Test of Symmetry",J. of c American Stat. Assoc., June,1977,Vol.72 No.358,p.387. c This compares the gaps in the ordered data set from the right side to c the left side of the median as a test of the symmetry. c c****************************************************************************** parameter (iii = 10000) real num,den,inv,yi,yn,vi,prod,tee,wii real xdata(n),wi(iii) num=0 den=0 c-- calculate the weights for each data point do 10 i=1,n/2 wii=0. do 20 j=i,n-i inv=1./float(j) wii=wii+inv 20 continue wi(i)=wii 10 continue c-- get the gaps in the data and calculte the statistic do 30 i=1,n/2 yi = xdata(i+1) - xdata(i) yn = xdata(n-i+1) - xdata(n-i) if (yi+yn.ne.0) vi = (yn-yi)/(yi+yn) prod = wi(i)*vi num = num + prod c-- have to normalize by the sum of the sqaure of the weights den = den + (wi(i)*wi(i)) 30 continue tee = num/(den**0.5) return end c------------------------------------------------------------------------------ SUBROUTINE TAIL (XDATA,N,TINDEX1,TINDEX2) c------------------------------------------------------------------------------ c--- TAIL provides a robust estimate of the weight of the tails c of a symmetric distribution (UREDA 322). Essentially, c this means we calculate: c c Q(.90) - Q(.10) c TINDEX1 = --------------- c Q(.75) - Q(.25) c c where Q is simply the rank statistic closest to the ideal c depth DEPTHS. c c TINDEX1 can then be compared to similar values as calculated c for a variety of interesting distributions. c c TINDEX2 is the scaled version of TINDEX1 relative to a gaussian c c--- NOTE -- TAIL assumes that XDATA is sorted low to high on input c c**************************************************************************** implicit real*4 (a-h,o-z) dimension xdata(n) data dt/0.33333/ c--- obtain estimate of Q parameter d90 = (n+dt)*.90+ dt d10 = (n+dt)*.10+ dt d75 = (n+dt)*.75+ dt d25 = (n+dt)*.25+ dt ig90 = int(d90) ig10 = int(d10) ig75 = int(d75) ig25 = int(d25) r90 = d90 - ig90 r10 = d10 - ig10 r75 = d75 - ig75 r25 = d25 - ig25 if (n .le. 7) then tindex1 = 0.0 tindex2 = 0.0 goto 10 endif q90 = (1.-r90)*xdata(ig90)+r90*xdata(ig90+1) q10 = (1.-r10)*xdata(ig10)+r10*xdata(ig10+1) q75 = (1.-r75)*xdata(ig75)+r75*xdata(ig75+1) q25 = (1.-r25)*xdata(ig25)+r25*xdata(ig25+1) c--- obtain estimate of TINDEX tindex1 = (q90-q10)/(q75-q25) tindex2 = tindex1/1.90 ! scaled to gaussian 10 return end c------------------------------------------------------------------------------ SUBROUTINE TRANS(N,NLV,XLV,POW) c------------------------------------------------------------------------------ c-- Trans gives the power for the data set to make a more symmetric c distribution (for a discussion see UREDA pp.105-111). c c--- CODE UNDER REVIEW ! USE WITH CAUTION. c********************************************************************** real*4 x,y,pp(8),pow,xm real*4 xlv(15,2),xxlv(15,2),xdepths(15) real*4 mmlv(8),sslv(15),rrlv(8),cclv(8) integer nlv,n,nnlv,nlvv c-- from the midsummaries a value pp will be found for each and c the power trans will be the median of those values xm = xlv(1,1) do 10 i=2,nlv x=((xlv(i,1)-xm)**2 + (xlv(i,2)-xm)**2)/(4.*xm) y=((xlv(i,1)+xlv(i,2))/2.)-xm pp(i-1)=1.-(y/x) 10 continue nlvv=nlv-1 c-- need to get the median of the pp(i) CALL MDIAN1(PP,NLVV,XMED) pow=xmed return end c------------------------------------------------------------------------------ SUBROUTINE TRIM (XDATA,N,ALPHA,TRIMMED) c------------------------------------------------------------------------------ c--- The Trimmed mean finds the mean of the set of N ordered statistics c after trimming the percentage ALPHA from both ends of the set. c The value TRIMMED is returned as the mean of the trimmed set c of ordered statistics. The weights given to the ordered stats c and the formula for computation of the trimmed mean can be c found on page 311 of UREDA. c c**************************************************************************** implicit real*4 (a-h,o-z) dimension xdata(n) data d1,d2,n1,n2,zero/1.0,2.0,1,2,0.0/ ig = int(alpha * n) r = (alpha*n) - float(ig) sum1 = zero do 11 i = ig+n2,n-ig-n1 sum1 = sum1 + xdata(i) 11 continue sum3 = (d1-r)*(xdata(ig+n1) + xdata(n-ig)) trimmed = (d1/(float(n)*(d1 -(d2*alpha))))*(sum3+sum1) return end c----------------------------------------------------------------------------- SUBROUTINE TRIMEAN (XDATA,N,XMED,XFL,XFU,TRIMN) c----------------------------------------------------------------------------- c--- The TRIMEAN is calculated for the set of N ORDERED statistics c passed in the array XDATA. The value of the TRIMEAN is returned c as TRIMN. The TRIMEAN is given by the formula c c TRIMEAN = .25(FL + 2M + FU) c c where FL,M,FU are the lower fourth, median, and upper fourth c respectively and are passed in the array XLETTER. For more c information see pages 313,314 in UREDA. c c***************************************************************************** implicit real*4 (a-h,o-z) dimension xdata(n) data do4,d2/0.25,2.0/ trimn = do4*(xfu + (d2*xmed) + xfl) return end c----------------------------------------------------------------------------- SUBROUTINE UPDATE (N,XDATA,M,XOUT) c----------------------------------------------------------------------------- real*4 xdata(n),xout(n) c--- only copy data from XDATA to XOUT if it is not = -32368 j=0 miss=0 do 10 i=1,n if (xdata(i).ne.-32368.) then j=j+1 xout(j)=xdata(i) else miss = miss + 1 endif 10 continue m = j return end c----------------------------------------------------------------------------- SUBROUTINE XAD (XDATA,N,XMED,XADM) c----------------------------------------------------------------------------- c--- The XAD subroutine calculates the Mean Absolute Deviation from c the sample median. The median, M , is subtracted from each c ORDERED statistic and then the absolute value is taken. c The AD is then defined to be the mean of these differences. c new set of statistics and is returned as XADM. The AD can c be defined: c c ADM = mean{ abs(x(i) - M) } c c where the x(i) are the values passed in the array XDATA, and c the median, M, is passed in the array XLETTER. c For more information see page 408 in UREDA. c c c**************************************************************************** implicit real*4 (a-h,o-z) dimension xdata(n) data zero/0.0/ xsum= zero do 11 i = 1,n xsum = abs(xdata(i) - xmed) + xsum 11 continue xadm = xsum/n return end c----------------------------------------------------------------------------- SUBROUTINE XBIWT (XDATA,N,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1) c----------------------------------------------------------------------------- c--- The subroutine XBIWT provides an estimator of the location and c scale of the data set XDATA. The scale uses the Biweight function c in the general formula of "A-estimators." This formula is given c on page of 416 in UREDA (formula 4). The BIWEIGHT scale estimate c is returned as the value XSBIWT. The BIWEIGHT function is given c by: c c u((1-u*u)**2) abs(u) <= 1 c f(u) = c 0 abs(u) > 1 c c where u is defined by c c u = (XDATA(I) - M) / c*MAD . c c M, MAD, and c are the median, the median absolute deviation from c the median, and the tuning constant respectively. The tuning c constant is a parameter which is chosen depending on the sample c size and the specific function being used for the scale estimate. c (See page 417 in UREDA). Here we take c = 9.0. c c--- The biweght location is found using the formula: c c T = M + (sums) c c where M is the sample median and sums are c as given on page 421 in UREDA c c the tuning constant c is set to 6.0 for calculation c of the location as reccommended by Tukey () c c--- NOTE that the biweight is meant to be an iterated estimator, but one c commonly only takes the first step of the iteration. Here we report c both the one-step estimators (XLBIWT1, XSBIWT1) and the preferred c fully iterated versions (XLBIWT, XSBIWT). c c**************************************************************************** parameter (iii = 10000) implicit real*4 (a-h,o-z) dimension xdata(n),u1(iii),u2(iii), + xlb(11),xsb(11) real depths(15),xletter(15,2),mletter(8),sletter(15) real rletter(8),cletter(8) data zero,d6,d1,d5,d9/0.0,6.0,1.0,5.0,9.0/ c--- sort the data and find the median CALL MDIAN1(XDATA,N,XM) c--- call xmad to find the median absolute deviation CALL XMAD(XDATA,N,XM,XMADM) c--- must choose value of the tuning constant "c" c here c = 6.0 for the location estimator and c 9.0 for the scale estimator c1 = d6 c2 = d9 if (xmadm.le..0001) then xlbiwt=xm xlbiwt1=xm xsbiwt=xmadm xsbiwt1=xmadm goto 20 endif do 11 i = 1,n u1(i) = (xdata(i) - xm)/(c1*xmadm) u2(i) = (xdata(i) - xm)/(c2*xmadm) 11 continue s1 = zero s2 = zero s3 = zero s4 = zero do 12 i = 1,n if (abs(u2(i)) .lt. d1) then s1 = s1+(((xdata(i)-xm)**2)*(d1-(u2(i)*u2(i)))**4) s2 = s2+((d1-u2(i)*u2(i))*(d1-(d5*u2(i)*u2(i)))) endif if (abs(u1(i)) .lt. d1) then s3 = s3+(xdata(i)-xm)*(d1-u1(i)*u1(i))**2 s4 = s4+(d1-u1(i)*u1(i))**2 endif 12 continue c--- here are the one-step estimators xlbiwt1 = xm+s3/s4 xsbiwt1 = float(n)/(float(n-1))**0.5*s1**0.5/abs(s2) c--- now obtain the fully-iterated versions c--- solve for new estimates of u1 and u2 xlb(1) = xlbiwt1 xsb(1) = xsbiwt1 do 15 j = 2,11 !assuming 10 iterations is sufficient xmm = xlb(j-1) do 13 i = 1,n u1(i) = (xdata(i) - xmm)/(c1*xmadm) u2(i) = (xdata(i) - xmm)/(c2*xmadm) 13 continue s1 = zero s2 = zero s3 = zero s4 = zero do 14 i = 1,n if (abs(u2(i)) .lt. d1) then s1 = s1+(((xdata(i)-xmm)**2)*(d1-(u2(i)*u2(i)))**4) s2 = s2+((d1-u2(i)*u2(i))*(d1-(d5*u2(i)*u2(i)))) endif if (abs(u1(i)) .lt. d1) then s3 = s3+(xdata(i)-xmm)*(d1-u1(i)*u1(i))**2 s4 = s4+(d1-u1(i)*u1(i))**2 endif 14 continue xlb(j) = xlb(j-1)+s3/s4 xsb(j) = float(n)/(float(n-1))**0.5*s1**0.5/abs(s2) 15 continue xlbiwt = xlb(11) xsbiwt = xsb(11) 20 continue return end c----------------------------------------------------------------------------- SUBROUTINE XMAD (XDATA,N,XMED,XMADM) c----------------------------------------------------------------------------- c--- The XMAD subroutine calculates the Median Absolute Deviation from c the sample median. The median, M , is subtracted from each c ORDERED statistic and then the absolute value is taken. This new c set of of statistics is then resorted so that they are ORDERED c statistics. The MAD is then defined to be the median of this c new set of statistics and is returned as XMADM. The MAD can c be defined: c c XMADM = median{ abs(x(i) - M) } c c where the x(i) are the values passed in the array XDATA, and c the median, M, is passed in the array XLETTER. The set of stats c in the brackets is assumed to be resorted. For more information c see page 408 in UREDA. c c c*******************************************************