SUBROUTINE BETRAN(N,ALPHA,BETA,ISEED,X)
C     ***** STILL NEEDS ALGORITHM WORK ******
C
C     PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N
C              FROM THE BETA DISTRIBUTION
C          WITH SINGLE PRECISION SHAPE
C          PARAMETERS = ALPHA AND BETA.
C              THE PROTOTYPE BETA DISTRIBUTION USED
C              HEREIN HAS MEAN = ALPHA/(ALPHA+BETA)
C              AND STANDARD DEVIATION =
C              SQRT((ALPHA*BETA) / ((ALPHA+BETA)**2)*(ALPHA+BETA+1))
C              THIS DISTRIBUTION IS DEFINED FOR ALL X
C              BETWEEN 0.0 (INCLUSIVELY) AND 1.0 (INCLUSIVELY).
C              AND HAS THE PROBABILITY DENSITY FUNCTION
C              F(X) = (1/CONSTANT) * X**(ALPHA-1) * (1.0-X)**(BETA-1)
C              WHERE THE CONSTANT = THE BETA FUNCTION EVALUATED
C              AT THE VALUES ALPHA AND BETA.
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER
C                                OF RANDOM NUMBERS TO BE
C                                GENERATED.
C                     --ALPHA  = THE SINGLE PRECISION VALUE OF THE
C                                FIRST  SHAPE PARAMETER.
C                                ALPHA SHOULD BE GREATER THAN
C                                OR EQUAL TO 1.0.
C                     --BETA   = THE SINGLE PRECISION VALUE OF THE
C                                SECOND SHAPE PARAMETER.
C                                BETA  SHOULD BE GREATER THAN
C                                OR EQUAL TO 1.0.
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                RANDOM SAMPLE WILL BE PLACED.
C     OUTPUT--A RANDOM SAMPLE OF SIZE N
C             FROM THE BETA DISTRIBUTION
C             WITH SHAPE PARAMETER VALUES = ALPHA AND BETA.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C                 --ALPHA SHOULD BE GREATER THAN
C                   OR EQUAL TO 1.0.
C                 --BETA  SHOULD BE GREATER THAN
C                   OR EQUAL TO 1.0.
C     OTHER DATAPAC   SUBROUTINES NEEDED--UNIRAN, NORRAN.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--GREENWOOD, 'A FAST GENERATOR FOR
C                 BETA-DISTRIBUTED RANDOM VARIABLES',
C                 COMPSTAT 1974, PROCEEDINGS IN
C                 COMPUTATIONAL STATISTICS, VIENNA,
C                 SEPTEMBER, 1974, PAGES 19-27.
C               --TOCHER, THE ART OF SIMULATION,
C                 1963, PAGES 24-27.
C               --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS,
C                 1964, PAGES 36-37.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGES 37-56.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 30-35.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 952.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-3651
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL BUREAU OF STANDARDS.
C           THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED,
C           MODIFIED, OR OTHERWISE USED IN A CONTEXT
C           OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM.
C     LANGUAGE--ANSI FORTRAN (1966)
C               EXCEPTION--HOLLARITH STRINGS IN FORMAT STATEMENTS
C                          DENOTED BY QUOTES RATHER THAN NH.
C     VERSION NUMBER--82.3
C     ORIGINAL VERSION--NOVEMBER  1975.
C     UPDATED         --FEBRUARY  1976.
C     UPDATED         --JUNE      1978.
C     UPDATED         --DECEMBER  1981.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
C
      DIMENSION U(10)
C
C---------------------------------------------------------------------
C
CCCCC CHARACTER*4 IFEEDB
CCCCC CHARACTER*4 IPRINT
C
CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
CCCCC COMMON /PRINT/IFEEDB,IPRINT
C
C-----DATA STATEMENTS-------------------------------------------------
C
      DATA ATHIRD/0.33333333/
      DATA SQRT3 /1.73205081/
C
      IPR=6
