File BDTR.FT (FORTRAN source file)

Directory of image this file is from
This file as a plain text file
C
C     ..................................................................
C
C        SUBROUTINE BDTR
C
C        PURPOSE
C           COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,
C           DISTRIBUTED ACCORDING TO THE BETA DISTRIBUTION WITH
C           PARAMETERS A AND B, IS LESS THAN OR EQUAL TO X.  F(A,B,X),
C           THE ORDINATE OF THE BETA DENSITY AT X, IS ALSO COMPUTED.
C
C        USAGE
C           CALL BDTR(X,A,B,P,D,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X   - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.
C           A   - BETA DISTRIBUTION PARAMETER (CONTINUOUS).
C           B   - BETA DISTRIBUTION PARAMETER (CONTINUOUS).
C           P   - OUTPUT PROBABILITY.
C           D   - OUTPUT DENSITY.
C           IER - RESULTANT ERROR CODE WHERE
C               IER= 0 --- NO ERROR
C               IER=-1,+1  CDTR HAS BEEN CALLED AND AN ERROR HAS
C                          OCCURRED.  SEE CDTR.
C               IER=-2 --- AN INPUT PARAMETER IS INVALID.  X IS LESS
C                          THAN 0.0 OR GREATER THAN 1.0, OR EITHER A OR
C                          B IS LESS THAN 0.5 OR GREATER THAN 10**(+5).
C                          P AND D ARE SET TO -1.7E38.                          0
C               IER=+2 --- INVALID OUTPUT.  P IS LESS THAN ZERO OR
C                          GREATER THAN ONE.  P IS SET TO 1.7E38.               0
C
C        REMARKS
C           SEE MATHEMATICAL DESCRIPTION.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           DLGAM
C           NDTR
C           CDTR
C
C        METHOD
C           REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL
C           DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,
C           IBM RESEARCH REPORT RC-1094, 1963.
C
C     ..................................................................
C
      SUBROUTINE BDTR(X,A,B,P,D,IER)
      DOUBLE PRECISION XX,DLXX,DL1X,AA,BB,G1,G2,G3,G4,DD,PP,XO,FF,FN,
     1XI,SS,CC,RR,DLBETA
C
C        TEST FOR VALID INPUT DATA
C
      IF(A-(.5-1.E-5)) 640,10,10
   10 IF(B-(.5-1.E-5)) 640,20,20
   20 IF(A-1.E+5) 30,30,640
   30 IF(B-1.E+5) 40,40,640
   40 IF(X) 640,50,50
   50 IF(1.-X) 640,60,60
C
C        COMPUTE LOG(BETA(A,B))
C
   60 AA=DBLE(A)
      BB=DBLE(B)
      CALL DLGAM(AA,G1,IOK)
      CALL DLGAM(BB,G2,IOK)
      CALL DLGAM(AA+BB,G3,IOK)
      DLBETA=G1+G2-G3
C
C        TEST FOR X NEAR 0.0 OR 1.0
C
      IF(X-1.E-8) 80,80,70
   70 IF((1.-X)-1.E-8) 130,130,140
   80 P=0.0
      IF(A-1.) 90,100,120
   90 D=1.7E38
      GO TO 660
  100 DD=-DLBETA
      IF(DD+1.68D02)  120,120,110
  110 DD=DEXP(DD)
      D=SNGL(DD)
      GO TO 660
  120 D=0.0
      GO TO 660
  130 P=1.0
      IF(B-1.) 90,100,120
C
C        SET PROGRAM PARAMETERS
C
  140 XX=DBLE(X)
      DLXX=DLOG(XX)
      DL1X=DLOG(1.D0-XX)
      XO=XX/(1.D0-XX)
      ID=0
C
C        COMPUTE ORDINATE
C
      DD=(AA-1.D0)*DLXX+(BB-1.D0)*DL1X-DLBETA
      IF(DD-1.68D02) 150,150,160
  150 IF(DD+1.68D02) 170,170,180
  160 D=1.7E38                                                                  0
      GO TO 190
  170 D=0.0
      GO TO 190
  180 DD=DEXP(DD)
      D=SNGL(DD)
C
C        A OR B OR BOTH WITHIN 1.E-8 OF 1.0
C
  190 IF(ABS(A-1.)-1.E-8)  200,200,210
  200 IF(ABS(B-1.)-1.E-8)  220,220,230
  210 IF(ABS(B-1.)-1.E-8)  260,260,290
  220 P=X
      GO TO 660
  230 PP=BB*DL1X
      IF(PP+1.68D02) 240,240,250
  240 P=1.0
      GO TO 660
  250 PP=DEXP(PP)
      PP=1.D0-PP
      P=SNGL(PP)
      GO TO 600
  260 PP=AA*DLXX
      IF(PP+1.68D02) 270,270,280
  270 P=0.0
      GO TO 660
  280 PP=DEXP(PP)
      P=SNGL(PP)
      GO TO 600
C
C        TEST FOR A OR B GREATER THAN 1000.0
C
  290 IF(A-1000.) 300,300,310
  300 IF(B-1000.) 330,330,320
  310 XX=2.D0*AA/XO
      XS=SNGL(XX)
      AA=2.D0*BB
      DF=SNGL(AA)
      CALL CDTR(XS,DF,P,DUMMY,IER)
      P=1.0-P
      GO TO 670
  320 XX=2.D0*BB*XO
      XS=SNGL(XX)
      AA=2.D0*AA
      DF=SNGL(AA)
      CALL CDTR(XS,DF,P,DUMMY,IER)
      GO TO 670
C
C        SELECT PARAMETERS FOR CONTINUED FRACTION COMPUTATION
C
  330 IF(X-.5) 340,340,380
  340 IF(AA-1.D0) 350,350,360
  350 RR=AA+1.D0
      GO TO 370
  360 RR=AA
  370 DD=DLXX/5.D0
      DD=DEXP(DD)
      DD=(RR-1.D0)-(RR+BB-1.D0)*XX*DD +2.D0
      IF(DD) 420,420,430
  380 IF(BB-1.D0) 390,390,400
  390 RR=BB+1.D0
      GO TO 410
  400 RR=BB
  410 DD=DL1X/5.D0
      DD=DEXP(DD)
      DD=(RR-1.D0)-(AA+RR-1.D0)*(1.D0-XX)*DD +2.D0
      IF(DD) 430,430,420
  420 ID=1
      FF=DL1X
      DL1X=DLXX
      DLXX=FF
      XO=1.D0/XO
      FF=AA
      AA=BB
      BB=FF
      G2=G1
C
C        TEST FOR A LESS THAN 1.0
C
  430 FF=0.D0
      IF(AA-1.D0) 440,440,470
  440 CALL DLGAM(AA+1.D0,G4,IOK)
      DD=AA*DLXX+BB*DL1X+G3-G2-G4
      IF(DD+1.68D02) 460,460,450
  450 FF=FF+DEXP(DD)
  460 AA=AA+1.D0
C
C        COMPUTE P USING CONTINUED FRACTION EXPANSION
C
  470 FN=AA+BB-1.D0
      RR=AA-1.D0
      II=80
      XI=DFLOAT(II)
      SS=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))
      SS=SS*XO
      DO 480 I=1,79
      II=80-I
      XI=DFLOAT(II)
      DD=(XI*(FN+XI))/((RR+2.D0*XI+1.D0)*(RR+2.D0*XI))
      DD=DD*XO
      CC=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))
      CC=CC*XO
      SS=CC/(1.D0+DD/(1.D0-SS))
  480 CONTINUE
      SS=1.D0/(1.D0-SS)
      IF(SS) 650,650,490
  490 CALL DLGAM(AA+BB,G1,IOK)
      CALL DLGAM(AA+1.D0,G4,IOK)
      CC=G1-G2-G4+AA*DLXX+(BB-1.D0)*DL1X
      PP=CC+DLOG(SS)
      IF(PP+1.68D02) 500,500,510
  500 PP=FF
      GO TO 520
  510 PP=DEXP(PP)+FF
  520 IF(ID) 540,540,530
  530 PP=1.D0-PP
  540 P=SNGL(PP)
C
C        SET ERROR INDICATOR
C
      IF(P) 550,570,570
  550 IF(ABS(P)-1.E-7) 560,560,650
  560 P=0.0
      GO TO 660
  570 IF(1.-P) 580,600,600
  580 IF(ABS(1.-P)-1.E-7) 590,590,650
  590 P=1.0
      GO TO 660
  600 IF(P-1.E-8) 610,610,620
  610 P=0.0
      GO TO 660
  620 IF((1.0-P)-1.E-8) 630,630,660
  630 P=1.0
      GO TO 660
  640 IER=-2
      D=-1.7E38                                                                 0
      P=-1.7E38                                                                 0
      GO TO 670
  650 IER=+2
      P= 1.7E38                                                                 0
      GO TO 670
  660 IER=0
  670 RETURN
      END
Feel free to contact me, David Gesswein djg@pdp8online.com with any questions, comments on the web site, or if you have related equipment, documentation, software etc. you are willing to part with. I am interested in anything PDP-8 related, computers, peripherals used with them, DEC or third party, or documentation.

PDP-8 Home Page PDP-8 Site Map PDP-8 Site Search