home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 473 / multi.arc / PC-MULTI.FOR < prev    next >
Text File  |  1985-10-16  |  18KB  |  604 lines
  1. C        
  2. C                                  PC-MULTI
  3. C                            Multiple Comparisons
  4. C                                Version 1.0
  5. C                              October 15, 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-MULTI and its author.  
  29. C
  30. C       Neither PC-Multi 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.       DIMENSION DATA(20), XNOBS(20), CIM(20,20), CIL(20,20), 
  36.      *   CIU(20,20), CILH(20), CIUH(20)
  37. C
  38. C        LINE IS A CHARACTER*>=NCELL VARIABLE
  39. C
  40.       CHARACTER*40 LINE*40, FNAME*50, GNAME(20)*6, QUERY*1
  41.       LOGICAL QFILE
  42.       DATA MAXGRP /20/, IIN /0/, IOUT /0/, IWOUT /11/, QFILE /.FALSE./
  43. C
  44.       WRITE(IOUT,1)
  45.     1 FORMAT (//36X,'PC-MULTI'/
  46.      *   30X,'Multiple Comparisons'/
  47.      *   34X,'Version 1.0'/32X,'October 15, 1985'//
  48.      *   32X,'Gerard E. Dallal'/
  49.      *   17X,'USDA Human Nutrition Research Center on Aging'/
  50.      *   30X,'at Tufts University'/29X,'711 Washington Street'/
  51.      *   31X,'Boston, MA  02111'///)
  52. c
  53.       WRITE(IOUT,10)
  54.    10 FORMAT(/' Do you wish to save the results of this session?  ',
  55.      *  '(Y or N):  '$)
  56.       READ (IIN,20) QUERY
  57.    20 FORMAT (A1)
  58.       IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QFILE = .TRUE.
  59.       IF (.NOT. QFILE) GOTO 90
  60.       WRITE (IOUT,30)
  61.    30 FORMAT (' Enter output filename:  '$)
  62.       READ (IIN,'(A)') FNAME
  63.       OPEN (IWOUT, FILE = FNAME, STATUS = 'NEW')
  64. C
  65.    90 CALL READ(IIN, IOUT, IWOUT, QFILE, DATA, ERROR, 
  66.      *   XNOBS, MAXGRP, NGRP, GNAME, V)
  67. C
  68.       R = NGRP
  69.       IF (V.GT.0.0) GOTO 120
  70. C
  71.       DO 100 I = 1, NGRP
  72.   100 V = V + XNOBS(I) - 1.0
  73. C
  74.   120 IF (QFILE) WRITE(IWOUT,130)
  75.       WRITE(IOUT,130)
  76.   130 FORMAT (' Enter the level of confidence (.GE. 0.90 .AND. ',
  77.      *   '.LE. 0.99):  '$)
  78.       READ (IIN,*) PVAL
  79.       IF (PVAL.LT.0.90 .OR. PVAL.GT.0.99) GOTO 120
  80.       IF(QFILE) WRITE(IWOUT,*) PVAL
  81.       WRITE(IOUT,140)
  82.   140 FORMAT(//' Confidence intervals based on Studentized range ',
  83.      *   'statistic...'///
  84.      *   T10,'VS',T23,'D-HSD',T35,'DIFF',T47,'D+HSD'/)
  85.       IF (QFILE) WRITE(IWOUT,140)
  86. C
  87. C        HM -- HARMONIC MEAN, USED LATER
  88. C
  89.       QCRIT = QTRNG(PVAL, V, R, IFAULT)
  90.       HM = 1.0 / XNOBS(NGRP)
  91. C
  92.       DO 200 I = 1, NGRP - 1
  93.       XNOBSI = 1.0 / XNOBS(I)
  94.       HM = HM + XNOBSI
  95. C
  96.       DO 200 K = 2, NGRP
  97.       PM = QCRIT * SQRT(ERROR / 2.0 * (XNOBSI + 1.0 / XNOBS(K)))
  98.       CIM(I,K) = DATA(I) - DATA(K)
  99.       CIL(I,K) = CIM(I,K) - PM
  100.       CIU(I,K) = CIM(I,K) + PM
  101.   200 CONTINUE
  102. C
  103.       HM = FLOAT(NGRP) / HM
  104. C
  105. C        GET MIN AND MAX
  106. C
  107.       CMIN = 0.0
  108.       CMAX = 0.0
  109.       DO 210 I = 1, NGRP-1
  110.       DO 210 K = I+1, NGRP
  111.       IF (CIL(I,K).LT.CMIN) CMIN = CIL(I,K)
  112.       IF (CIU(I,K).GT.CMAX) CMAX = CIU(I,K)
  113.   210 CONTINUE
  114. C
  115. C        SET UP GRAPH
  116. C
  117.       NCELL = 20
  118.       XNCELL = NCELL
  119.       XRNG = CMAX - CMIN
  120.       NZERO = 1.0 + (-CMIN) / XRNG * XNCELL
  121.       NZERO = MIN0(NZERO, NCELL)
  122.       DO 220 I = 1, NCELL
  123.   220 LINE(I:I) = ' '
  124. C
  125.       DO 500 I = 1, NGRP-1
  126.       DO 500 K = I+1, NGRP
  127.       NL = 1.0 + (CIL(I,K) - CMIN) / XRNG * XNCELL
  128.       NL = MIN0(NL, NCELL)
  129. C
  130.       NM = 1.0 + (CIM(I,K) - CMIN) / XRNG * XNCELL
  131.       NM = MIN0(NM, NCELL)
  132. C
  133.       NU = 1.0 + (CIU(I,K) - CMIN) / XRNG * XNCELL
  134.       NU = MIN0(NU, NCELL)
  135. C
  136. C        insert -'s
  137. C
  138.       DO 240 J = NL+1, NU-1
  139.   240 LINE(J:J) = '-'
  140. C
  141. C        insert 0 bar
  142. C
  143.       LINE(NZERO:NZERO) = '|'
  144. C
  145. C        insert ()[]
  146. C
  147.       IF (CIL(I,K).GT.0.0) LINE(NL:NL) = '('
  148.       IF (CIL(I,K).LE.0.0) LINE(NL:NL) = '['
  149.       IF (CIU(I,K).LT.0.0) LINE(NU:NU) = ')'
  150.       IF (CIU(I,K).GE.0.0) LINE(NU:NU) = ']'
  151.       LINE(NM:NM) = 'X'
  152. C
  153.       WRITE (IOUT,250) GNAME(I),GNAME(K),CIL(I,K),
  154.      *   CIM(I,K),CIU(I,K),LINE
  155.       IF (QFILE) WRITE (IWOUT,250) GNAME(I),GNAME(K),CIL(I,K),
  156.      *   CIM(I,K),CIU(I,K),LINE
  157.   250 FORMAT (1X,2A8,2X,3G12.5,4X,A40)
  158. C
  159. C        clear
  160. C
  161.       DO 280 J = NL, NU
  162.   280 LINE(J:J) = ' '
  163. C
  164.   500 CONTINUE
  165. C
  166.       WRITE (IOUT,560) ('.',J=1,NCELL)
  167.       IF (QFILE) WRITE(IWOUT,560) ('.',J=1,NCELL)
  168.   560 FORMAT (59X,40A1)
  169.       WRITE (IOUT,560) '|',(' ',J=1,NCELL-2),'|'
  170.       IF (QFILE) WRITE(IWOUT,560) '|',(' ',J=1,NCELL-2),'|'
  171.       WRITE (IOUT,562) CMIN, CMAX
  172.       IF (QFILE) WRITE(IWOUT,562) CMIN, CMAX
  173.   562 FORMAT (53X,G12.5,8X,G12.5)
  174. C
  175. C        USE HARMONIC MEAN
  176. C
  177. C        HHSD -- HALF AN HONEST SIGNIFICANT DIFFERENCE
  178. C
  179.       HHSD =  0.5 * QCRIT * SQRT(ERROR / HM)
  180.       WRITE (IOUT,610) HM
  181.   610 FORMAT (///' Using harmonic mean (',G12.5,')...'///
  182.      *   T13,'M-.5*HSD',T27,'MEAN',T37,'M+.5*HSD'/)
  183.       IF (QFILE) WRITE (IWOUT,610) HM
  184. C
  185.       DO 620 I = 1, NGRP
  186.       CILH(I) = DATA(I) - HHSD
  187.       CIUH(I) = DATA(I) + HHSD
  188.   620 CONTINUE
  189. C
  190. C        GET MIN AND MAX
  191. C
  192.       CMIN = CILH(1)
  193.       CMAX = CIUH(1)
  194.       DO 630 I = 2, NGRP
  195.       IF (CILH(I).LT.CMIN) CMIN = CILH(I)
  196.       IF (CIUH(I).GT.CMAX) CMAX = CIUH(I)
  197.   630 CONTINUE
  198. C
  199. C        SET UP GRAPH
  200. C
  201.       NCELL = 28
  202.       XNCELL = NCELL
  203.       XRNG = CMAX - CMIN
  204.       NZERO = 1.0 + (-CMIN) / XRNG * XNCELL
  205.       NZERO = MIN0(NZERO, NCELL)
  206.       DO 640 I = 1, NCELL
  207.   640 LINE(I:I) = ' '
  208. C
  209.       DO 700 I = 1, NGRP
  210.       NL = 1.0 + (CILH(I) - CMIN) / XRNG * XNCELL
  211.       NL = MIN0(NL, NCELL)
  212. C
  213.       NM = 1.0 + (DATA(I) - CMIN) / XRNG * XNCELL
  214.       NM = MIN0(NM, NCELL)
  215. C
  216.       NU = 1.0 + (CIUH(I) - CMIN) / XRNG * XNCELL
  217.       NU = MIN0(NU, NCELL)
  218. C
  219. C        insert -'s
  220. C
  221.       DO 650 J = NL+1, NU-1
  222.   650 LINE(J:J) = '-'
  223. C
  224. C        insert 0 bar
  225. C
  226.       LINE(NZERO:NZERO) = '|'
  227. C
  228. C        insert ()[]
  229. C
  230.       LINE(NL:NL) = '('
  231.       LINE(NU:NU) = ')'
  232.       LINE(NM:NM) = 'X'
  233. C
  234.       WRITE (IOUT,655) GNAME(I),CILH(I),DATA(I),CIUH(I),LINE
  235.   655 FORMAT (1X,A8,2X,3G12.5,4X,A40)
  236.       IF (QFILE) WRITE (IWOUT,655) 
  237.      *   GNAME(I),CILH(I),DATA(I),CIUH(I),LINE
  238. C
  239. C
  240. C        clear
  241. C
  242.       DO 660 J = NL, NU
  243.   660 LINE(J:J) = ' '
  244. C
  245.   700 CONTINUE
  246. C
  247.       WRITE (IOUT,760) ('.',J=1,NCELL)
  248.   760 FORMAT (51X,40A1)
  249.       IF (QFILE) WRITE(IWOUT,760) ('.',J=1,NCELL)
  250.       WRITE (IOUT,760) '|',(' ',J=1,NCELL-2),'|'
  251.       IF (QFILE) WRITE(IWOUT,760) '|',(' ',J=1,NCELL-2),'|'
  252.       WRITE (IOUT,762) CMIN, CMAX
  253.   762 FORMAT (45X,G12.5,16X,G12.5)
  254.       IF (QFILE) WRITE(IWOUT,762) CMIN, CMAX
  255. C
  256.       WRITE (IOUT,3020)
  257.  3020 FORMAT (///' Do you wish to continue?  (Y or N):  '$)
  258.       READ (IIN,20) QUERY
  259.       IF (QUERY.NE.'Y' .AND. QUERY.NE.'y') STOP ' '
  260.       IF (QFILE) WRITE (IWOUT, 3030)
  261.  3030 FORMAT(1H1)
  262.       GOTO 90
  263.  
  264.       END
  265. C
  266.       SUBROUTINE READ(IIN, IOUT, IWOUT, QFILE, DATA, ERROR, XNOBS, 
  267.      *   MAXGRP, NGRP, GNAME, DFE)
  268. C
  269.       LOGICAL QEXT, QFILE
  270.       CHARACTER FNAME*50, QUERY*1, GNAME(MAXGRP)*6, D3*6
  271.       DIMENSION DATA(MAXGRP), XNOBS(MAXGRP)
  272.  
  273.       DATA IFIN /10/, QEXT/.FALSE./
  274. C
  275.   200 WRITE (IOUT,210)
  276.   210 FORMAT(/' Will data be entered through the keyboard (K)'/
  277.      *        '             or read from an external file (F)?  '$)
  278.       READ (IIN,220) QUERY
  279.   220 FORMAT(A1)
  280.       IF (QUERY.EQ.'K' .OR. QUERY.EQ.'k') GOTO 400
  281.       IF (QUERY.NE.'F' .AND. QUERY.NE.'f') GOTO 200
  282. C
  283.       IF (QFILE) WRITE(IWOUT,240)
  284.   230 WRITE (IOUT,240)
  285.   240 FORMAT (/' Enter filename:  '$)
  286.       READ (IIN,'(A)') FNAME
  287.       IF (QFILE) WRITE(IWOUT,250) FNAME
  288.   250 FORMAT(1X,A50)
  289.       OPEN (IFIN, FILE = FNAME, STATUS = 'OLD', ERR = 230)
  290. C
  291.       READ (IFIN,260) QUERY, ERROR
  292.   260 FORMAT(BN,A1,F18.0)
  293.       IF((QUERY.EQ.'E' .OR. QUERY.EQ.'e')
  294.      *   .AND. ERROR.GT.0.0) GOTO 280
  295. C
  296.       WRITE (IOUT,270)
  297.   270 FORMAT (/' ** ** ** Error term misspecified...')
  298.       STOP ' '
  299. C
  300.   280 READ (IFIN,260) QUERY, DFE
  301.       IF((QUERY.EQ.'D' .OR. QUERY.EQ.'d')
  302.      *   .AND. DFE.GE.0.0) GOTO 300
  303. C
  304.       WRITE (IOUT,290)
  305.   290 FORMAT (/' ** ** ** Degrees of freedom misspecified...')
  306.       STOP ' '
  307. C
  308.   300 NGRP = 0
  309.       DO 340 KOUNT = 1, MAXGRP+1
  310.       READ (IFIN,*,END=350) D1, D2, D3
  311.       IF (D2.LE.0.0) GOTO 340
  312.       NGRP = NGRP + 1
  313.       DATA(NGRP) = D1
  314.       XNOBS(NGRP) = D2
  315.       GNAME(NGRP) = D3
  316.   340 CONTINUE
  317.       WRITE( IOUT,341)
  318.   341 FORMAT (/' *** Too many means in the data file')
  319.       STOP ' '
  320. C
  321.   350 CLOSE (IFIN, STATUS = 'KEEP')
  322.       IF (NGRP.GT.1) RETURN
  323.       WRITE(IOUT,360)
  324.   360 FORMAT (/' *** Too few means in the data file')
  325.       STOP ' '
  326. C
  327.   400 WRITE (IIN,410)
  328.       IF (QFILE) WRITE(IWOUT,410)
  329.   410 FORMAT (/' Enter error variance :  '$)
  330.       READ (IIN,*) ERROR
  331.       IF (ERROR.LE.0.0) GOTO 400
  332.       IF (QFILE) WRITE(IWOUT,*) ERROR      
  333. C
  334.       IF (QFILE) WRITE(IWOUT,430)
  335.   420 WRITE (IIN,430)
  336.   430 FORMAT (' Enter degrees of freedom for error variance :  '$)
  337.       READ (IIN,*) DFE
  338.       IF (DFE.LT.0.0) GOTO 420
  339.       IF (QFILE) WRITE(IWOUT,*) ERROR      
  340. C
  341.       WRITE (IOUT,440)
  342.       IF (QFILE) WRITE(IWOUT,440)
  343.   440 FORMAT (' Enter the number of groups:  '$)
  344.       READ (IIN,*) NGRP
  345.       IF (QFILE) WRITE(IWOUT,*) NGRP
  346.       IF (NGRP.LE.MAXGRP) GOTO 450
  347. C
  348.       WRITE(IOUT,445)
  349.   445 FORMAT(/' *** Too many means in the data file')
  350.       STOP ' '
  351. C
  352.   450 WRITE (IOUT,460)
  353.       IF (QFILE) WRITE(IWOUT,460)
  354.   460 FORMAT(' Enter triples of mean, sample size, and group name:')
  355.       READ (IIN,*) (DATA(I),XNOBS(I),GNAME(I),I=1,NGRP)
  356.       IF (.NOT. QFILE) RETURN
  357. C
  358.       DO 500 I = 1, NGRP
  359.   500 WRITE(IWOUT,*) DATA(I), XNOBS(I)
  360. C
  361.       RETURN
  362.       END
  363. C
  364.       FUNCTION ALNORM(X,UPPER)
  365. C
  366. C        ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3
  367. C
  368. C        EVALUATES THE TAIL AREA OF THE STANDARD NORMAL CURVE
  369. C        FROM X TO INFINITY IF UPPER IS .TRUE. OR
  370. C        FROM MINUS INFINITY TO X IF UPPER IS .FALSE.
  371. C
  372.       REAL LTONE, UTZERO, ZERO, HALF, ONE, CON, Z, Y, X
  373.       LOGICAL UPPER, UP
  374. C
  375. C        LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR
  376. C        COMPUTER (SEE INTRODUCTORY TEXT)
  377. C
  378.       DATA LTONE, UTZERO/7.0, 18.66/
  379.       DATA ZERO, HALF, ONE, CON/0.0, 0.5, 1.0, 1.28/
  380.       UP = UPPER
  381.       Z = X
  382.       IF (Z .GE. ZERO) GOTO 10
  383.       UP = .NOT. UP
  384.       Z = -Z
  385.    10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20
  386.       ALNORM = ZERO
  387.       GOTO 40
  388.    20 Y = HALF * Z * Z
  389.       IF (Z .GT. CON) GOTO 30
  390. C
  391.       ALNORM = HALF - Z * (0.398942280444 - 0.399903438504 * Y /
  392.      1   (Y + 5.75885480458 - 29.8213557808 /
  393.      2   (Y + 2.62433121679 + 48.6959930692 /
  394.      3   (Y + 5.92885724438))))
  395.       GOTO 40
  396. C
  397.    30 ALNORM = .398942280385 * EXP(-Y) /
  398.      1   (Z - 3.8052E-8 + 1.00000615302 /
  399.      2   (Z + 3.98064794E-4 + 1.98615381364 /
  400.      3   (Z - 0.151679116635 + 5.29330324926 /
  401.      4   (Z + 4.8385912808 - 15.1508972451 /
  402.      5   (Z + 0.742380924027 + 30.789933034 / (Z + 3.99019417011))))))
  403. C
  404.    40 IF (.NOT. UP) ALNORM = ONE - ALNORM
  405.       RETURN
  406.       END
  407.       FUNCTION GAUINV(P,IFAULT)
  408. C
  409. C              ALGORITHM AS 70 APPL. STATIST. (1974) VOL. 23, NO.1
  410. C              GAUINV FINDS PERCENTAGE POINTS OF THE NORMAL DISTRIBUTION
  411. C
  412.       DATA ZERO,ONE,HALF,ALIMIT/0.0, 1.0, 0.5, 1.0E-20/
  413. C
  414.       DATA             P0,  P1,             P2,                P3
  415.      1   / -.322232431088, -1., -.342242088547, -.204231210245E-1/
  416. C
  417.       DATA                P4,               Q0,            Q1
  418.      1   / -.453642210148E-4, .993484626060E-1, .588581570495/
  419. C
  420.       DATA            Q2,            Q3,              Q4
  421.      1   / .531103462366, .103537752850, .38560700634E-2/
  422. C
  423.       IFAULT=1
  424.       GAUINV=ZERO
  425.       PS=P
  426.       IF (PS.GT.HALF) PS=ONE-PS
  427.       IF (PS.LT.ALIMIT) RETURN
  428.       IFAULT=0
  429.       IF (PS.EQ.HALF) RETURN
  430.       YI=SQRT(ALOG(ONE/(PS*PS)))
  431.       GAUINV=YI+((((YI*P4+P3)*YI+P2)*YI+P1)*YI+P0)
  432.      1         /((((YI*Q4+Q3)*YI+Q2)*YI+Q1)*YI+Q0)
  433.       IF (P.LT.HALF) GAUINV=-GAUINV
  434.       RETURN
  435.       END
  436. C
  437.       FUNCTION PRTRNG(Q, V, R, IFAULT)
  438. C
  439. C        ALGORITHM AS 190  APPL. STATIST. (1983) VOL.32, NO.2
  440. C        INCLUDES MODIFICATION  FROM VOL.34, NO.1
  441. C
  442. C        EVALUATES THE PROBABILITY FORM 0 TO Q FOR A STUDENTIZED
  443. C        RANGE HAVING V DEGREES OF FREEDOM AND R SAMPLES.
  444. C
  445.       REAL Q, V, R, VW(30), QW(30), PCUTJ, PCUTK, STEP, VMAX, ZERO,
  446.      *   FIFTH, HALF, ONE, TWO, CV1, CV2, CVMAX, CV(4)
  447.       DATA PCUTJ, PCUTK, STEP, VMAX /0.00003E0, 0.0001E0, 0.45E0,
  448.      *  120.0E0/, ZERO, FIFTH, HALF, ONE, TWO /0.0E0, 0.2E0, 0.5E0,
  449.      *  1.0E0, 2.0E0/, CV1, CV2, CVMAX /0.193064705E0, 0.293525326E0,
  450.      *  0.39894228E0/, CV(1), CV(2), CV(3), CV(4) /0.318309886E0,
  451.      *  -0.268132716E-2, 0.347222222E-2, 0.833333333E-1/
  452.       DATA JMIN, JMAX, KMIN, KMAX /3, 15, 7, 15/
  453. C
  454. C        CHECK INITIAL VALUES
  455. C
  456.       PRTRNG = ZERO
  457.       IFAULT = 0
  458.       IF (V .LT. ONE .OR. R .LT. TWO) IFAULT = 1
  459.       IF (Q .LE. ZERO .OR. IFAULT .EQ. 1) GOTO 99
  460. C
  461.       G = STEP * R ** (-FIFTH)
  462.       GMID = HALF * ALOG(R)
  463.       R1 = R - ONE
  464.       C = ALOG(R * G * CVMAX)
  465.       IF ( V .GT. VMAX) GOTO 20
  466.  
  467. C
  468.       H = STEP * V ** (-HALF)
  469.       V2 = V * HALF
  470.       IF (V .EQ. ONE) C = CV1
  471.       IF (V .EQ. TWO) C = CV2
  472.       IF (.NOT.(V .EQ. ONE .OR. V .EQ. TWO)) C = SQRT(V2) * CV(1)
  473.      *   / (ONE + ((CV(2) / V2 + CV(3)) / V2 + CV(4)) / V2)
  474.       C = ALOG(C * R * G * H)
  475. C
  476.    20 GSTEP = G
  477.       QW(1) = -ONE
  478.       QW(JMAX + 1) = -ONE
  479.       PK1 = ONE
  480.       PK2 = ONE
  481.       DO 28 K = 1, KMAX
  482.       GSTEP = GSTEP - G
  483.    21 GSTEP = -GSTEP
  484.       GK = GMID + GSTEP
  485.       PK = ZERO
  486.       IF (PK2 .LE. PCUTK .AND. K .GT. KMIN) GOTO 26
  487.       W0 = C - GK * GK * HALF
  488.       PZ = ALNORM(GK, .TRUE.)
  489.       X = ALNORM(GK - Q, .TRUE.) - PZ
  490.       IF (X .GT. ZERO) PK = EXP(W0 + R1 * ALOG(X))
  491.       IF (V .GT. VMAX) GOTO 26
  492. C
  493.       JUMP = -JMAX
  494.    22 JUMP = JUMP + JMAX
  495.       DO 24 J = 1, JMAX
  496.       JJ = J + JUMP
  497.       IF (QW(JJ) .GT. ZERO) GOTO 23
  498.       HJ = H * FLOAT(J)
  499.       IF (J .LT. JMAX) QW(JJ + 1) = -ONE
  500.       EHJ = EXP(HJ)
  501.       QW(JJ) = Q * EHJ
  502.       VW(JJ) = V * (HJ + HALF - EHJ * EHJ * HALF)
  503. C
  504.    23 PJ = ZERO
  505.       X = ALNORM(GK - QW(JJ), .TRUE.) - PZ
  506.       IF (X .GT. ZERO) PJ = EXP(W0 + VW(JJ) + R1 * ALOG(X))
  507.       PK = PK + PJ
  508.       IF (PJ .GT. PCUTJ) GOTO 24
  509.       IF (JJ .GT. JMIN .OR. K .GT. KMIN) GOTO 25
  510.    24 CONTINUE
  511.    25 H = -H
  512.       IF (H .LT. ZERO) GOTO 22
  513. C
  514.    26 PRTRNG = PRTRNG + PK
  515.       IF (K .GT. KMIN .AND. PK .LE. PCUTK .AND. PK1 .LE. PCUTK) GOTO 99
  516.       PK2 = PK1
  517.       PK1 = PK
  518.       IF (GSTEP .GT. ZERO) GOTO 21
  519.    28 CONTINUE
  520. C
  521.    99 RETURN
  522.       END
  523. C
  524.       FUNCTION QTRNG(P, V, R, IFAULT)
  525. C
  526. C        WRITTEN TO REPLACE...
  527. C        ALGORITHM AS 190.1  APPL. STATIST.  (1983) VOL.32, NO.2
  528. C        INCLUDES MODIFICATION  FROM VOL.34, NO.1
  529. C
  530. C        APPROXIMATES THE QUANTILE P FOR A STUDENTIZED RANGE
  531. C        DISTRIBUTION HAVING V DEGREES OF FREEDOM AND R SAMPLES
  532. C        FOR PROBABILITY P, P.GE.0.90 .AND. P.LE.0.99
  533. C
  534.       REAL P, V, R, PCUT, P75, P80, P90, P99, P995, P175, ONE, TWO, FIVE
  535.       DATA JMAX, PCUT, P75, P80, P90, P99, P995, P175, ONE, TWO, FIVE /
  536.      *  8, 0.0001E0, 0.75E0, 0.80E0, 0.90E0, 0.99E0, 0.995E0, 1.75E0,
  537.      *  1.0E0, 2.0E0, 5.0E0/
  538. C
  539. C        CHECK INPUT PARAMETERS
  540. C
  541.       IFAULT = 0
  542.       NFAULT = 0
  543.       IF (V .LT. ONE .OR. R .LT. TWO) IFAULT = 1
  544.       IF (P .LT. P90 .OR. P .GT. P99) IFAULT = 2
  545.       IF (IFAULT .NE. 0) GOTO 99
  546. C
  547. C        OBTAIN INITIAL VALUES
  548. C
  549.       Q1 = QTRNG0(P, V, R, NFAULT)
  550.       IF (NFAULT .NE. 0) GOTO 99
  551.       P1 = PRTRNG(Q1, V, R, NFAULT)
  552.       IF (NFAULT .NE. 0) GOTO 99
  553.       IF(ABS(P1-P).GT.PCUT) GOTO 10
  554.       QTRNG = Q1
  555.       RETURN
  556. C
  557.    10 SIGN = 1.0
  558.       IF (P1.GT.P) SIGN = -1.0
  559. C
  560.    90 Q2 = Q1 + SIGN
  561.       P2 = PRTRNG(Q2, V, R, NFAULT)
  562.       IF (NFAULT .NE. 0) GOTO 99
  563.       IF(ABS(P2-P).GT.PCUT) GOTO 110
  564.       QTRNG = Q2
  565.       RETURN
  566. C
  567.   110 IF ((P1-P) * (P2-P) .LT. 0.0 ) GOTO 120
  568.       Q1 = Q2
  569.       P1 = P2
  570.       GOTO 90
  571. C
  572. C        BINARY SEARCH
  573. C
  574.   120 QMIN = AMIN1(Q1, Q2)
  575.       QMAX = AMAX1(Q1, Q2)
  576. C
  577.   140 QTRNG = (QMIN + QMAX) / 2.0
  578.       PNEW = PRTRNG(QTRNG, V, R, NFAULT)
  579.       IF (NFAULT .NE. 0) GOTO 99
  580.       IF (ABS(PNEW-P) .LE. PCUT) RETURN
  581.       IF (PNEW.GT.P) QMAX = QTRNG
  582.       IF (PNEW.LT.P) QMIN = QTRNG
  583.       GOTO 140
  584. C
  585.    99 IF (NFAULT .NE .0) IFAULT = 9
  586.       RETURN
  587.       END
  588. C
  589.       FUNCTION QTRNG0(P, V, R, IFAULT)
  590. C
  591. C        ALGORITHM 190.2  APPL. STATIST.  (1983) VOL.32, NO.2
  592. C
  593.       REAL P, V, R, Q, T, VMAX, HALF, ONE, FOUR, C1, C2, C3, C4, C5
  594.       DATA VMAX, HALF, ONE, FOUR, C1, C2, C3, C4, C5 /120.0E0, 0.5E0,
  595.      *  1.0E0, 4.0E0, 0.8843E0, 0.2368E0, 1.214E0, 1.208E0, 1.4142E0/
  596. C
  597.       T = GAUINV(HALF + HALF * P, IFAULT)
  598.       IF (V .LT. VMAX) T = T + (T * T * T + T) / V / FOUR
  599.       Q = C1 - C2 * T
  600.       IF (V .LT. VMAX) Q = Q - C3 / V + C4 * T / V
  601.       QTRNG0 = T * (Q * ALOG(R - ONE) + C5)
  602.       RETURN
  603.       END
  604.