C
C
C-----START POINT-----------------------------------------------------
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(ALPHA.LT.1.0)GOTO60
      IF(BETA.LT.1.0)GOTO65
      GOTO90
   50 WRITE(IPR, 5)
      WRITE(IPR,47)N
      RETURN
   60 WRITE(IPR,16)
      WRITE(IPR,46)ALPHA
      RETURN
   65 WRITE(IPR,26)
      WRITE(IPR,46)BETA
      RETURN
   90 CONTINUE
    5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 BETRAN SUBROUTINE IS NON-POSITIVE *****)
   16 FORMAT(1H , 95H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 BETRAN SUBROUTINE IS SMALLER THAN 1.0 *****)
   26 FORMAT(1H , 95H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 BETRAN SUBROUTINE IS SMALLER THAN 1.0 *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C     GENERATE N BETA RANDOM NUMBERS
C     BY USING THE FACT THAT
C     IF X1 IS A GAMMA VARIATE WITH PARAMETER ALPHA
C     AND IF X2 IS A GAMMA VARIATE WITH PARAMETER BETA,
C     THEN THE RATIO X1/(X1+X2) IS A BETA VARIATE
C     WITH PARAMETERS ALPHA AND BETA.
C
C     TO GENERATE N GAMMA DISTRIBUTION RANDOM NUMBERS,
C     USE GREENWOOD'S REJECTION ALGORITHM--
C     1) GENERATE A NORMAL RANDOM NUMBER;
C     2) TRANSFORM THE NORMAL VARIATE TO AN APPROXIMATE
C        GAMMA VARIATE USING THE WILSON-HILFERTY
C        APPROXIMATION (SEE THE JOHNSON AND KOTZ
C        REFERENCE, PAGE 176);
C     3) FORM THE REJECTION FUNCTION VALUE, BASED
C        ON THE PROBABILITY DENSITY FUNCTION VALUE
C        OF THE ACTUAL DISTRIBUTION OF THE PSEUDO-GAMMA
C        VARIATE, AND THE PROBABILITY DENSITY FUNCTION VALUE
C        OF A TRUE GAMMA VARIATE.
C     4) GENERATE A UNIFORM RANDOM NUMBER;
C     5) IF THE UNIFORM RANDOM NUMBER IS LESS THAN
C        THE REJECTION FUNCTION VALUE, THEN ACCEPT
C        THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE;
C        IF THE UNIFORM RANDOM NUMBER IS LARGER THAN
C        THE REJECTION FUNCTION VALUE, THEN REJECT
C        THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE.
C
      A1=1.0/(9.0*ALPHA)
      B1=SQRT(A1)
      XN01=-SQRT3+B1
      XG01=ALPHA*(1.0-A1+B1*XN01)**3
      A2=1.0/(9.0*BETA)
      B2=SQRT(A2)
      XN02=-SQRT3+B2
      XG02=BETA*(1.0-A2+B2*XN02)**3
C
      DO100I=1,N
C
  150 CALL NORRAN(1,ISEED,XN)
      XG=ALPHA*(1.0-A1+B1*XN)**3
      IF(XG.LT.0.0)GOTO150
      TERM=(XG/XG01)**(ALPHA-ATHIRD)
      ARG=0.5*XN*XN-XG-0.5*XN01*XN01+XG01
      FUNCT=TERM*EXP(ARG)
      CALL UNIRAN(1,ISEED,U)
      IF(U(1).LE.FUNCT)GOTO170
      GOTO150
  170 XG1=XG
C
  250 CALL NORRAN(1,ISEED,XN)
      XG=BETA*(1.0-A2+B2*XN)**3
      IF(XG.LT.0.0)GOTO250
      TERM=(XG/XG02)**(BETA-ATHIRD)
      ARG=0.5*XN*XN-XG-0.5*XN02*XN02+XG02
      FUNCT=TERM*EXP(ARG)
      CALL UNIRAN(1,ISEED,U)
      IF(U(1).LE.FUNCT)GOTO270
      GOTO250
  270 XG2=XG
C
      X(I)=XG1/(XG1+XG2)
C
  100 CONTINUE
C
      RETURN
      END