home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 474 / ems.arc / PITMAN.FOR < prev    next >
Text File  |  1985-12-12  |  23KB  |  870 lines
  1. C        
  2. C                                 PC-PITMAN
  3. C                            Randomization Tests
  4. C                                Version 2.0
  5. C                             December 12, 1985
  6. C
  7. C
  8. C                              Gerard E. Dallal
  9. C
  10. C               USDA Human Nutrition Research Center on Aging
  11. C                            at Tufts University
  12. C                           711 Washington Street
  13. C                             Boston, MA  02111
  14. C
  15. C                                    and
  16. C
  17. C                    Tufts University School of Nutrition
  18. C                             132 Curtis Street
  19. C                             Medford, MA 02155
  20. C        
  21. C
  22. C
  23. C                                  NOTICE
  24. C
  25. C       Documentation and original code copyright 1985 by  Gerard  E.  
  26. C       Dallal.  Reproduction of material for non-commercial purposes 
  27. C       is   permitted,   without  charge,   provided  that  suitable 
  28. C       reference is made to PC-PITMAN and its author.  
  29. C
  30. C       Neither PC-PITMAN nor its documentation should be modified in 
  31. C       any way without permission from the author,  except for those 
  32. C       changes  that are  essential  to move  PC-PITMAN  to  another 
  33. C       computer.  
  34. C
  35. C       Please acknowledge  PC-PITMAN in any manuscript that uses its 
  36. C       calculations.  
  37. C
  38. C
  39.       DIMENSION H(36000), DATA(100), DATAI(100), DCOPY(100), 
  40.      *   DATAG(100)
  41.       CHARACTER*1 QUERY
  42.       DATA MAXTOT /36000/, MAXOBS /100/
  43.       DATA IIN /0/, IOUT /0/, IFIN/10/
  44. C
  45.       WRITE(IOUT,1)
  46.     1 FORMAT (//35X,'PC-PITMAN'/
  47.      *   30X,'Randomization Tests'/
  48.      *   34X,'Version 2.0'/31X,'December 12, 1985'//
  49.      *   32X,'Gerard E. Dallal'/
  50.      *   17X,'USDA Human Nutrition Research Center on Aging'/
  51.      *   30X,'at Tufts University'/29X,'711 Washington Street'/
  52.      *   31X,'Boston, MA  02111'///36X,'NOTICE'//9X,'Please '
  53.      *   'acknowledge PC-PITMAN in any manuscript that uses its',/
  54.      *   33X,'calculations.'/)
  55. C
  56.    10 CALL READ(IIN,IOUT,IFIN,DATA,DATAG,MAXOBS,NOBS,M,N,NGRP)
  57.       XNOBS = NOBS
  58. C
  59.   110 WRITE(IOUT,120)
  60.   120 FORMAT (' Multiply data by what integer power of 10?  '$)
  61.       READ(IIN,*,ERR=110) DEC
  62.       NDEC = DEC
  63.       IF(FLOAT(NDEC).NE.DEC) GOTO 110
  64. C
  65. C        1.2 ADDITION... ROUND RATHER THAN TRUNCATE
  66. C
  67.       DO 130 I = 1, NOBS
  68.   130 DCOPY(I) = AINT(DATA(I) * 10.0**NDEC + SIGN(0.5,DATA(I)))
  69. C
  70.       IF (NGRP.EQ.2) GOTO 200
  71. C
  72.       CALL PITMN1(DATAI,DCOPY,NOBS,H,MAXTOT,IOUT)
  73. C
  74. C        T-TEST
  75. C
  76.       IF (NOBS.LE.1) STOP 'Fewer than 2 observations.'
  77.       AVG = 0.0
  78.       DO 140 I = 1, NOBS
  79.   140 AVG = AVG + DATA(I)
  80.       AVG = AVG / XNOBS
  81.       SSQ = 0.0
  82.       DO 150 I = 1, NOBS
  83.   150 SSQ = SSQ + (DATA(I) - AVG)**2
  84.       IF (SSQ.GT.0.0) GOTO 180
  85. C
  86.       WRITE(IOUT,170)
  87.   170 FORMAT (' All data values equal.  No t-test.')
  88.       GOTO 300
  89. C
  90.   180 T = AVG / SQRT(SSQ / ((XNOBS - 1) * XNOBS))
  91.       NDF = NOBS - 1
  92.       PVALT = 2.0 * PROBST(SIGN(T,-1.0),NDF)
  93.       WRITE(IOUT,190) T, NDF, PVALT
  94.   190 FORMAT(' Student''s t-statistic = ',G11.4,' with',
  95.      *   I3,' d.f.  (P-value =',F7.4,')')
  96.       GOTO 300
  97. C
  98.   200 CALL PITMN2(DATAI,DCOPY,M,N,NOBS,H,MAXTOT,IOUT)
  99. C
  100. C        T-TEST
  101. C
  102.       IF (NOBS.LE.2 .OR. M.EQ.0 .OR. N.EQ.0) 
  103.      *   STOP 'Degenerate data set.'
  104.       AVG1 = 0.0
  105.       DO 230 I = 1, M
  106.   230 AVG1 = AVG1 + DATA(I)
  107.       AVG1 = AVG1 / FLOAT(M)
  108.       AVG2 = 0.0
  109.       DO 240 I = 1, N
  110.   240 AVG2 = AVG2 + DATA(M+I)
  111.       AVG2 = AVG2 / FLOAT(N)
  112.       SSQ = 0.0
  113.       DO 250 I = 1, M
  114.   250 SSQ = SSQ + (DATA(I) - AVG1)**2
  115.       DO 260 I = 1, N
  116.   260 SSQ = SSQ + (DATA(M+I) - AVG2)**2
  117.       IF (SSQ.GT.0.0) GOTO 280
  118. C
  119.       WRITE(IOUT,270)
  120.   270 FORMAT (' Within group variance zero.  No t-test.')
  121.       GOTO 300
  122. C
  123.   280 T = (AVG1 - AVG2) / SQRT((SSQ / FLOAT(NOBS-2)) *
  124.      *   (1.0 / FLOAT(M) + 1.0 / FLOAT(N)))
  125.       NDF = NOBS - 2
  126.       PVALT = 2.0 * PROBST(SIGN(T,-1.0),NDF)
  127.       WRITE(IOUT,190) T, NDF, PVALT
  128. C
  129.   300 WRITE (IOUT,310)
  130.   310 FORMAT (//' Do you wish to continue?  (Y or N):  '$)
  131.       READ (IIN,320) QUERY
  132.   320 FORMAT (A1)
  133.       IF (QUERY.NE.'Y' .AND. QUERY.NE.'y') STOP ' '
  134.       GOTO 10
  135.       END
  136. C
  137.       FUNCTION ALGAMA(S)
  138. C
  139. C        CACM ALGORITHM 291
  140. C        LOGARITHM OF GAMMA FUNCTION
  141. C        BY PIKE, M.C. AND HILL, I.D.
  142. C
  143. C        TRANSLATED INTO FORTRAN BY
  144. C                   Gerard E. Dallal
  145. C                   USDA-Human Nutrition Research Center on Aging
  146. C                              at Tufts University
  147. C                   711 Washington Street
  148. C                   Boston, MA  02111
  149. C
  150. C
  151.       X = S
  152.       ALGAMA = 0.0
  153.       IF (X .LE. 0.0) RETURN
  154.       F = 0.0
  155.       IF (X .GE. 7.0) GOTO 30
  156. C
  157.       F = 1.0
  158.       Z = X - 1.0
  159. C
  160.    10 Z = Z + 1.0
  161.       IF (Z .GE. 7.0) GOTO 20
  162.       X = Z
  163.       F = F * Z
  164.       GOTO 10
  165. C
  166.    20 X = X + 1
  167.       F = -ALOG(F)
  168. C
  169.    30 Z = 1.0 / X ** 2
  170.       ALGAMA = F + ( X - 0.5) * ALOG(X) - X + 0.918938533204673 +
  171.      1   (((-0.000595238095238 * Z + 0.000793650793651) * Z
  172.      2   - 0.002777777777778) * Z + 0.083333333333333) / X
  173.       RETURN
  174.       END
  175. C
  176.       SUBROUTINE BIGGIE(IOBS, IMAX, ICAUSE, IOUT)
  177. C
  178.       GOTO (100,200), ICAUSE
  179. C
  180.   100 WRITE(IOUT,110) IOBS, IMAX
  181.   110 FORMAT(/' ** Too many observations.'/
  182.      *        '    Your data set contains',I4/
  183.      *        '    The program is configured to handle,'I4/)
  184.       RETURN
  185. C
  186.   200 WRITE(IOUT,210) IOBS, IMAX
  187.   210 FORMAT(/
  188.      * ' ** Your data set requires a storage array of',I8,' elements.'/
  189.      * '    At present, the array contains',I7,' elements.'/)
  190.       RETURN
  191.       END
  192. C
  193.       SUBROUTINE PITMN1(DATAI ,DCOPY, NOBSX, H, MAXTOT, 
  194.      *   IOUT)
  195.       DOUBLE PRECISION PLOW, PHIGH
  196.       DIMENSION DATAI(NOBSX), H(MAXTOT), DCOPY(NOBSX)
  197. C
  198.       NOBS = 0
  199.       DO 20 I = 1, NOBSX
  200.       IF (DCOPY(I).EQ.0) GOTO 20
  201.       NOBS = NOBS + 1
  202.       DCOPY(NOBS) = DCOPY(I)
  203.    20 CONTINUE
  204.       IF (NOBS.EQ.0) STOP 'All observations are 0.'
  205. C
  206. C         ITEST = 1 -- Sign test
  207. C                 2 -- Wilcoxon signed-rank test
  208. C                 3 -- Pitman randomization test
  209. C
  210.       DO 1300 ITEST = 1, 3
  211. C
  212.       GOTO (60,50,30), ITEST
  213. C
  214.    30 DO 40 I = 1, NOBS
  215.    40 DATAI(I) = DCOPY(I)
  216.       GOTO 80
  217. C
  218.    50 CALL RANK(DATAI ,DCOPY, NOBS)
  219.       GOTO 80
  220. C
  221. C        SIGN TEST
  222. C
  223.    60 NSTAT = 0
  224.       DO 65 I = 1, NOBS
  225.       IF(DCOPY(I).GT.0) NSTAT = NSTAT + 1
  226.    65 CONTINUE
  227.       NSTAT = MIN0(NSTAT,NOBS - NSTAT)
  228. C
  229.       PVAL = 0.0
  230.       XNOBS = NOBS
  231.       XNOBS1 = NOBS + 1
  232.       DO 70 I = 0, NSTAT
  233.       TERM = ALGAMA(XNOBS1) - ALGAMA(FLOAT(I+1)) - 
  234.      *   ALGAMA(XNOBS1-FLOAT(I)) - XNOBS * ALOG(2.0)
  235.       IF (TERM.GT.-78.0) PVAL = PVAL + EXP(TERM)
  236.    70 CONTINUE
  237.       PVAL = AMIN1(1.0, 2.0 * PVAL)
  238.       GOTO 1205
  239. C
  240.    80 ISUM = 1
  241.       DO 100 I = 1, NOBS
  242.   100 ISUM = ISUM + IFIX(ABS(DATAI(I)))
  243.       IF (ISUM.LE.MAXTOT) GOTO 110
  244.       CALL BIGGIE(ISUM,MAXTOT,2,IOUT)
  245.       GOTO 1300
  246. C
  247.   110 NPDF = 1
  248.       H(1) = 1
  249. C
  250.       DO 200 K = 1, NOBS
  251. C
  252. C
  253.       INC = ABS(DATAI(K))
  254.       DO 150 I = 1, INC
  255.   150 H(NPDF+I) = 0
  256.       DO 160 I = NPDF, 1, -1
  257.   160 H(INC+I) = H(INC+I) + H(I)
  258.       NPDF = NPDF + INC
  259. C
  260.   200 CONTINUE
  261. C
  262. C        CALCULATE P-VALUE (TWO-TAILED)
  263. C
  264.       ISTAT = 0
  265.       DO 1100 I = 1, NOBS
  266.  1100 ISTAT = ISTAT + MAX1(0.0,DATAI(I))
  267.  
  268.       PLOW = 0.0
  269.       DO 1200 I = 1, ISTAT
  270.  1200 PLOW = PLOW + DBLE(H(I))
  271.       PHIGH = 2.0D0**NOBS - PLOW
  272.       PLOW = PLOW + DBLE(H(ISTAT + 1))
  273. C 
  274.       PVAL = 2.0D0 * DMIN1(PLOW,PHIGH) / 2.0D0**NOBS
  275.       PVAL = AMIN1(PVAL,1.0)
  276. C
  277.  1205 IF (ITEST.EQ.1) WRITE(IOUT,1210) PVAL
  278.  1210 FORMAT(/' Sign test:  P-value =',F7.4)
  279.       IF (ITEST.EQ.2) WRITE(IOUT,1220) PVAL
  280.  1220 FORMAT(' Wilcoxon signed-rank test:  P-value =',F7.4)
  281.       IF (ITEST.EQ.3) WRITE(IOUT,1230) PVAL
  282.  1230 FORMAT(' Pitman randomization test:  P-value =',F7.4)
  283. C
  284.  1300 CONTINUE
  285.       RETURN
  286.       END
  287.  
  288.       SUBROUTINE PITMN2(DATAI, DCOPY, M, N, NOBS, H, MAXTOT,
  289.      *   IOUT)
  290.       DOUBLE PRECISION PLOW, PHIGH
  291.       DIMENSION DATAI(NOBS), DCOPY(NOBS), H(MAXTOT)
  292. C
  293. C        ADJUST DATA
  294. C
  295.       MIN = DCOPY(1)
  296.       DO 10 I = 1, NOBS
  297.    10 IF (DCOPY(I).LT.MIN) MIN = DCOPY(I)
  298.       DO 20 I = 1, NOBS
  299.    20 DCOPY(I) = DCOPY(I) - MIN
  300. C
  301. C        ITEST = 1 -- Wilcoxon-Mann-Whitney test
  302. C                2 -- Pitman randomization test
  303. C
  304.       DO 2500 ITEST = 1, 2
  305. C
  306.       GOTO (80,110), ITEST
  307. C
  308.    80 CALL RANK(DATAI, DCOPY, NOBS)
  309. C
  310. C        SUBTRACT MINIMUM RANK
  311. C        (MIN RANK CAN BE ANYTHING DEPENDING ON TIES)
  312. C
  313.       MIN = DATAI(1)
  314.       DO 90 I = 1, NOBS
  315.    90 IF (MIN.GT.IFIX(DATAI(I))) MIN = DATAI(I)
  316.       DO 100 I = 1, NOBS
  317.   100 DATAI(I) = DATAI(I) - FLOAT(MIN)
  318.       GOTO 130
  319. C
  320.   110 DO 120 I = 1, NOBS
  321.   120 DATAI(I) = DCOPY(I)
  322.  
  323. C
  324. C        CALCULATE TEST STATISTIC BEFORE DATA ARE SORTED
  325. C
  326.   130 ISTAT = 0
  327.       DO 140 I = 1, M
  328.   140 ISTAT = ISTAT + IFIX(DATAI(I))
  329. C
  330. C        PP3 -- LINE 1
  331. C
  332.       MOBS = MIN0(M,N)
  333. C
  334. C        PP3 -- LINES 2,3
  335. C      
  336.       CALL SHELL(DATAI, NOBS)
  337. C
  338. C           IQ -- 1 + SUM OF M-LARGEST
  339. C           IU -- SUM OF M-SMALLEST
  340. C
  341.       IQ = 1
  342.       IU = 0
  343.       DO 190 K = 1, M
  344.       IQ = IQ + IFIX(DATAI(NOBS+1-K))
  345.       IU = IU + IFIX(DATAI(K))
  346.   190 CONTINUE
  347.       ISUM = IQ + IFIX(DATAI(NOBS))
  348.       ISUM = ISUM * (MOBS + 3)
  349.       IF (ISUM.LE.MAXTOT) GOTO 200
  350.       CALL BIGGIE(ISUM,MAXTOT,2,IOUT)
  351.       GOTO 2500
  352. C
  353. C           NPDFR(C) -- NUMBER OF ROWS(COLUMNS) IN H, THE PDF
  354. C
  355.   200 NPDFR = 1
  356.       NPDFC = 1
  357.       H(1) = 1
  358. C
  359. C        PP3 -- LINE 4
  360. C
  361.       DO 1000 K = 1, NOBS
  362. C
  363. C        PP3 -- LINE 5
  364. C      
  365.       INCR = DATAI(K)
  366.       INCC = 1
  367.       NPDFR2 = NPDFR + INCR
  368.       NPDFC2 = NPDFC + INCC
  369. C
  370. C        ADD ROWS
  371. C
  372.       INCR = NPDFR2 - NPDFR
  373.       IF (INCR.EQ.0) GOTO 250
  374.       DO 240 J = NPDFC, 1, -1
  375.       INDEX0 = (J-1) * NPDFR2
  376.       DO 220 I = NPDFR, 1, -1
  377.   220 H(INDEX0+I) = H((J-1)*NPDFR+I)
  378.       DO 230 I = INCR, 1, -1
  379.   230 H(INDEX0 + NPDFR + I) = 0
  380.   240 CONTINUE
  381. C
  382. C        ADD COLUMN
  383. C
  384.   250 INDEX0 = NPDFR2 * NPDFC
  385.       DO 260 I = 1, NPDFR2
  386.   260 H(INDEX0 + I) = 0
  387. C
  388.       DO 280 J = NPDFC, 1, -1
  389.       INDEX1 = (J - 1) * NPDFR2
  390.       DO 280 I = NPDFR, 1, -1
  391.       INDEX0 = (J + INCC - 1) * NPDFR2 + I + INCR
  392.       H(INDEX0) = H(INDEX0) + H(INDEX1 + I)
  393.   280 CONTINUE
  394. C
  395.       NPDFR = NPDFR2
  396.       NPDFC = NPDFC2
  397.  
  398. C
  399. C        PP3 -- LINE 6
  400. C
  401.       NPDFR2 = MIN0(NPDFR,IQ)
  402.       NPDFC2 = MIN0(NPDFC,MOBS+2)
  403. C
  404.       IF (NPDFR.EQ.NPDFR2) GOTO 300
  405.       DO 290 J = 1, NPDFC2
  406.       INDEX0 = (J - 1) * NPDFR2
  407.       INDEX1 = (J - 1) * NPDFR
  408.       DO 290 I = 1, NPDFR2
  409.   290 H(INDEX0 + I) = H(INDEX1 + I)
  410.       NPDFR = NPDFR2
  411.   300 NPDFC = NPDFC2
  412. C
  413. C        PP3 -- LINE 7
  414. C
  415.       IF (K .LE. MOBS) GOTO 1000
  416. C
  417. C        PP3 -- LINE 8
  418. C
  419.       NPDFC = NPDFC - 1
  420.  
  421. C
  422. C        PP3 -- LINE 9
  423. C
  424.       IF (K .LE. NOBS-MOBS) GOTO 1000
  425. C
  426. C        PP3 -- LINE 10
  427. C
  428.       NPDFC2 = NPDFC - 1
  429.       INCR = DATAI(K+MOBS-NOBS)
  430.       NPDFR2 = NPDFR - INCR
  431.       DO 310 J = 1, NPDFC2
  432.       INDEX0 = (J - 1) * NPDFR2
  433.       INDEX1 = J * NPDFR + INCR
  434.       DO 310 I = 1, NPDFR2
  435.   310 H(INDEX0 + I) = H(INDEX1 + I)
  436.       NPDFR = NPDFR2
  437.       NPDFC = NPDFC2
  438.  
  439. C
  440. C        PP3 -- LINE 11
  441. C
  442.  1000 CONTINUE
  443. C
  444. C        PP3 -- LINE 12
  445. C
  446.       IF (NPDFC .NE. 1) STOP 12
  447. C
  448.       DO 1100 I = 1, IU
  449.  1100 H(NPDFR + I) = 0
  450.       DO 1110 I = NPDFR, 1, -1
  451.  1110 H(I + IU) = H(I)
  452.       NPDFR = IU + NPDFR
  453. C
  454.       DO 1120 I = 1, IU
  455.  1120 H(I) = 0
  456. C
  457.       INCR = IQ - NPDFR
  458.       IF (INCR .LE. 0) GOTO 1160
  459.       DO 1150 I = 1, INCR
  460.  1150 H(NPDFR + I) = 0
  461. C
  462.  1160 NPDFR = IQ
  463. C
  464. C        PP3 -- LINE 13
  465. C
  466.       IF (M .LE. N) GOTO 2000
  467. C
  468. C        PP3 -- LINE 14
  469. C
  470. C          DROP
  471. C
  472.       NPDFR = NPDFR - IU
  473.       DO 1200 I = 1, NPDFR
  474.  1200 H(I) = H(I + IU)
  475. C
  476. C          ROTATE
  477. C
  478.       NPDFR2 = NPDFR / 2
  479.       DO 1210 I = 1, NPDFR2
  480.       DUM = H(I)
  481.       H(I) = H(NPDFR+1-I)
  482.  1210 H(NPDFR+1-I) = DUM
  483. C
  484. C          TAKE
  485. C
  486.       NPDFR2 = IQ + IFIX(DATAI(MOBS+1))
  487.       INCR = NPDFR2 - NPDFR
  488.       IF (INCR) 1220, 2000, 1250
  489. C
  490.  1220 INCR = -INCR
  491.       NPDFR = NPDFR2
  492.       DO 1230 I = 1, NPDFR
  493.  1230 H(I) = H(I+INCR)
  494.       GOTO 2000
  495. C
  496.  1250 DO 1260 I = NPDFR, 1, -1
  497.  1260 H(INCR+I) = H(I)
  498.       DO 1270 I = 1, INCR
  499.  1270 H(I) = 0
  500.       NPDFR = NPDFR2
  501. C
  502. C        PP3 -- ALL DONE.
  503. C
  504.  2000 ITOT = 0
  505.       DO 2100 I = 1, NOBS
  506.  2100 ITOT = ITOT + IFIX(DATAI(I))
  507. C
  508. C        TACK ZEROS ONTO END OF PDF
  509. C
  510.       INCR = ITOT - IQ + 1
  511.       IF (INCR.EQ.0) GOTO 2200
  512. C
  513.       DO 2120 I = 1, INCR
  514.  2120 H(NPDFR+I) = 0
  515.       NPDFR = NPDFR + INCR
  516. C
  517. C        ADJUST LENGTH
  518. C
  519.  2200 NPDFR2 = ITOT+1
  520.       INCR = NPDFR2 - NPDFR
  521.       IF (INCR) 2210, 2300, 2250
  522. C
  523.  2210 INCR = -INCR
  524.       NPDFR = NPDFR2
  525.       DO 2220 I = 1, NPDFR
  526.  2220 H(I) = H(I+INCR)
  527.       GOTO 2300
  528. C
  529.  2250 DO 2260 I = NPDFR, 1, -1
  530.  2260 H(INCR+I) = H(I)
  531.       DO 2270 I = 1, INCR
  532.  2270 H(I) = 0
  533.       NPDFR = NPDFR2
  534. C
  535.  2300 ISTAT1 = ISTAT+1
  536.       PLOW = 0.0
  537.       DO 2310 I = 1, ISTAT1
  538.  2310 PLOW = PLOW + DBLE(H(I))
  539. C
  540.       PHIGH = 0.0
  541.       DO 2320 I = ISTAT1, NPDFR
  542.  2320 PHIGH = PHIGH + DBLE(H(I))
  543. C
  544.       PVAL = 2.0D0 * DMIN1(PLOW, PHIGH)
  545.       PVAL = AMIN1(1.0, PVAL / EXP(ALGAMA(FLOAT(NOBS + 1)) - 
  546.      *   ALGAMA(FLOAT(M + 1)) - ALGAMA(FLOAT(N + 1))))
  547. C
  548.       IF (ITEST.EQ.1) WRITE(IOUT,2420) PVAL
  549.  2420 FORMAT(/' Wilcoxon-Mann-Whitney test:  P-value =',F7.4)
  550.       IF (ITEST.EQ.2) WRITE(IOUT,2430) PVAL
  551.  2430 FORMAT(' Pitman randomization test:  P-value =',F7.4)
  552.  2500 CONTINUE
  553.       RETURN   
  554.       END
  555. C
  556.       FUNCTION PROBST(T,IDF)
  557. C
  558. C         ALGORITHM AS 3  J.R.S.S. C, (1968) VOL. 17, NO. 2.
  559. C
  560. C         STUDENT T PROBABILITY (LOWER TAIL)
  561. C         USES METHOD DUE TO D.B. OWEN (BIOMETRIKA VOL. 52 1965 DEC. P438)
  562. C
  563.       DATA G1/0.3183098862/
  564. C
  565.       F=IDF
  566.       A=T/SQRT(F)
  567.       B=F/(F+T**2)
  568.       IM2=IDF-2
  569.       IOE=IDF-2*(IDF/2)
  570.       S=1.0
  571.       C=1.0
  572.       KS=2+IOE
  573.       FK=KS
  574.       IF(IM2-2) 6,7,7
  575.     7 DO 8 K=KS,IM2,2
  576.       C=C*B*(FK-1.0)/FK
  577.       S=S+C
  578.     8 FK=FK+2.
  579.     6 IF(IOE) 1,1,2
  580.     1 PROBST=.5+.5*A*SQRT(B)*S
  581.       GO TO 3
  582.     2 IF (IDF-1) 4,4,5
  583.     4 S=0.0
  584.     5 PROBST=.5+(A*B*S+ATAN(A))*G1
  585.     3 RETURN
  586.       END
  587.  
  588.       SUBROUTINE RANK(DATAI, DCOPY, N)
  589. C
  590. C        DCOPY:  DATA
  591. C        DATA:  INPUT --  SCRATCH
  592. C                 OUTPUT --  RANKS
  593.  
  594.       DIMENSION DATAI(N), DCOPY(N)
  595. C
  596.       DO 200 I = 1,N
  597.       NTIE = 0
  598.       NRANK = 0
  599.       IDUM = ABS(DCOPY(I))
  600.  
  601.       DO 100 K = 1,N
  602.       KDUM = ABS(DCOPY(K))
  603.       IF (IDUM.EQ.KDUM) NTIE = NTIE + 1
  604.       IF (IDUM.GT.KDUM) NRANK = NRANK + 1
  605.   100 CONTINUE
  606.  
  607.       DATAI(I) = 2 * NRANK + 2 + (NTIE - 1)
  608.   200 CONTINUE
  609.  
  610.       DO 210 I = 1, N
  611.   210 DATAI(I) = SIGN(DATAI(I),DCOPY(I))
  612. C
  613. C        IF ALL RANKS ARE EVEN DIVIDE BY TWO
  614. C
  615.       DO 220 I = 1, N
  616.       IF (MOD(IFIX(DATAI(I)),2).EQ.1) RETURN
  617.   220 CONTINUE
  618.       DO 230 I = 1,N
  619.   230 DATAI(I) = DATAI(I) / 2.0
  620.  
  621.       RETURN
  622.       END
  623.  
  624.       SUBROUTINE READ(IIN, IOUT, IFIN, DATA, DATAG, MAXOBS, NOBS, 
  625.      *   M, N, NGRP)
  626. C
  627.       CHARACTER FNAME*50, QUERY*1
  628.       DIMENSION DATA(MAXOBS), DATAG(MAXOBS)
  629. C
  630.   100 WRITE(IOUT,110)
  631.   110 FORMAT(/' Does this problem involve paired data (P)'/
  632.      *        '            or two independent samples (I)?  ',$)
  633.       READ (IIN,220) QUERY
  634.       NGRP = 0
  635.       IF (QUERY.EQ.'P' .OR. QUERY.EQ.'p') NGRP = 1
  636.       IF (QUERY.EQ.'I' .OR. QUERY.EQ.'i') NGRP = 2
  637.       IF (NGRP.EQ.0) GOTO 100
  638. C
  639.   200 WRITE (IOUT,210)
  640.   210 FORMAT(/' Will data be entered through the keyboard (K)'/
  641.      *        '             or read from an external file (F)?  '$)
  642.       READ (IIN,220) QUERY
  643.   220 FORMAT(A1)
  644.       IF (QUERY.EQ.'K' .OR. QUERY.EQ.'k') GOTO 400
  645.       IF (QUERY.NE.'F' .AND. QUERY.NE.'f') GOTO 200
  646. C
  647.   230 WRITE (IOUT,240)
  648.   240 FORMAT (/' Enter filename:  '$)
  649.       READ (IIN,'(A)') FNAME
  650.       OPEN (IFIN, FILE = FNAME, STATUS = 'OLD', IOSTAT = IOCHK)
  651.       IF (IOCHK.NE.0) GOTO 230
  652. C
  653.   245 WRITE (IOUT,250)
  654.   250 FORMAT (' Enter the number of variables in the file:  '$)
  655.       READ (IIN,*,ERR=245) NVAR
  656.       KVARA = 1
  657.       KVARG = 0
  658.       IF (NVAR.EQ.1) GOTO 280
  659. C
  660. C        1.2 ADDITION...  DIFFERENCES ORIGINAL VARIABLES
  661. C
  662.       IF (NGRP.EQ.2) GOTO 260
  663.   251 WRITE(IOUT,252)
  664.   252 FORMAT (' Does the files contain paired values (P)',
  665.      *   ' or differences (D)?:  '$)
  666.       READ (IIN,220) QUERY
  667.       IF (QUERY.EQ.'D' .OR. QUERY.EQ.'d') GOTO 260
  668.       IF (QUERY.NE.'P' .AND. QUERY.NE.'p') GOTO 251
  669.   253 WRITE(IOUT,254)
  670.   254 FORMAT(' Which variables are paired?:  '$)
  671.       READ(IIN,*,ERR=253) KVARA,KVARG
  672.       IF(KVARA.LT.1 .OR. KVARA.GT.NVAR .OR. KVARA.EQ.KVARG .OR.
  673.      *   KVARG.LT.1 .OR. KVARG.GT.NVAR) GOTO 253
  674. C
  675.       KVAR1 = MIN0(KVARA,KVARG)
  676.       KVAR2 = MAX0(KVARA,KVARG) 
  677.       NOBS = 0
  678.   255 READ (IFIN,*,END=256) (DUMMY,I=1,KVAR1-1),DATAX,
  679.      *   (DUMMY,I=KVAR1+1,KVAR2-1),DATAY,(DUMMY,I=KVAR2+1,NVAR)
  680.       NOBS = NOBS + 1
  681.       IF (NOBS.GT.MAXOBS) CALL BIGGIE(NOBS, MAXOBS, 1, IOUT)
  682. C
  683.       DATA(NOBS) = DATAX - DATAY
  684.       GOTO 255
  685. C
  686.   256 CLOSE (IFIN, STATUS = 'KEEP')
  687.       RETURN
  688. C
  689.   260 WRITE (IOUT,270)
  690.   270 FORMAT (' Enter the number of the variable to be analyzed:  '$)
  691.       READ (IIN,*,ERR=260) KVARA
  692.       IF (KVARA.GT.NVAR) GOTO 260
  693.       IF (NGRP.EQ.1) GOTO 280
  694.   271 WRITE (IOUT,272)
  695.   272 FORMAT (' Enter the number of the grouping variable ',
  696.      *   '(0, if none):  '$)
  697.       READ (IIN,*,ERR=271) KVARG
  698.       IF (KVARG.GT.NVAR) GOTO 271
  699.       IF (KVARG.EQ.0) GOTO 280
  700. C
  701. C        GROUPING VARIABLE PRESENT
  702. C
  703.       KVAR1 = MIN0(KVARA,KVARG)
  704.       KVAR2 = MAX0(KVARA,KVARG) 
  705.       NOBS = 0
  706.   273 READ (IFIN,*,END=274) (DUMMY,I=1,KVAR1-1),DATAX,
  707.      *   (DUMMY,I=KVAR1+1,KVAR2-1),DATAY,(DUMMY,I=KVAR2+1,NVAR)
  708.       NOBS = NOBS + 1
  709.       IF (NOBS.GT.MAXOBS) CALL BIGGIE(NOBS, MAXOBS, 1, IOUT)
  710. C
  711. C        BLOCK IF BELOW
  712. C
  713.       IF (KVARA.LT.KVARG) THEN
  714.          DATA(NOBS) = DATAX
  715.          DATAG(NOBS) = DATAY
  716.       ELSE
  717.          DATA(NOBS) = DATAY
  718.          DATAG(NOBS) = DATAX
  719.       ENDIF
  720. C
  721.       GOTO 273
  722. C
  723.   274 CLOSE (IFIN, STATUS = 'KEEP')
  724.       CALL SHELL2(DATAG,DATA,NOBS)
  725.       IF(DATAG(NOBS).EQ.DATAG(1)) STOP ' *** Only one group...'
  726. C
  727.       M = 1
  728.       DO 275 I = 2, NOBS-1
  729.       IF (DATAG(I).EQ.DATAG(1)) M = M + 1
  730.       IF (DATAG(I).NE.DATAG(1) .AND. DATAG(I).NE.DATAG(NOBS))
  731.      *   STOP ' *** More than two groups...'
  732.   275 CONTINUE
  733.       N = NOBS - M
  734.       RETURN
  735. C
  736.   280 NOBS = 0
  737.       IF (NGRP.EQ.1) GOTO 300
  738.   285 WRITE (IOUT,290)
  739.   290 FORMAT (' Enter number of observations in each group:  '$) 
  740.       READ(IIN,*,ERR=285) M, N
  741.       NOBSX = M + N
  742. C
  743. C        SORRY ABOUT STATEMENT 300.  IT USES FORTRAN-77
  744. C        CONVENTIONS THAT IF 1>KVAR-1 OR KVAR+1>NVAR THEN
  745. C        THAT LOOP WILL NOT BE EXECUTED.
  746. C
  747. C        THIS CODE IS NEEDED TO INSURE THE THE APPROPRIATE 
  748. C        NUMBER OF RECORDS ARE SKIPPED IF THERE ARE MORE THAN 
  749. C        ONE RECORD PER CASE.
  750. C
  751.   300 READ (IFIN,*,END=350) (DUMMY,I=1,KVARA-1),DATAX,
  752.      *   (DUMMY,I=KVARA+1,NVAR)
  753.       NOBS = NOBS + 1
  754.       IF (NOBS.GT.MAXOBS) CALL BIGGIE(NOBS, MAXOBS, 1, IOUT)
  755.       DATA(NOBS) = DATAX
  756.       GOTO 300
  757. C
  758.   350 CLOSE (IFIN, STATUS = 'KEEP')
  759.       IF (NGRP.EQ.1 .OR. NOBS.EQ.NOBSX) RETURN
  760.       WRITE(IOUT,360) NOBSX, NOBS
  761.   360 FORMAT(/' ** Sum of group sizes is',I4,'.  File contains',I4,
  762.      *   ' observations.')
  763.       STOP ' '
  764. C
  765.   400 IF (NGRP.EQ.2) GOTO 630
  766. C
  767. C        1.2 ADDITION... COMPUTE DIFFERENCES
  768. C
  769.   405 WRITE(IOUT,408)
  770.   408 FORMAT(/' Will you be entering pairs (P) or ',
  771.      *   'differences (D)?:  '$)
  772.       READ (IIN,220) QUERY
  773.       IF(QUERY.EQ.'D' .OR. QUERY.EQ.'d') GOTO 500
  774.       IF(QUERY.NE.'P' .AND. QUERY.NE.'p') GOTO 405
  775. C
  776.   410 WRITE(IOUT,420)
  777.   420 FORMAT(/' Enter the number of pairs:  '$)
  778.       READ (IIN,*,ERR=410) NOBS
  779.       IF (NOBS.GT.MAXOBS) CALL BIGGIE(NOBS, MAXOBS, 1, IOUT)
  780.       WRITE(IOUT,425)
  781.   425 FORMAT(/' Enter data, one pair at a time:')
  782. C
  783.       DO 460 I = 1, NOBS
  784.       GOTO 450
  785.   430 WRITE (IOUT,440)
  786.   440 FORMAT (' Error... re-enter pair:')
  787.   450 READ (IIN,*,ERR=430) DATAX, DATAY
  788.       DATA(I) = DATAX - DATAY
  789.   460 CONTINUE
  790.       RETURN
  791. C  
  792.   500 WRITE (IOUT,510)
  793.   510 FORMAT (/' Enter the number of observations:  '$)
  794.       READ (IIN,*,ERR=500) NOBS
  795.       IF (NOBS.GT.MAXOBS) CALL BIGGIE(NOBS, MAXOBS, 1, IOUT)
  796.       GOTO 650
  797. C
  798.   630 WRITE (IOUT,290)
  799.       READ(IIN,*,ERR=630) M, N
  800.       NOBS = M + N
  801.  
  802.   650 WRITE (IOUT,660)
  803.   660 FORMAT(/' Enter data:  '$)
  804.       READ (IIN,*,ERR=650) (DATA(I),I=1,NOBS)
  805.       RETURN
  806.       END
  807.  
  808.       SUBROUTINE SHELL(X, N)
  809. C        SHELL SORT WITH DECREMENTS SUGGESTED BY KNUTH
  810.       INTEGER H0(12), HP1, H, S, T
  811.       DIMENSION X(N)
  812.       DATA H0(1), H0(2), H0(3), H0(4), H0(5), H0(6)/
  813.      *         1,     4,    13,    40,   121,   364/
  814.       DATA H0(7), H0(8), H0(9), H0(10), H0(11), H0(12)/
  815.      *      1093,  3280,  9841,  29524,  88573, 265720/
  816.       NN = N
  817.       DO 1 T = 1, 10
  818.       IF (H0(T+2).GE.NN) GOTO 2
  819.     1 CONTINUE
  820.       T = 11
  821.     2 DO 6 S = 1, T
  822.       H = H0(T+1-S)
  823.       HP1 = H + 1
  824.       DO 5 J = HP1, NN
  825.       I = J - H
  826.       XM = X(J)
  827.     3 IF (XM.GE.X(I)) GOTO 4
  828.       X(I+H) = X(I)
  829.       I = I - H
  830.       IF (I.GT.0) GOTO 3
  831.     4 X(I+H) = XM
  832.     5 CONTINUE
  833.     6 CONTINUE
  834.       RETURN
  835.       END
  836. C
  837.       SUBROUTINE SHELL2(X, Y, N)
  838. C        SHELL SORT WITH DECREMENTS SUGGESTED BY KNUTH
  839. C        Y GOES ALONG FOR THE RIDE
  840.       INTEGER H0(12), HP1, H, S, T
  841.       DIMENSION X(N), Y(N)
  842.       DATA H0(1), H0(2), H0(3), H0(4), H0(5), H0(6)/
  843.      *         1,     4,    13,    40,   121,   364/
  844.       DATA H0(7), H0(8), H0(9), H0(10), H0(11), H0(12)/
  845.      *      1093,  3280,  9841,  29524,  88573, 265720/
  846.       NN = N
  847.       DO 1 T = 1, 10
  848.       IF (H0(T+2).GE.NN) GOTO 2
  849.     1 CONTINUE
  850.       T = 11
  851.     2 DO 6 S = 1, T
  852.       H = H0(T+1-S)
  853.       HP1 = H + 1
  854.       DO 5 J = HP1, NN
  855.       I = J - H
  856.       XM = X(J)
  857.       YM = Y(J)
  858.     3 IF (XM.GE.X(I)) GOTO 4
  859.       X(I+H) = X(I)
  860.       Y(I+H) = Y(I)
  861.       I = I - H
  862.       IF (I.GT.0) GOTO 3
  863.     4 X(I+H) = XM
  864.       Y(I+H) = YM
  865.     5 CONTINUE
  866.     6 CONTINUE
  867.  
  868.       RETURN
  869.       END
  870.