home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 473 / multi.arc / PC-SIZE.FOR < prev    next >
Text File  |  1985-11-21  |  44KB  |  1,431 lines
  1. C        
  2. C                                  PC-SIZE
  3. C                  A Program for Sample Size Determinations
  4. C                                Version 2.0
  5. C                              August 8, 1985
  6. C
  7. C                              Gerard E. Dallal
  8. C        
  9. C               USDA Human Nutrition Research Center on Aging
  10. C                            at Tufts University
  11. C                           711 Washington Street
  12. C                             Boston, MA  02111
  13. C
  14. C                                    and
  15. C        
  16. C                    Tufts University School of Nutrition
  17. C                             132 Curtis Street
  18. C                             Medford, MA  02155
  19. C        
  20. C
  21. C
  22. C                                  NOTICE
  23. C
  24. C       Documentation and original code copyright 1985 by  Gerard  E.  
  25. C       Dallal.  Reproduction of material for non-commercial purposes 
  26. C       is   permitted,   without  charge,   provided  that  suitable 
  27. C       reference is made to PC-SIZE and its author.  
  28. C
  29. C       Neither PC-SIZE nor its documentation should be  modified  in 
  30. C       any way without permission from the author,  except for those 
  31. C       changes  that  are  essential  to  move  PC-SIZE  to  another 
  32. C       computer.  
  33. C
  34. C       Please acknowledge PC-SIZE in any manuscript  that  uses  its 
  35. C       calculations.  
  36. C
  37. C
  38. C        IIN   -- INPUT UNIT NUMBER (SCREEN)
  39. C        IOUT  -- OUTPUT UNIT NUMBER (SCREEN)
  40. C        IWOUT -- SAVED OUTPUT UNIT NUMBER 
  41. C        NMAX0 -- LARGEST PERMISSIBLE INTEGER
  42. C        NEMAX -- LARGEST NUMBER OF EFFECTS THAT CAN BE ENTERED
  43. C                    INDIVIDUALLY  (DIMENSION OF 'EFFECT')
  44. C
  45. C        OPEN STATEMENT LOCATED JUST BEFORE STATEMENT 10
  46. C
  47.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  48.       LOGICAL QFILE, QAB
  49.       CHARACTER*1 QUERY, BLANK
  50.       CHARACTER FNAME*40, CF2C0*18, CEFECT*18
  51.       DIMENSION EFFECT(50)
  52. C
  53.       DATA IIN /0/, IOUT /0/, IWOUT /11/, NMAX0 /2147483647/
  54.       DATA QFILE /.FALSE./, BLANK /' '/, NEMAX /50/
  55. C
  56. C        INITIALIZE CHOICES
  57. C
  58.       MODE = 1
  59.       ALPHA = 0.05
  60.       RPOWER = 0.80
  61.       NG = 2
  62.       NROW = 2
  63.       NCOL = 2
  64.       EDIFF = 0.0
  65.       WSD = 0.0
  66.       F1 = 0.0
  67.       F2C1 = 0.0
  68.       F2C0 = 0.0
  69.       XLAM0 = 0.0
  70.       XLAM01 = 0.0
  71.       IMETH = 3
  72.       RANGE = 0.0
  73.       RHO = 0.0
  74.       P1 = 0.0
  75.       P2 = 0.0
  76. C
  77.       WRITE(IOUT,1)
  78.     1 FORMAT (//36X,'PC-SIZE'/
  79.      *   20X,'A Program for Sample Size Determinations'/
  80.      *   34X,'Version 2.0'/32X,'August 8, 1985'//32X,'Gerard E. Dallal'/
  81.      *   17X,'USDA Human Nutrition Research Center on Aging'/
  82.      *   30X,'at Tufts University'/29X,'711 Washington Street'/
  83.      *   31X,'Boston, MA  02111'///37X,'NOTICE'//9X,'Please '
  84.      *   'acknowledge  PC-SIZE  in any manuscript  that uses its',/
  85.      *   9X,'calculations.')
  86.  
  87.       WRITE(IOUT,2)
  88.     2 FORMAT(//' Do you wish to save the results of this session?  ',
  89.      *  '[N]:  '$)
  90.       READ (IIN,3) QUERY
  91.     3 FORMAT (A1)
  92.       IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QFILE = .TRUE.
  93.       IF (.NOT. QFILE) GOTO 10
  94.       WRITE (IOUT,4)
  95.     4 FORMAT (' Enter output filename:  '$)
  96.       READ (IIN,'(A)') FNAME
  97.       OPEN (IWOUT, FILE = FNAME, STATUS = 'NEW')
  98.       WRITE (IWOUT,1)
  99.  
  100.    10 WRITE(IOUT,20) 
  101.       IF (QFILE) WRITE (IWOUT,20) 
  102.    20 FORMAT(///' This program operates in seven modes:'//
  103.      *   '     1  --  single factor design'/
  104.      *   '     2  --  two factor design'/
  105.      *   '     3  --  randomized blocks design'/
  106.      *   '     4  --  paired t-test'/
  107.      *   '     5  --  generic F'/
  108.      *   '     6  --  correlation coefficient'/
  109.      *   '     7  --  proportions'/)
  110.       IF (QFILE) WRITE (IWOUT,22) MODE
  111.    21 WRITE(IOUT,22) MODE
  112.    22 FORMAT(' Choose a mode [',I1,']:  '$)
  113.       READ (IIN,50,ERR=21) XMODE
  114.       IF (XMODE.NE.DINT(XMODE) .OR. XMODE.LT.0.0D0
  115.      *   .OR. XMODE.GT.7.0D0) GOTO 21
  116.       IF (XMODE.NE.0.0D0) MODE = XMODE
  117.       IF (QFILE) WRITE(IWOUT,29) MODE
  118.    29 FORMAT (1X,I6)
  119.       IF (QFILE) WRITE (IWOUT,40) ALPHA
  120.    30 WRITE (IOUT,40) ALPHA
  121.    40 FORMAT (/' Enter the level of the test [',F5.3,']:  '$)
  122.       READ (IIN,50,ERR=30) DALPHA
  123.       IF (DALPHA.LT.0.0D0 .OR. DALPHA.GE.1.0D0) GOTO 30
  124.       IF (DALPHA.NE.0.0D0) ALPHA = DALPHA
  125.       IF (QFILE) WRITE (IWOUT,49) ALPHA
  126.    49 FORMAT (1X,G12.5)
  127.    50 FORMAT (BN,F18.0)
  128.  
  129.       IF (RPOWER.LT.1.0D0 .AND. QFILE) WRITE (IWOUT,61) RPOWER
  130.       IF (RPOWER.GE.1.0D0 .AND. QFILE) WRITE (IWOUT,62) RPOWER
  131.    60 IF (RPOWER.LT.1.0D0) WRITE (IOUT,61) RPOWER
  132.    61 FORMAT (' Enter the required power [',F5.3,']:  '$)
  133.       IF (RPOWER.GE.1.0D0) WRITE (IOUT,62) RPOWER
  134.    62 FORMAT (' Enter the required power [',F8.0,']:  '$)
  135.       READ (IIN,50,ERR=60) DPOWER
  136.       IF (DPOWER.LT.0.0D0 .OR. 
  137.      *   (DPOWER.GT.1.D0 .AND. DPOWER.NE.DINT(DPOWER))) GOTO 60
  138.       IF (MODE.EQ.6 .AND. DPOWER.GE.1.0D0 .AND. DPOWER.LT.3.0D0)
  139.      *   GOTO 60
  140.       IF (DPOWER.NE.0.0D0) RPOWER = DPOWER
  141.       IF (QFILE) WRITE (IWOUT,49) RPOWER
  142.       IF (RPOWER.LT.1.0D0) GOTO 79
  143.  
  144.       NSTART = RPOWER
  145.       IF (QFILE) WRITE (IWOUT,69) NSTART
  146.       WRITE (IOUT,69) NSTART
  147.    69 FORMAT (' Power calculations will start with a sample size of:',
  148.      *   I6)
  149.       IF (QFILE) WRITE (IWOUT,71) 
  150.    70 WRITE (IOUT,71) 
  151.    71 FORMAT(' Enter the largest sample size to be considered:  '$)
  152.       READ (IIN,50,ERR=70) XNMAX
  153.       IF (XNMAX.NE.DINT(XNMAX) .OR. XNMAX.LT.RPOWER) GOTO 70
  154.       NMAX = XNMAX
  155.       IF (QFILE) WRITE (IWOUT,29) NMAX
  156.       NSTEP = 1
  157.       IF (NMAX.EQ.NSTART .OR. NMAX.EQ.NSTART+1) GOTO 79
  158. C
  159.       IF (QFILE) WRITE (IWOUT,73) 
  160.    72 WRITE (IOUT,73) 
  161.    73 FORMAT (' Enter increment:  '$)
  162.       READ (IIN,50,ERR=72) XNSTEP
  163.       IF(XNSTEP.LT.1.0D0 .OR. XNSTEP.NE.DINT(XNSTEP)) GOTO 72
  164.       NSTEP = XNSTEP
  165.       IF (QFILE) WRITE (IWOUT,29) NSTEP
  166. C
  167. C        GET NUMERATOR DEGREES OF FREEDOM
  168. C
  169.    79 GOTO (80, 110, 210, 230, 300, 2000, 3000), MODE
  170. C
  171. C        SINGLE FACTOR DESIGN
  172. C
  173.    80 IF (QFILE) WRITE (IWOUT,90) NG
  174.    81 WRITE (IOUT,90) NG
  175.    90 FORMAT (' Enter the number of groups [',I2,']:  '$)
  176.       READ (IIN,50,ERR=81) XNG
  177.       IF (XNG.LT.0.0D0 .OR. XNG.NE.DINT(XNG)) GOTO 81
  178.       IF (XNG.NE.0.0D0) NG = XNG
  179.       XNG = NG
  180.       IF (NG.LT.2) GOTO 81
  181.       IF (QFILE) WRITE (IWOUT,29) NG
  182.       F1 = NG - 1
  183.       GOTO 400
  184. C
  185. C        TWO FACTOR DESIGN
  186. C
  187.   110 IF (QFILE) WRITE (IWOUT,120) NROW
  188.   111 WRITE (IOUT,120) NROW
  189.   120 FORMAT (/' Enter the number of levels for Factor A [',I2,']:  '$)
  190.       READ (IIN,50,ERR=111) XNROW
  191.       IF (XNROW.LT.0.0D0 .OR. XNROW.NE.DINT(XNROW)) GOTO 111
  192.       IF (XNROW.NE.0.0D0) NROW = XNROW
  193.       IF (NROW.LT.2) GOTO 111
  194.       IF (QFILE) WRITE (IWOUT,29) NROW
  195.       XNROW = NROW
  196.  
  197.       IF (QFILE) WRITE (IWOUT,140) NCOL
  198.   130 WRITE (IOUT,140) NCOL
  199.   140 FORMAT (' Enter the number of levels for Factor B [',I2,']:  '$)
  200.       READ (IIN,50,ERR=130) XNCOL
  201.       IF (XNCOL.LT.0.0D0 .OR. XNCOL.NE.DINT(XNCOL)) GOTO 130
  202.       IF (XNCOL.NE.0.0D0) NCOL = XNCOL
  203.       IF (NCOL.LT.2) GOTO 130
  204.       IF (QFILE) WRITE (IWOUT,29) NCOL
  205.       XNCOL = NCOL
  206.       NG = NROW
  207.       XNG = NG
  208.       F1 = NG - 1
  209.  
  210.       IF (QFILE) WRITE (IWOUT,160)
  211.   150 WRITE (IOUT,160)
  212.   160 FORMAT (' Do interactions appear in the model and the ANOVA ',
  213.      *   'table?  (Y or N):  '$)
  214.       READ (IIN,3) QUERY
  215.       IF (QUERY.NE.'Y' .AND. QUERY.NE.'y' .AND. QUERY.NE.'N'
  216.      *   .AND. QUERY.NE.'n') GOTO 150
  217.       QAB = .FALSE.
  218.       IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QAB = .TRUE.
  219.       IF (QFILE) WRITE (IWOUT,161) QUERY
  220.   161 FORMAT(1X,A1)
  221.       GOTO 400
  222. C
  223. C        RANDOMIZED BLOCKS DESIGN
  224. C      
  225.   210 IF (QFILE) WRITE (IWOUT,220) NG
  226.   211 WRITE (IOUT,220) NG
  227.   220 FORMAT(' Enter the number of levels for the treatment factor [',
  228.      *   I2,']:  '$)
  229.       READ (IIN,50,ERR=211) XNG
  230.       IF (XNG.LT.0.0D0 .OR. XNG.NE.DINT(XNG)) GOTO 211
  231.       IF (XNG.NE.0.0) NG = XNG
  232.       XNG = NG
  233.       IF (NG.LT.2) GOTO 211
  234.       IF (QFILE) WRITE (IWOUT,29) NG
  235.       F1 = NG - 1
  236.       GOTO 400
  237. C
  238. C        PAIRED T-TEST
  239. C
  240.   230 IF (EDIFF.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,240) EDIFF
  241.   240 FORMAT (' Enter expected difference [',G12.5,']:  '$)
  242.       IF (EDIFF.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,241) 
  243.   241 FORMAT (' Enter expected difference [0.00]:  '$)
  244.   242 IF (EDIFF.NE.0.0D0) WRITE (IOUT,240) EDIFF
  245.       IF (EDIFF.EQ.0.0D0) WRITE (IOUT,241) 
  246.       READ (IIN,50,ERR=242) XDIFF
  247.       IF (XDIFF.NE.0.0D0) EDIFF = XDIFF
  248.       IF (EDIFF.EQ.0.0D0) GOTO 242
  249.       IF (QFILE) WRITE (IWOUT,49) EDIFF
  250.       F1 = 1
  251.       NG = 2
  252.       XNG = 2.0
  253.  
  254.       IF (WSD.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,250) WSD
  255.   250 FORMAT (' Enter estimate of standard deviation of difference [',
  256.      *   G12.5,']:  '$)
  257.       IF (WSD.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,255) 
  258.   255 FORMAT (' Enter estimate of standard deviation of difference [',
  259.      *   '0.00]:  '$)
  260.   256 IF (WSD.NE.0.0D0) WRITE (IOUT,250) WSD
  261.       IF (WSD.EQ.0.0D0) WRITE (IOUT,255) 
  262.       READ (IIN,50,ERR=256)  XWSD
  263.       IF (XWSD.LT.0.0D0) GOTO 256
  264.       IF (XWSD.GT.0.0D0) WSD = XWSD
  265.       IF (QFILE) WRITE (IWOUT,49) WSD
  266.       IF (WSD.EQ.0.0D0) GOTO 256
  267.       XLAM0 =  (EDIFF / WSD) ** 2
  268.       GOTO 900
  269. C
  270. C        GENERIC
  271. C
  272.  
  273.   300 NF1 = F1
  274.       IF (QFILE) WRITE (IWOUT,301) NF1
  275.   301 FORMAT (' Enter the numerator degrees of freedom [',I2,']:  '$)
  276.   310 NF1 = F1
  277.   311 WRITE (IOUT,301) NF1
  278.       READ (IIN,50,ERR=311) XF1
  279.       IF (XF1.LT.0.0D0 .OR. XF1.NE.DINT(XF1)) GOTO 311
  280.       IF (XF1.GT.0.0D0) F1 = XF1
  281.       IF (F1.LE.0.0D0) GOTO 310
  282.       IF (QFILE) WRITE (IWOUT,49) F1
  283.  
  284.       WRITE (IOUT,313) 
  285.       IF (QFILE) WRITE (IWOUT,313) 
  286.   313 FORMAT (' The denominator degrees of freedom are calculated as'/
  287.      *   '   the linear function  C1 * N - C0  of the sample size, N.')
  288.       NF2C1 = F2C1
  289.       IF (QFILE) WRITE (IWOUT,317) NF2C1
  290.   314 NF2C1 = F2C1
  291.   316 WRITE (IOUT,317) NF2C1
  292.   317 FORMAT(' Enter C1 [',I3,']:  '$)
  293.       READ (IIN,50,ERR=316) XF2C1
  294.       IF (XF2C1.LT.0.0D0 .OR. XF2C1.NE.DINT(XF2C1)) GOTO 316
  295.       IF (XF2C1.NE.0.0D0) F2C1 = XF2C1
  296.       IF (F2C1.LE.0.0D0) GOTO 314
  297.       IF (QFILE) WRITE (IWOUT,49) F2C1
  298.  
  299.       NF2C0 = F2C0
  300.       IF (QFILE) WRITE (IWOUT,322) NF2C0
  301.   320 NF2C0 = F2C0
  302.   321 WRITE (IOUT,322) NF2C0
  303.   322 FORMAT (' Enter C0 [',I3,']:  '$)
  304.       READ (IIN,'(A)') CF2C0
  305.  
  306. C
  307. C        CHECK FOR EMPTY
  308. C
  309.       DO 326 I = 1,18
  310.       IF (CF2C0(I:I).NE.BLANK) GOTO 327
  311.   326 CONTINUE
  312.       GOTO 330
  313.  
  314.   327 READ (CF2C0,50,ERR=321) XF2C0
  315.       IF (XF2C0.NE.DINT(XF2C0)) GOTO 320
  316.       F2C0 = XF2C0
  317.   330 IF (QFILE) WRITE (IWOUT,49) F2C0
  318.  
  319.       WRITE (IOUT,340) 
  320.       IF (QFILE) WRITE (IWOUT,340) 
  321.   340 FORMAT (' The noncentrality parameter '
  322.      *   'is calculated as  C * N.')
  323.       IF (XLAM0.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,350) 
  324.   350 FORMAT (' Enter  C [0.00]:  '$) 
  325.       IF (XLAM0.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,360) XLAM0
  326.   360 FORMAT (' Enter  C [',G12.5,']:  '$) 
  327.   365 IF (XLAM0.EQ.0.0D0) WRITE (IOUT,350) 
  328.       IF (XLAM0.NE.0.0D0) WRITE (IOUT,360) XLAM0
  329.       READ (IIN,50,ERR=365) XXLAM0
  330.       IF (XXLAM0.LT.0.0D0) GOTO 365
  331.       IF (XXLAM0.NE.0.0D0) XLAM0 = XXLAM0
  332.       IF (XLAM0.LE.0.0D0) GOTO 365
  333.       IF (QFILE) WRITE (IWOUT,49) XLAM0
  334.  
  335.       GOTO 900
  336.  
  337.   400 WRITE (IOUT,410) 
  338.       IF (QFILE) WRITE (IWOUT,410) 
  339.   410 FORMAT (/' There are three ways to specify the alternative:'//
  340.      *   '    1 -- Specify individual effects.'/
  341.      *   '    2 -- Specify a range throughout which group'/
  342.      *   '            means are uniformly distrbuted.'/
  343.      *   '    3 -- Specify average squared effect divided'/
  344.      *   '            by the error variance.'/)
  345. C    
  346.       IF (QFILE) WRITE (IWOUT,420) IMETH
  347.   415 WRITE (IOUT,420) IMETH
  348.   420 FORMAT (' Choose a method [',I1,']:  '$)  
  349. C
  350.       READ (IIN,50,ERR=415) XIMETH
  351.       IF (XIMETH.LT.0.0D0 .OR. 
  352.      *   XIMETH.GT.3.0D0 .OR. XIMETH.NE.DINT(XIMETH)) GOTO 415
  353.       IF (XIMETH.NE.0.0D0) IMETH = XIMETH
  354.       IF (QFILE) WRITE (IWOUT,29) IMETH
  355.  
  356.       GOTO (500,600,800),IMETH
  357.  
  358.   500 EAVG = 0.
  359.       IF (NG.LE.NEMAX) GOTO 502
  360.       WRITE(IOUT,501) NEMAX
  361.   501 FORMAT (/' NO MORE THAN',I3,' EFFECTS MAY BE SPECIFIED')
  362.       GOTO 400
  363.  
  364.   502 DO 520  I = 1, NG
  365.       IF (MODE.EQ.1 .AND. QFILE) WRITE (IWOUT,510) I,EFFECT(I)
  366.   510 FORMAT (' Enter estimate of effect for group ',I2,' [',
  367.      *   G12.5,']:  '$)
  368.       IF (MODE.EQ.2 .AND. QFILE) WRITE (IWOUT,511) I,EFFECT(I)
  369.   511 FORMAT (' Enter estimate of effect for Factor A, level',I2,' [',
  370.      *   G12.5,']:  '$)
  371.       IF (MODE.EQ.3 .AND. QFILE) WRITE (IWOUT,512) I,EFFECT(I)
  372.   512 FORMAT (' Enter estimate of effect for treatment ',I2,' [',
  373.      *   G12.5,']:  '$)
  374.   513 IF (MODE.EQ.1) WRITE (IOUT,510) I,EFFECT(I)
  375.       IF (MODE.EQ.2) WRITE (IOUT,511) I,EFFECT(I)
  376.       IF (MODE.EQ.3) WRITE (IOUT,512) I,EFFECT(I)
  377.  
  378.       READ (IIN,'(A)')  CEFECT
  379. C
  380. C        CHECK FOR EMPTY
  381. C
  382.       DO 514 II = 1, 18
  383.       IF (CEFECT(II:II).NE.BLANK) GOTO 515
  384.   514 CONTINUE
  385.       GOTO 519
  386.  
  387.   515 READ (CEFECT,50,ERR=513) EFFECT(I)
  388.  
  389.   519 IF (QFILE) WRITE (IWOUT,49) EFFECT(I)
  390.       EAVG = EAVG + EFFECT(I)
  391.   520 CONTINUE
  392.       EAVG = EAVG / XNG
  393.       XLAM0 = 0.
  394.       DO 530 I = 1, NG
  395.   530 XLAM0 = XLAM0 + (EFFECT(I) - EAVG) ** 2
  396.       IF (XLAM0.GT.0.0D0) GOTO 700
  397.       WRITE (IOUT,540)
  398.       IF (QFILE) WRITE (IWOUT,540)
  399.   540 FORMAT (/' EFFECTS ARE ALL EQUAL.  PLEASE RESPECIFY.'/)
  400.       GOTO 500
  401.  
  402.   600 IF (RANGE.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,610) 
  403.   610 FORMAT (' Enter range [0.00]:  '$)
  404.       IF (RANGE.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,620) RANGE
  405.   620 FORMAT (' Enter range [',G12.5,']:  '$)
  406.   630 IF (RANGE.EQ.0.0D0) WRITE (IOUT,610) 
  407.       IF (RANGE.NE.0.0D0) WRITE (IOUT,620) RANGE
  408.       READ (IIN,50,ERR=630) XRANGE
  409.       IF (XRANGE.LT.0.0) GOTO 630
  410.       IF (XRANGE.GT.0.0D0) RANGE = XRANGE
  411.       IF (RANGE.LE.0.0D0) GOTO 630
  412.       IF (QFILE) WRITE (IWOUT,49) RANGE
  413.       DELTA = RANGE / (XNG - 1.0D0)
  414.       XLAM0 = DELTA ** 2 * (XNG ** 3 - XNG) / 12.0D0 
  415.  
  416.   700 IF (WSD.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,720) WSD
  417.   720 FORMAT (' Enter estimate of within cell standard deviation [',
  418.      *   G12.5,']:  '$)
  419.       IF (WSD.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,730) 
  420.   730 FORMAT (' Enter estimate of within cell standard deviation [',
  421.      *   '0.00]:  '$)
  422.   740 IF (WSD.NE.0.0D0) WRITE (IOUT,720) WSD
  423.       IF (WSD.EQ.0.0D0) WRITE (IOUT,730) 
  424.       READ (IIN,50,ERR=740) XWSD
  425.       IF (XWSD.LT.0.0D0) GOTO 740
  426.       IF (XWSD.GT.0.0D0) WSD = XWSD
  427.       IF (WSD.LE.0.0D0) GOTO 740
  428.       IF (QFILE) WRITE (IWOUT,49) WSD
  429.       XLAM0 =  XLAM0 / WSD ** 2
  430.       GOTO 900
  431.      
  432.   800 IF (XLAM01.GT.0.0D0 .AND. QFILE) WRITE (IWOUT,810) XLAM01
  433.   810 FORMAT (' Specify average squared effect divided'/
  434.      *   '    by the error variance [',G12.5,']:  '$)
  435.       IF (XLAM01.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,820) 
  436.   820 FORMAT (' Specify average squared effect divided'/
  437.      *   '    by the error variance [0.00]:  '$)
  438.   830 IF (XLAM01.GT.0.0D0) WRITE (IOUT,810) XLAM01
  439.       IF (XLAM01.EQ.0.0D0) WRITE (IOUT,820) 
  440.       READ (IIN,50,ERR=830) XXLAM0
  441.       IF (XXLAM0.LT.0.0D0) GOTO 830
  442.       IF (XXLAM0.GT.0.0D0) XLAM01 = XXLAM0
  443.       IF (XLAM01.EQ.0.0D0) GOTO 830
  444.       IF (QFILE) WRITE (IWOUT,49) XLAM01
  445.       XLAM0 = XLAM01 * XNG
  446. C
  447.   900 IF (MODE.EQ.2) XLAM0 =  XNCOL * XLAM0
  448.       WRITE (IOUT,910) ALPHA,XLAM0
  449.       IF (QFILE) WRITE (IWOUT,910) ALPHA,XLAM0
  450.   910 FORMAT(//' Level of the test:  ',F6.4/
  451.      *   ' Non-centrality parameter: ',G12.5,' *  sample size')
  452. C
  453. C        LARGE SAMPLE APPROXIMATION
  454. C
  455.       IF (RPOWER .GE. 1.0D0) GOTO 979
  456.       ALPHAC = 1.0D0 - ALPHA
  457.       CHISQ = PPCHI2(ALPHAC, F1)
  458.       CALL CHISQL(CHISQ, F1, 1.0D0 - RPOWER, FLLRGE)
  459.       XNLRGE = FLLRGE / XLAM0 
  460.       NLRGE = XNLRGE
  461.       NLRGE1 = 0.05D0 * XNLRGE
  462. C
  463.       NSTART = MAX0(NLRGE - 10, 1)
  464.       IF (NLRGE1.GT.10) NSTART = NLRGE - NLRGE1
  465.       NMAX = NMAX0
  466.       NSTEP = DLOG10(XNLRGE) - DLOG10(50.D0)
  467.       IF (NSTEP.LT.0) NSTEP = 0
  468.       NSTEP = 10 ** NSTEP
  469.       NSTART = NSTART / NSTEP * NSTEP
  470. C
  471.   979 WRITE (IOUT,980)
  472.       IF (QFILE) WRITE (IWOUT,980)
  473.   980 FORMAT (//'     SAMPLE'/'      SIZE     POWER'/)
  474.  
  475.   989 DO 1000 N = NSTART, NMAX, NSTEP
  476.       XN = N
  477.       XLAM = XN * XLAM0
  478.       IF (MODE.EQ.1) F2 = XNG * (XN - 1.0D0)
  479.       IF (MODE.EQ.2 .AND. QAB) F2 = XNROW * XNCOL * (XN - 1.0D0)
  480.       IF (MODE.EQ.2 .AND. .NOT.QAB)
  481.      *    F2 = XNROW * XNCOL * XN - XNROW - XNCOL + 1.0D0
  482.       IF (MODE.EQ.3 .OR. MODE.EQ.4) F2 = (XNG - 1.0D0) * (XN - 1.0D0)
  483.       IF (MODE.EQ.5) F2 = F2C1 * XN - F2C0
  484.       IF (F2 .LT. 1.0D0) GOTO 1000
  485.       AA = F2 / 2.0D0
  486.       BB = F1 / 2.0D0
  487.       XALPHA = XINBTA(AA, BB, ALPHA)
  488.       FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
  489.       POWER = FNCDFC(FALPHA, F1, F2, XLAM)
  490.       WRITE (IOUT,990) N, POWER
  491.       IF (QFILE) WRITE (IWOUT,990) N, POWER
  492.   990 FORMAT (1X, I8, 3X, F8.5)
  493.       IF (POWER .GT. RPOWER) GOTO 1010
  494.  1000 CONTINUE
  495.  1010 WRITE (IOUT,1020)
  496.  1020 FORMAT (///' Do you wish to continue?  [Y]:  '$)
  497.       READ (IIN,3) QUERY
  498.       IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
  499.       GOTO 10
  500. C
  501. C        CORRELATION
  502. C
  503.  2000 IF (QFILE) WRITE (IWOUT,2080) RHO
  504.  2080 FORMAT (/' Enter the correlation coefficient'/
  505.      *   '    at the alternative [',F5.3,']:  '$)
  506.  2081 WRITE (IOUT,2080) RHO
  507.       READ (IIN,50,ERR=2081) DRHO
  508.       IF (DABS(DRHO).GE.1.0D0) GOTO 2081
  509.       IF (DRHO.NE.0.0D0) RHO = DRHO
  510.       IF (RHO.EQ.0.0D0) GOTO 2081
  511.       IF (QFILE) WRITE (IWOUT,49) RHO
  512.       IF (RPOWER.LT.1.0D0) GOTO 2200
  513. C
  514.       WRITE (IOUT,980)
  515.       IF (QFILE) WRITE (IWOUT,980)
  516.  
  517.       DO 2100 N = NSTART, NMAX, NSTEP
  518.       XN = N
  519.       R = RCRIT(ALPHA,XN - 2.0D0)
  520.       POWER = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
  521.      *   RDF(N, R, RHO, IFAULT) 
  522.       WRITE (IOUT,990) N, POWER
  523.       IF (QFILE) WRITE (IWOUT,990) N, POWER
  524.  2100 CONTINUE
  525.       WRITE (IOUT,1020)
  526.       READ (IIN,3) QUERY
  527.       IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
  528.       GOTO 10
  529.  
  530.  
  531. C
  532. C        N = 3
  533. C
  534.  2200 WRITE (IOUT,980)
  535.       IF (QFILE) WRITE (IWOUT,980)
  536.       N = 3
  537.       XN = N
  538.       R = RCRIT(ALPHA,XN - 2.0D0)
  539.       P3 = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
  540.      *   RDF(N, R, RHO, IFAULT)
  541.       WRITE (IOUT,990) N, P3
  542.       IF (QFILE) WRITE (IWOUT,990) N, P3
  543.       IF (P3.GT.RPOWER) GOTO 230
  544. C
  545. C        N = 4
  546. C
  547.       N = 4
  548.       XN = N
  549.       R = RCRIT(ALPHA,XN - 2.0D0)
  550.       P4 = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
  551.      *   RDF(N, R, RHO, IFAULT) 
  552. C
  553.       WRITE (IOUT,990) N, P4
  554.       IF (QFILE) WRITE (IWOUT,990) N, P4
  555.       IF (P4.GT.RPOWER) GOTO 2230
  556. C
  557.  2210 N = 2 * N
  558.       XN = N
  559.       R = RCRIT(ALPHA,XN - 2.0D0)
  560.       POWER = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
  561.      *   RDF(N, R, RHO, IFAULT) 
  562.       WRITE (IOUT,990) N, POWER
  563.       IF (QFILE) WRITE (IWOUT,990) N, POWER
  564.       IF (POWER.LT.RPOWER) GOTO 2210
  565. C
  566.       NHI = N
  567.       HIPOWR = POWER
  568.       NLO = N / 2
  569.  2220 NM = (NLO + NHI) / 2
  570.       XN = NM
  571.       R = RCRIT(ALPHA,XN - 2.0D0)
  572.       POWER = RDF(NM, -R, RHO, IFAULT) + 1.0D0 -
  573.      *   RDF(NM, R, RHO, IFAULT) 
  574.       IF (POWER.LT.RPOWER) NLO = NM
  575.       IF (POWER.GE.RPOWER) NHI = NM
  576.       IF (POWER.GE.RPOWER) HIPOWR = POWER
  577.       IF (NHI - NLO.GT.1) GOTO 2220
  578.       WRITE (IOUT,990) NHI, HIPOWR
  579.       IF (QFILE) WRITE (IWOUT,990) NHI, HIPOWR
  580.  
  581.  2230 WRITE (IOUT,1020)
  582.       READ (IIN,3) QUERY
  583.       IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
  584.       GOTO 10
  585. C
  586. C        PROPORTIONS -- TWO BY TWO TABLES
  587. C
  588.  3000 IF (QFILE) WRITE (IWOUT,3080) P1
  589.  3080 FORMAT (/' Enter the first proportion under'/
  590.      *   '    the alternative to equality [',F5.3,']:  '$)
  591.  3081 WRITE (IOUT,3080) P1
  592.       READ (IIN,50,ERR=3081) DP1
  593.       IF (DP1.LT.0.0D0 .OR. DP1.GE.1.0D0) GOTO 3081
  594.       IF (DP1.NE.0.0D0) P1 = DP1
  595.       IF (P1.EQ.0.0D0) GOTO 3081
  596.       IF (QFILE) WRITE (IWOUT,49) P1
  597. C
  598.       IF (QFILE) WRITE (IWOUT,3090) P2
  599.  3090 FORMAT (/' Enter the second proportion under'/
  600.      *   '    the alternative to equality [',F5.3,']:  '$)
  601.  3091 WRITE (IOUT,3090) P2
  602.       READ (IIN,50,ERR=3091) DP2
  603.       IF (DP2.LT.0.0D0 .OR. DP2.GE.1.0D0) GOTO 3091
  604.       IF (DP2.NE.0.0D0) P2 = DP2
  605.       IF (P2.EQ.0.0D0) GOTO 3091
  606.       IF (QFILE) WRITE (IWOUT,49) P2
  607. C
  608.       IF (P1.NE.P2) GOTO 3095
  609.       WRITE(IOUT,3094)
  610.  3094 FORMAT(/' The two proportions are equal.  Please respecify.'/)
  611.       GOTO 3000
  612. C
  613.  3095 PBAR = (P1 + P2) / 2.0D0
  614.       DELTA = DABS(P2 - P1)
  615. C
  616.       IF (RPOWER.LT.1.0D0) GOTO 3200
  617. C
  618.       WRITE (IOUT,980)
  619.       IF (QFILE) WRITE (IWOUT,980)
  620.  
  621.       DO 3100 N = NSTART, NMAX, NSTEP
  622.       XN = N
  623.       XM = XN - 2.0D0 / DELTA + 1.0D0 / (XN * DELTA**2)
  624.       PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
  625.       PC2 = (PC1 * DSQRT(2.0D0 * PBAR * (1.0D0 - PBAR)) - DELTA *
  626.      *   DSQRT(XM)) / DSQRT(P1 * (1.0D0 - P1) + P2 * (1.0D0 - P2))
  627.       POWER = ALNORM(PC2,.TRUE.)
  628.       WRITE (IOUT,990) N, POWER
  629.       IF (QFILE) WRITE (IWOUT,990) N, POWER
  630.  3100 CONTINUE
  631.       WRITE (IOUT,1020)
  632.       READ (IIN,3) QUERY
  633.       IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
  634.       GOTO 10
  635. C
  636.  3200 PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
  637.       PC2 = GAUINV(1.0D0 - RPOWER,IFAULT) 
  638.       XN = ((PC1 * DSQRT(2.0D0 * 
  639.      *   PBAR * (1.0D0 - PBAR)) - PC2 * DSQRT(P1 *
  640.      *   (1.0D0 - P1) + P2 * (1.0D0 - P2))) / DELTA)**2
  641.       XN = XN / 4.0D0 * (1.0D0 + DSQRT(1.0D0 + 4.0D0 /
  642.      *   (XN * DELTA)))**2 + 0.9999999999D0
  643. C
  644.       N = XN
  645.       XN = N
  646.       XM = XN - 2.0D0 / DELTA + 1.0D0 / (XN * DELTA**2)
  647.       PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
  648.       PC2 = (PC1 * DSQRT(2.0D0 * PBAR * (1.0D0 - PBAR)) - DELTA *
  649.      *   DSQRT(XM)) / DSQRT(P1 * (1.0D0 - P1) + P2 * (1.0D0 - P2))
  650.       POWER = ALNORM(PC2,.TRUE.)
  651.  
  652.       WRITE (IOUT,980)
  653.       IF (QFILE) WRITE (IWOUT,980)
  654.       WRITE (IOUT,990) N, POWER
  655.       IF (QFILE) WRITE (IWOUT,990) N, POWER
  656.  
  657.       WRITE (IOUT,1020)
  658.       READ (IIN,3) QUERY
  659.       IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
  660.       GOTO 10
  661.  
  662.       END
  663. C
  664.       DOUBLE PRECISION FUNCTION ALGAMA(S)
  665.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  666. C
  667. C        CACM ALGORITHM 291
  668. C        LOGARITHM OF GAMMA FUNCTION
  669. C        BY PIKE, M.C. AND HILL, I.D.
  670. C
  671. C        TRANSLATED INTO FORTRAN BY
  672. C                   Gerard E. Dallal
  673. C                   USDA-Human Nutrition Research Center on Aging
  674. C                              at Tufts University
  675. C                   711 Washington Street
  676. C                   Boston, MA  02111
  677. C
  678. C
  679.       X = S
  680.       ALGAMA = 0.0
  681.       IF (X .LE. 0.0D0) RETURN
  682.       F = 0.0
  683.       IF (X .GE. 7.0D0) GOTO 30
  684. C
  685.       F = 1.0
  686.       Z = X - 1.0D0
  687. C
  688.    10 Z = Z + 1.0D0
  689.       IF (Z .GE. 7.0D0) GOTO 20
  690.       X = Z
  691.       F = F * Z
  692.       GOTO 10
  693. C
  694.    20 X = X + 1.0D0
  695.       F = -DLOG(F)
  696. C
  697.    30 Z = 1.0D0 / X ** 2
  698.       ALGAMA = F + ( X - 0.5D0) * DLOG(X) - X + 0.918938533204673D0 +
  699.      1   (((-0.000595238095238D0 * Z + 0.000793650793651D0) * Z
  700.      2   - 0.002777777777778D0) * Z + 0.083333333333333D0) / X
  701.       RETURN
  702.       END
  703.  
  704.       DOUBLE PRECISION FUNCTION ALNORM(X,UPPER)
  705. C
  706. C        ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3
  707. C
  708. C        EVALUATES THE TAIL AREA OF THE STANDARD NORMAL CURVE
  709. C        FROM X TO INFINITY IF UPPER IS .TRUE. OR
  710. C        FROM MINUS INFINITY TO X IF UPPER IS .FALSE.
  711. C
  712.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  713.       DOUBLE PRECISION LTONE
  714.       LOGICAL UPPER, UP
  715. C
  716. C        LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR
  717. C        COMPUTER (SEE INTRODUCTORY TEXT)
  718. C
  719.       DATA LTONE, UTZERO/7.0D0, 18.66D0/
  720.       DATA ZERO, HALF, ONE, CON/0.0D0, 0.5D0, 1.0D0, 1.28D0/
  721.       UP = UPPER
  722.       Z = X
  723.       IF (Z .GE. ZERO) GOTO 10
  724.       UP = .NOT. UP
  725.       Z = -Z
  726.    10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20
  727.       ALNORM = ZERO
  728.       GOTO 40
  729.    20 Y = HALF * Z * Z
  730.       IF (Z .GT. CON) GOTO 30
  731. C
  732.       ALNORM = HALF - Z * (0.398942280444D0 - 0.399903438504D0 * Y /
  733.      1   (Y + 5.75885480458D0 - 29.8213557808D0 /
  734.      2   (Y + 2.62433121679D0 + 48.6959930692D0 /
  735.      3   (Y + 5.92885724438D0))))
  736.       GOTO 40
  737. C
  738.    30 ALNORM = .398942280385D0 * EXP(-Y) /
  739.      1   (Z - 3.8052D-8 + 1.00000615302D0 /
  740.      2   (Z + 3.98064794D-4 + 1.98615381364D0 /
  741.      3   (Z - 0.151679116635D0 + 5.29330324926D0 /
  742.      4   (Z + 4.8385912808D0 - 15.1508972451D0 /
  743.      5   (Z + 0.742380924027D0 + 30.789933034D0 / 
  744.      6   (Z + 3.99019417011D0))))))
  745. C
  746.    40 IF (.NOT. UP) ALNORM = ONE - ALNORM
  747.       RETURN
  748.       END
  749. C
  750.       DOUBLE PRECISION FUNCTION BETAIN(X,P,Q)
  751.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  752. C
  753. C        ALGORITHM AS 63  APPL.STATIST. (1973), VOL.22, NO.3
  754. C        MODIFIED AS PER REMARK ASR 19  (1977), VOL.26, NO.1
  755. C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
  756. C
  757. C        COMPUTES INCOMPLETE BETA FUNCTION RATIO FOR ARGUMENTS
  758. C        X BETWEEN ZERO AND ONE, P AND Q POSITIVE.
  759. C        LOG OF COMPLETE BETA FUNCTION, BETA, ASSUMED TO BE KNOWN.
  760. C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
  761. C
  762.       LOGICAL INDEX
  763. C
  764. C        DEFINE ACCURACY AND INITIALIZE
  765. C
  766.       DATA ACU /0.1D-7/
  767.       BETAIN = X 
  768. C
  769. C        TEST FOR ADMISIBILITY OF AGRUMENTS
  770. C
  771.       IF (P .LE. 0.0D0 .OR. Q .LE. 0.0D0) STOP 40
  772.       IF (X .LT. 0.0D0 .OR .X .GT. 1.0D0) STOP 41
  773.       IF (X .EQ .0.0D0 .OR .X .EQ. 1.0D0) RETURN
  774. C
  775. C     CHANGE TAIL IF NECESSARY AND DETERMINE S
  776. C
  777.       PSQ = P + Q
  778.       CX = 1.0D0 - X
  779.       IF (P .GE .PSQ * X) GOTO 1
  780.       XX = CX
  781.       CX = X
  782.       PP = Q
  783.       QQ = P
  784.       INDEX = .TRUE.
  785.       GOTO 2
  786.     1 XX = X
  787.       PP = P
  788.       QQ = Q
  789.       INDEX = .FALSE.
  790.     2 TERM = 1.0
  791.       AI = 1.0
  792.       BETAIN = 1.0
  793.       NS = QQ + CX * PSQ
  794. C
  795. C        USE SOPER'S REDUCTION FORMULAE
  796. C
  797.       RX = XX / CX
  798.     3 TEMP = QQ - AI
  799.       IF (NS .EQ. 0) RX = XX
  800.     4 TERM = TERM * TEMP * RX / (PP + AI)
  801.       BETAIN = BETAIN + TERM
  802.       TEMP = ABS(TERM)
  803.       IF (TEMP .LE. ACU .AND. TEMP .LE. ACU * BETAIN) GOTO 5
  804.       AI = AI + 1.0D0
  805.       NS = NS - 1
  806.       IF (NS .GE. 0) GOTO 3
  807.       TEMP = PSQ
  808.       PSQ = PSQ + 1.0D0
  809.       GOTO 4
  810. C
  811. C        CALCULATE RESULT
  812. C
  813.     5 BETA = ALGAMA(P) + ALGAMA(Q) - ALGAMA(P+Q)
  814.       BETAIN = BETAIN * 
  815.      *     DEXP(PP * DLOG(XX) + (QQ - 1.0D0) * DLOG(CX) - BETA) / PP
  816.       IF (INDEX) BETAIN = 1.0D0 - BETAIN
  817.       RETURN
  818.       END
  819. C
  820.       SUBROUTINE CHISQL(X, DF, FX, FL)
  821.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  822. C
  823. C
  824. C        DEFINE ACCURACY AND INITIALIZE
  825. C
  826. C        N SHOULD BE SECIFIED SUCH THAT ACC IS GREATER THAN 
  827. C        OR EQUAL TO (AU - AL) / 2 ** N
  828. C
  829. C
  830.       DATA ACC, N /1.0D-6, 30/
  831.       AL = 0.0
  832.       AINC = 80.0
  833.       AU = 80.0
  834. C
  835.       IF (DF .LE. 0.0D0) STOP 71
  836.       IF (X .LT. 0.0D0) STOP 72
  837.       IF (FX .LE. 0.0D0) STOP 73
  838. C
  839.     1 APROX = CHI2NC(X, DF, AU)
  840.       IF (FX .GT. APROX) GOTO 2
  841.       IF (FX .LT. APROX) AL = AU
  842.       AU = AU + AINC
  843.       GOTO 1
  844.     2 DO 3 J = 1, N
  845.       FL = 0.5D0 * (AL + AU)
  846.       APROX = CHI2NC(X, DF, FL)
  847.       IF ((AU - AL)/AU .LT. ACC) GOTO 4
  848.       IF (FX .LT. APROX) AL = FL
  849.       IF (FX .GE. APROX) AU = FL
  850.     3 CONTINUE
  851.     4 RETURN
  852.       END
  853. C
  854.       DOUBLE PRECISION FUNCTION CHI2NC(FQUAN, DFN, XLAM)
  855.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  856. C
  857. C        THE CDF OF THE NON-CENTRAL CHI-SQUARE DISTRIBUTION
  858. C        WITH DFN DEGREES OF FREEDOM
  859. C        AND NON-CENTRALITY PARAMETER XLAM
  860. C        EVALUATED AT FQUAN.
  861. C
  862. C        CALCULATED BY SUMMING AN 'INFINITE SERIES' 
  863. C
  864. C        THE SERIES IS A WEIGHTED SUM OF CENTRAL CHI-SQUARE
  865. C        PROBABILITIES 
  866. C        WITH WEIGHTS GIVEN BY THE POISSON DISTRIBUTION WITH 
  867. C        MEAN XLAM / 2.  SUMMATION STARTS AT XLAM / 2 AND CONTINUES
  868. C        UP AN DOWN UNTIL THE CONTRIBUTION IN BOTH DIRECTIONS IS LESS 
  869. C        THAN ACU.
  870. C
  871. C
  872.       DATA ACU / 1.0D-10/
  873. C
  874.       XLAMH = XLAM / 2.0D0
  875.       NSPLIT = XLAMH + 0.499999D0
  876.       IF (NSPLIT .LT. 1) NSPLIT = 1
  877.    
  878.       CHI2NC = 0.0
  879.       ESQ = FQUAN / 2.0D0
  880.  
  881.       XI = NSPLIT
  882.       DO 100 I = 1, NSPLIT
  883.       XI = XI - 1.0D0
  884.       TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
  885.       IF (TERM .LT. -80.0D0) GOTO 200
  886.       TERM = DEXP(TERM)
  887.       TERM = TERM * GAMAIN (ESQ, DFN/2.0D0 + XI)
  888.       CHI2NC = CHI2NC + TERM
  889.       IF (TERM .LT. ACU) GOTO 200
  890.   100 CONTINUE
  891.  
  892.   200 XI = NSPLIT - 1
  893.  
  894.   300 XI = XI + 1.0D0
  895.       TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
  896.       IF (TERM .LT. -80.0D0) RETURN
  897.       TERM = DEXP(TERM)
  898.       TERM = TERM * GAMAIN (ESQ, DFN / 2.0D0 + XI)
  899.       CHI2NC = CHI2NC + TERM
  900.       IF (TERM .LT. ACU) RETURN
  901.       GOTO 300
  902.       END
  903. C
  904.       DOUBLE PRECISION FUNCTION FNCDFC(FQUAN, DFN, DFD, XLAM)
  905.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  906. C
  907. C        THE COMPLEMENT OF CDF FOR THE NON-CENTRAL F DISTRIBUTION
  908. C        WITH DFN NUMERATOR DEGREES OF FREEDOM, DFD DENOMINATOR 
  909. C        DEGREES OF FREEDOM, AND NON-CENTRALITY PARAMETER XLAM
  910. C        EVALUATED AT FQUAN.
  911. C
  912. C        CALCULATED BY SUMMING AN 'INFINITE SERIES' 
  913. C        (GRAYBILL (1961), AN INTRODUCTION TO LINEAR STATISTICAL
  914. C        MODELS, VOLUME 1, P.79, EQU.4.5.).  
  915. C
  916. C        THE SERIES IS A WEIGHTED SUM OF CENTRAL F PROBABILITIES 
  917. C        WITH WEIGHTS GIVEN BY THE POISSON DISTRIBUTION WITH 
  918. C        MEAN XLAM / 2.  SUMMATION STARTS AT XLAM / 2 AND CONTINUES
  919. C        UP AN DOWN UNTIL THE CONTRIBUTION IN BOTH DIRECTIONS IS LESS 
  920. C        THAN ACU.
  921. C
  922. C        NMAX -- MAXIMUM NUMBER OF TERMS SUMMED
  923. C
  924.       DATA ACU / 1.0D-10/
  925. C
  926.       XLAMH = XLAM / 2.0D0
  927.       NSPLIT = XLAMH + 0.499999D0
  928.       IF (NSPLIT .LT. 1) NSPLIT = 1
  929.    
  930.       FNCDFC = 1.0
  931.       ESQ = (DFN * FQUAN) / (DFD + DFN * FQUAN)
  932.  
  933.       XI = NSPLIT
  934.       DO 100 I = 1, NSPLIT
  935.       XI = XI - 1.0
  936.       TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
  937.       IF (TERM .LT. -80.0D0) GOTO 200
  938.       TERM = DEXP(TERM)
  939.       TERM = TERM * BETAIN (ESQ, DFN / 2.0D0 + XI, DFD / 2.0D0)
  940.       FNCDFC = FNCDFC - TERM
  941.       IF (TERM .LT. ACU) GOTO 200
  942.   100 CONTINUE
  943.  
  944.   200 XI = NSPLIT - 1
  945.  
  946.   300 XI = XI + 1.0D0
  947.       TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
  948.       IF (TERM .LT. -80.0D0) RETURN
  949.       TERM = DEXP(TERM)
  950.       TERM = TERM * BETAIN (ESQ, DFN / 2.0D0 + XI, DFD / 2.0D0)
  951.       FNCDFC = FNCDFC - TERM
  952.       IF (TERM .LT. ACU) RETURN
  953.       GOTO 300
  954.       END
  955. C
  956.       DOUBLE PRECISION FUNCTION GAMAIN(X,P)
  957.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  958. C
  959. C        ALGORITHM AS 32 J.R.STATIST.SOC. C. (1970) VOL.19 NO.3
  960. C
  961. C        COMPUTES INCOMPLETE GAMMA RATIO FOR POSITIVE VALUES OF
  962. C        ARGUMENTS X AND P.  G MUST BE SUPPLIED AND SHOULD BE EQUAL TO
  963. C        LN(GAMMA(P)).
  964. C        IFAULT = 1 IF P.LE.0 ELSE 2 IF X.LT.0 ELSE 0
  965. C        USES SERIES EXPANSION IF P.GT.X OR X.LE.1, OTHERWISE A 
  966. C        CONTINUED FRACTION APPROXIMATION.C
  967. C
  968. C        [FAULT INDICATOR REMOVED BY G.E.D
  969. C        CALCULATION OF G INSERTED BY G.E.D.]
  970. C
  971.       DIMENSION PN(6)
  972. C
  973. C        DEFINE ACCURACY AND INITIALIZE
  974. C
  975.       ACU = 1.0D-8
  976.       OFLO = 1.0D30
  977.       GIN = 0.0
  978. C
  979. C        TEST FOR ADMISSIBILITY OF ARGUMENTS
  980. C
  981.       IF (P .LE. 0.0D0) STOP 21
  982.       IF (X .LT. 0.0D0) STOP 22
  983.       IF (X .EQ. 0.0D0) GOTO 50
  984.       G = ALGAMA(P)
  985.       FACTOR = DEXP(P * DLOG(X) - X - G)
  986.       IF (X .GT. 1.0D0 .AND. X .GE. P) GOTO 30
  987. C
  988. C        CALCULATION BY SERIES EXPANSION
  989. C
  990.       GIN = 1.0
  991.       TERM = 1.0
  992.       RN = P
  993.    20 RN = RN + 1.0D0
  994.       TERM = TERM * X / RN
  995.       GIN = GIN + TERM
  996.       IF (TERM .GT. ACU) GOTO 20
  997.       GIN = GIN * FACTOR / P
  998.       GOTO 50
  999. C
  1000. C        CALCULATION BY CONTINUED FRACTION
  1001. C
  1002.    30 A = 1.0D0 - P
  1003.       B = A + X + 1.0D0
  1004.       TERM = 0.0
  1005.       PN(1) = 1.0D0
  1006.       PN(2) = X
  1007.       PN(3) = X + 1.0D0
  1008.       PN(4) = X * B
  1009.       GIN = PN(3) / PN(4)
  1010.    32 A = A + 1.0D0
  1011.       B = B + 2.0D0
  1012.       TERM = TERM + 1.0D0
  1013.       AN = A * TERM
  1014.       DO 33 I = 1, 2
  1015.    33 PN(I+4) = B * PN(I+2) - AN * PN(I)
  1016.       IF (PN(6) .EQ. 0.0D0) GOTO 35
  1017.       RN = PN(5) / PN(6)
  1018.       DIF = ABS(GIN - RN)
  1019.       IF (DIF .GT. ACU) GOTO 34
  1020.       IF (DIF .LE. ACU * RN) GOTO 42
  1021.    34 GIN = RN
  1022.    35 DO 36 I = 1, 4
  1023.    36 PN(I) = PN(I+2)
  1024.       IF (ABS(PN(5)).LT.OFLO) GOTO 32
  1025.       DO 41 I = 1, 4
  1026.    41 PN(I) = PN(I) / OFLO
  1027.       GOTO 32
  1028.    42 GIN = 1.0D0 - FACTOR * GIN
  1029. C
  1030.    50 GAMAIN = GIN
  1031.       RETURN
  1032.       END
  1033. C
  1034.       DOUBLE PRECISION FUNCTION GAUINV(P,IFAULT)
  1035.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  1036. C
  1037. C              ALGORITHM AS 70 APPL. STATIST. (1974) VOL. 23, NO.1
  1038. C              GAUINV FINDS PERCENTAGE POINTS OF THE NORMAL DISTRIBUTION
  1039. C
  1040.       DATA ZERO,ONE,HALF,ALIMIT/0.0D0, 1.0D0, 0.5D0, 1.0D-20/
  1041. C
  1042.       DATA               P0,    P1,               P2,                P3
  1043.      1   / -.322232431088D0, -1.D0, -.342242088547D0, -.204231210245D-1/
  1044. C
  1045.       DATA                P4,               Q0,              Q1
  1046.      1   / -.453642210148D-4, .993484626060D-1, .588581570495D0/
  1047. C
  1048.       DATA              Q2,              Q3,              Q4
  1049.      1   / .531103462366D0, .103537752850D0, .38560700634D-2/
  1050. C
  1051.       IFAULT=1
  1052.       GAUINV=ZERO
  1053.       PS=P
  1054.       IF (PS.GT.HALF) PS=ONE-PS
  1055.       IF (PS.LT.ALIMIT) RETURN
  1056.       IFAULT=0
  1057.       IF (PS.EQ.HALF) RETURN
  1058.       YI=DSQRT(DLOG(ONE/(PS*PS)))
  1059.       GAUINV=YI+((((YI*P4+P3)*YI+P2)*YI+P1)*YI+P0)
  1060.      1         /((((YI*Q4+Q3)*YI+Q2)*YI+Q1)*YI+Q0)
  1061.       IF (P.LT.HALF) GAUINV=-GAUINV
  1062.       RETURN
  1063.       END
  1064. C
  1065.       DOUBLE PRECISION FUNCTION PPCHI2(P,V)
  1066.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  1067. C
  1068. C        ALGORITHM AS 91  APPL. STATIST. (1975) VOL.24, NO.3
  1069. C
  1070. C        TO EVALUATE THE PECENTAGE POINTS OF THE CHI-SQUARED
  1071. C        PROBABILITY DISTRIBUTION FUNCTION.
  1072. C        P MUST LIE IN THE RANGE 0.000002 TO 0.999998, V MUST BE POSITIVE,
  1073. C        G MUST BE SUPPLIED AND SHOULD BE EQUAL TO LN(GAMMA(V/2.0))
  1074. C
  1075. C        [FAULT INDICATOR REMOVED BY G.E.D.  
  1076. C        CALCULATION OF G ADDED BY G.E.D.]
  1077. C
  1078.       DATA E, AA /0.5D-6, 0.6931471805D0/
  1079. C
  1080. C        AFTER DEFINING THE ACCURACY AND LN(2), TEST ARGUMENTS AND INITIALIZE
  1081. C
  1082.       PPCHI2 = -1.0
  1083.       IF (P .LT. 0.000002D0 .OR. P .GT. 0.999998D0) STOP 31
  1084.       IF (V .LE. 0.0D0) STOP 32
  1085.       G = ALGAMA(V / 2.0D0)
  1086.       XX = 0.5D0 * V
  1087.       C = XX - 1.0D0
  1088. C
  1089. C        STARTING APPROXIMATION FOR SMALL CHI-SQUARED
  1090. C
  1091.       IF (V .GE. -1.24D0 * DLOG(P)) GOTO 1
  1092.       CH = (P * XX * DEXP(G + XX * AA)) ** (1.0D0 / XX)
  1093.       IF (CH - E) 6, 4, 4
  1094. C
  1095. C        STARTING APPROXIMATION FOR V LESS THAN OR EQUAL TO 0.32
  1096. C
  1097.     1 IF (V .GT. 0.32D0) GOTO 3
  1098.       CH = 0.4
  1099.       A = DLOG(1.0D0 - P)
  1100.     2 Q = CH
  1101.       P1 = 1.0D0 + CH * (4.67D0 + CH)
  1102.       P2 = CH * (6.73D0 + CH * (6.66D0 + CH))
  1103.       T = -0.5D0 + (4.67D0 + 2.0D0 * CH) / P1 -
  1104.      *   (6.73D0 + CH * (13.32D0 + 3.0D0 * CH)) / P2
  1105.       CH = CH - (1.0D0 - DEXP(A + G + 0.5D0 * CH + C * AA) * P2 / P1)/T
  1106.       IF (ABS(Q / CH - 1.0D0) - 0.01D0) 4, 4, 2
  1107. C
  1108. C        CALL TO ALGORITHM AS 70 - NOTE THAT P HAS BEEN TESTED ABOVE
  1109. C
  1110.     3 X = GAUINV(P,IF1)
  1111. C
  1112. C        STARTING APPROXIMATION USING WILSON AND HILFERTY ESTIMATE
  1113. C
  1114.       P1 = (2.0D0 / 9.0D0) / V
  1115.       CH = V * (X * DSQRT(P1) + 1.0D0 - P1) ** 3
  1116. C
  1117. C        STARTING APPROXIMATION FOR P TENDING TO 1
  1118. C
  1119.       IF (CH .GT. 2.2D0 * V + 6.0D0)
  1120.      *   CH = -2.0D0 * (DLOG(1.0D0 - P) - C * DLOG(0.5D0 * CH) + G)
  1121. C
  1122. C        CALL TO ALGORITHM AS 32 AND CALCULATION OF SEVEN 
  1123. C        TERM TAYLOR SERIES
  1124. C
  1125.     4 Q = CH
  1126.       P1 = 0.5D0 * CH
  1127.       P2 = P - GAMAIN(P1, XX)
  1128.       T = P2 * DEXP(XX * AA + G + P1 - C * DLOG(CH))
  1129.       B = T / CH
  1130.       A = 0.5D0 * T - B * C
  1131.       S1 = (210.0D0+A*(140.0D0+A*(105.0D0+A*
  1132.      *   (84.0D0+A*(70.0D0+60.0D0*A))))) / 420.0D0
  1133.       S2 = (420.0D0+A*(735.0D0+A*(966.0D0+A*(1141.0D0+1278.0D0
  1134.      *   *A)))) / 2520.0D0
  1135.       S3 = (210.0D0 + A * (426.0D0 + A * (707.0D0 + 932.0D0 
  1136.      *   * A))) / 2520.0D0
  1137.       S4 =(252.0D0+A*(672.0D0+1182.0D0*A)+C*(294.0D0+A*
  1138.      *   (889.0D0+1740.0D0*A)))/5040.0D0
  1139.       S5 = (84.0D0 + 264.0D0 * A + C * (175.0D0 
  1140.      *   + 606.0D0 * A)) / 2520.0D0
  1141.       S6 = (120.0D0 + C * (346.0D0 + 127.0D0 * C)) / 5040.0D0
  1142.       CH = CH+T*(1.0+0.5D0*T*S1-B*C*(S1-B*(S2-B
  1143.      *   *(S3-B*(S4-B*(S5-B*S6))))))
  1144.       IF (ABS(Q / CH - 1.0D0) .GT. E) GOTO 4
  1145. C
  1146.     6 PPCHI2 = CH
  1147.       RETURN
  1148.       END
  1149.  
  1150.       DOUBLE PRECISION FUNCTION RCRIT(ALPHA,F2)
  1151.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  1152.       F1 = 1.0
  1153.       AA = F2 / 2.0D0
  1154.       BB = F1 / 2.0D0
  1155.       XALPHA = XINBTA(AA, BB, ALPHA)
  1156.       FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
  1157.       RCRIT = DSQRT(FALPHA / (FALPHA + F2))
  1158.       RETURN
  1159.       END
  1160.  
  1161.       DOUBLE PRECISION FUNCTION RDF(NOBS, R, RHO, IFAULT)
  1162. C
  1163. C        THE INTEGRAL OF THE DISTRIBUTION FUNCTION OF THE CORRELATION 
  1164. C        COEFFICIENT IN A SAMPLE OF SIZE 'N' AND POPULATION VALUE 'RHO'
  1165. C        FROM -1 TO 'R'
  1166. C
  1167. C        INCORPORATES A BINARY SEARCH TO SHORTEN THE LIMITS OF NUMERICAL
  1168. C        INTEGRATION
  1169. C
  1170.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  1171.       DATA PI /3.141592653589793238D0/, NBIN /4/, SMALL /1.0D-10/
  1172. C
  1173. C        CHECK ARGUMENTS
  1174. C
  1175.       IFAULT = 0
  1176.       IF (NOBS.LE.2) IFAULT = 1
  1177.       IF (DABS(R).GT.1.0D0) IFAULT = 2
  1178.       IF (DABS(RHO).GE.1.0D0) IFAULT = 3
  1179.       IF (IFAULT.GT.0) RETURN
  1180. C
  1181. C        PATHOLOGICAL CASES
  1182. C
  1183.       RDF = (R + 1.0D0) / 2.0D0
  1184.       IF (DABS(R).EQ.1.0D0) RETURN
  1185. C
  1186.       IF (NOBS - 4) 10, 20, 30
  1187. C
  1188. C         THREE OBSERVATIONS
  1189. C
  1190.    10 RDF = (DACOS(-R) - RHO * DSQRT((1.0D0 - R**2) / 
  1191.      *   (1.0D0 - RHO**2 * R**2)) * DACOS(-RHO * R)) / PI
  1192.       RETURN
  1193. C  
  1194. C         FOUR OBSERVATIONS
  1195. C
  1196.    20 RDF = (DACOS(RHO) + R * DACOS(-RHO * R) *
  1197.      *   ((1.0D0 - RHO**2) / (1.0D0 - RHO**2 * R**2))**1.5) / PI
  1198.       IF (RHO.NE.0.0D0) RDF = RDF + ((1.0D0 - RHO**2)**1.5 /
  1199.      *   (1.0D0 - RHO**2 * R**2) - DSQRT(1.0D0 - RHO**2)) /
  1200.      *   (PI * RHO)
  1201.       RETURN
  1202. C
  1203. C         FIVE OR MORE OBSERVATIONS
  1204. C
  1205.    30 IF (R.GT.RHO) GOTO 100
  1206.       RDF = 0.0
  1207.       IF (RORD(NOBS, R, RHO).LE.SMALL) RETURN
  1208. C
  1209. C         LOWER TAIL
  1210. C
  1211. C         BINARY SEARCH
  1212. C
  1213.       XLO = -1.0
  1214.       XHI = R
  1215. C
  1216.       DO 40 I = 1, NBIN
  1217.       XM = (XLO + XHI) / 2.0D0
  1218.       YM = RORD(NOBS, XM, RHO)
  1219.       IF (YM.GT.SMALL) XHI = XM
  1220.       IF (YM.LE.SMALL) XLO = XM
  1221.    40 CONTINUE
  1222. C
  1223.       RDF = SIMPSN(XLO,R,NOBS,RHO)
  1224.       RETURN
  1225. C
  1226. C         UPPER TAIL
  1227. C
  1228. C         BINARY SEARCH
  1229. C
  1230.   100 RDF = 1.0
  1231.       IF (RORD(NOBS, R, RHO).LE.SMALL) RETURN
  1232.       XLO = R
  1233.       XHI = 1.0
  1234. C
  1235.       DO 140 I = 1, NBIN
  1236.       XM = (XLO + XHI) / 2.0D0
  1237.       YM = RORD(NOBS, XM, RHO)
  1238.       IF (YM.GT.SMALL) XLO = XM
  1239.       IF (YM.LE.SMALL) XHI = XM
  1240.   140 CONTINUE
  1241. C
  1242.       RDF = 1.0D0 - SIMPSN(R,XHI,NOBS,RHO)
  1243.       RETURN
  1244.       END
  1245.  
  1246.       DOUBLE PRECISION FUNCTION RORD(NOBS, R, RHO)
  1247.  
  1248. C
  1249. C        ORDINATE OF THE DISTRIBUTION OF THE CORRELATION COEFFICIENT
  1250. C        IN A SAMPLE OF SIZE  'N' AND POPULATION VALUE 'RHO' AT 'R'
  1251. C
  1252. C        F.N. DAVID RECURRENCE FORMULA
  1253. C
  1254.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  1255.       DATA PI /3.141592653589793238D0/
  1256. C
  1257.       RORD = 0.0
  1258.       IF (R.EQ.-1.0D0 .OR. R.EQ.1.0D0) RETURN
  1259. C
  1260.       R2 = R ** 2
  1261.       RHO2 = RHO ** 2
  1262.       RRHO = R * RHO
  1263.       R2RHO2 = R2 * RHO2
  1264.       X1 = RRHO * DSQRT((1.0D0 - RHO2) * (1.0D0 - R2)) /
  1265.      *   (1.0D0 - R2RHO2)
  1266.       X2 = (1.0D0 - RHO2) * (1.0D0 - R2) / (1.0D0 - R2RHO2)
  1267.  
  1268. C
  1269.       YNM2 = (1.0D0 - RHO2) / (PI * DSQRT(1.0D0 - R2)) *
  1270.      *   (1.0D0 / (1.0D0 - R2RHO2) + RRHO * DACOS(-RRHO) /
  1271.      *   (1.0D0 - R2RHO2)**1.5)
  1272.       YNM1 = (1.0D0 - RHO2)**1.5 / PI * (3.0D0 * RRHO /
  1273.      *   (1.0D0 - R2RHO2)**2 + (1.0D0 + 2.0D0 * R2RHO2) * 
  1274.      *   DACOS(-RRHO) / (1.0D0 - R2RHO2)**2.5)
  1275. C
  1276.       DO 100 I = 5, NOBS
  1277.       XI = I
  1278.       YN = (2.0D0 * XI - 5.0D0) / (XI - 3.0D0) * X1 * YNM1
  1279.      *   + (XI - 3.0D0) / (XI - 4.0D0) * X2 * YNM2
  1280.       YNM2 = YNM1
  1281.       YNM1 = YN
  1282.   100 CONTINUE
  1283. C
  1284.       RORD = YN
  1285.       RETURN
  1286.       END
  1287.  
  1288.       DOUBLE PRECISION FUNCTION SIMPSN(XLO,XHI,NOBS,RHO)
  1289. C
  1290. C        SIMPSON'S RULE... THE INTEGRAL OF F(X) FROM
  1291. C           X = XLO TO X = XHI
  1292. C
  1293.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  1294.       DATA ERR /1.0D-6/, RERR /1.0D-4/
  1295. C
  1296. C        AT EACH STEP ADD NEW POINTS BETWEEN THE OLD POINTS,
  1297. C        WHEN A SET OF POINTS FIRST ENTERS IT COMES IN WITH
  1298. C           A COEFFICIENT OF 4 .  AT THE NEXT AND ALL 
  1299. C           SUBSEQUENT STEPS IT HAS A COEFFICIENT OF 2.
  1300. C
  1301.       LOOP = -1
  1302.       RANGE = XHI - XLO
  1303.       H = RANGE
  1304.       SUM = RORD(NOBS,XLO,RHO) + RORD(NOBS,XHI,RHO)
  1305.       SIMPSN = 0.0
  1306.       IF (SUM.LT.ERR) RETURN
  1307.       SIMP0 = H / 3.0D0 * SUM
  1308. C
  1309.    10 LOOP = LOOP + 1
  1310.       NSTOP = 2 ** LOOP
  1311.       X2STOP = 2 * NSTOP
  1312.       H = H / 2.0D0
  1313.       TERM = 0.0
  1314. C
  1315.       DO 100 I = 1, NSTOP
  1316.       XORD = XLO + DBLE(2 * I - 1) / X2STOP * RANGE
  1317.       TERM = TERM + RORD(NOBS,XORD,RHO)
  1318.   100 CONTINUE
  1319. C
  1320.       SUM = SUM + 4.0D0 * TERM
  1321.       SIMPSN = H / 3.0D0 * SUM
  1322.       DIFF = DABS(SIMPSN - SIMP0)
  1323.       IF (DIFF.LE.ERR .AND. DIFF/SIMPSN.LE.RERR) RETURN
  1324. C
  1325.       SIMP0 = SIMPSN
  1326.       SUM = SUM - 2.0D0 * TERM
  1327.       GOTO 10
  1328.       END
  1329.  
  1330. C
  1331.       DOUBLE PRECISION FUNCTION XINBTA(P, Q, ALPHA)
  1332.       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  1333. C
  1334. C        ALGORITHM AS 109 APPL.STATIST. (1977), VOL.26, NO.1
  1335. C        (REPLACING ALGORITHM AS 64  APPL. STATIST. (1973),
  1336. C        VOL.22, NO.3)
  1337. C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
  1338. C
  1339. C        COMPUTES INVERSE OF INCOMPLETE BETA FUNCTION
  1340. C        RATIO FOR GIVEN POSITIVE VALUES OF THE ARGUMENTS
  1341. C        P AND Q, ALPHA BETWEEN ZERO AND ONE.
  1342. C        LOG OF COMPLETE BETA FUNCTION, BETA, IS ASSUMED TO BE KNOWN.
  1343. C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
  1344. C
  1345.       LOGICAL INDEX
  1346. C
  1347. C        DEFINE ACCURACY AND INITIALIZE
  1348. C
  1349.       DATA ACU /1.0D-14/
  1350.       XINBTA = ALPHA
  1351. C
  1352. C        TEST FOR ADMISIBILITY OF PARAMETERS
  1353. C
  1354.       IF (P.LE.0.0D0.OR.Q.LE.0.0D0) STOP 50
  1355.       IF (ALPHA.LT.0.0D0.OR.ALPHA.GT.1.0D0) STOP 51
  1356.       IF (ALPHA.EQ.0.0D0.OR.ALPHA.EQ.1.0D0) RETURN
  1357. C
  1358.       BETA=ALGAMA(P)+ALGAMA(Q)-ALGAMA(P+Q)
  1359. C
  1360. C        CHANGE TAIL IF NECESSARY
  1361. C
  1362.       IF (ALPHA .LE. 0.5D0) GOTO 1
  1363.       A = 1.0D0 - ALPHA
  1364.       PP = Q
  1365.       QQ = P
  1366.       INDEX = .TRUE.
  1367.       GOTO 2
  1368.     1 A = ALPHA
  1369.       PP = P 
  1370.       QQ = Q
  1371.       INDEX = .FALSE.
  1372. C
  1373. C        CALCULATE INITIAL APPROXIMATION
  1374. C
  1375.     2 R = DSQRT(-DLOG(A * A))
  1376.       Y = R - (2.30753D0 + 0.27061D0 * R) / 
  1377.      *             (1.0D0 + (0.99229D0 + 0.04481D0 * R) * R)
  1378.       IF (PP .GT. 1.0D0 .AND. QQ .GT. 1.0D0) GOTO 5
  1379.       R = QQ + QQ
  1380.       T = 1.0D0 / (9.0D0 * QQ)
  1381.       T = R * (1.0D0 - T + Y * DSQRT(T)) ** 3
  1382.       IF (T .LE. 0.0D0) GOTO 3
  1383.       T = (4.0D0 * PP + R - 2.0D0) / T
  1384.       IF (T .LE. 1.0D0) GOTO 4
  1385.       XINBTA = 1.0D0 - 2.0D0 / (T + 1.0D0)
  1386.       GOTO 6
  1387.     3 XINBTA = 1.0D0 - DEXP((DLOG((1.0D0 - A) * QQ) + BETA) / QQ)
  1388.       GOTO 6
  1389.     4 XINBTA = DEXP((DLOG(A * PP) + BETA) / PP)
  1390.       GOTO 6
  1391.     5 R = (Y * Y - 3.0D0) / 6.0D0
  1392.       S = 1.0D0 / (PP + PP - 1.0D0)
  1393.       T = 1.0D0 / (QQ + QQ - 1.0D0)
  1394.       H = 2.0D0 / (S + T)
  1395.       W = Y * DSQRT(H + R) / H - (T - S) * 
  1396.      *   (R + 5.0D0 / 6.0D0 - 2.0D0 / (3.0D0 * H))
  1397.       XINBTA = PP / (PP + QQ * DEXP(W + W))
  1398. C
  1399. C        SOLVE FOR X BY A MODIFIED NEWTON-RAPHSON METHOD,
  1400. C        USING THE FUNCTION BETAIN.
  1401. C
  1402.     6 R = 1.0D0 - PP
  1403.       T = 1.0D0 - QQ
  1404.       YPREV = 0.0
  1405.       SQ = 1.0
  1406.       PREV = 1.0
  1407.       IF (XINBTA .LT. 0.0001D0) XINBTA = 0.0001D0
  1408.       IF (XINBTA .GT. 0.9999D0) XINBTA = 0.9999D0
  1409.     7 Y = BETAIN(XINBTA, PP, QQ)
  1410.       Y = (Y - A) * DEXP(BETA +
  1411.      *             R * DLOG(XINBTA) + T * DLOG(1.0D0 - XINBTA))
  1412.       IF (Y * YPREV .LE. 0.0D0) PREV = SQ
  1413.       G = 1.0
  1414.     9 ADJ = G * Y
  1415.       SQ = ADJ * ADJ
  1416.       IF (SQ .GE. PREV) GOTO 10
  1417.       TX = XINBTA - ADJ
  1418.       IF (TX .GE. 0.0D0 .AND. TX .LE. 1.0D0) GOTO 11
  1419.    10 G = G / 3.0D0
  1420.       GOTO 9
  1421.    11 IF (PREV .LE. ACU) GOTO 12
  1422.       IF (Y * Y .LE. ACU) GOTO 12
  1423.       IF (TX .EQ. 0.0D0 .OR. TX .EQ. 1.0D0) GOTO 10
  1424.       IF (TX .EQ. XINBTA) GOTO 12
  1425.       XINBTA = TX
  1426.       YPREV = Y
  1427.       GOTO 7
  1428.    12 IF (INDEX) XINBTA = 1.0D0 - XINBTA
  1429.       RETURN
  1430.       END
  1431.