Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/bdtr.ssp
There are 2 other files named bdtr.ssp in the archive. Click here to see a list.
C                                                                       BDTR  10
C     ..................................................................BDTR  20
C                                                                       BDTR  30
C        SUBROUTINE BDTR                                                BDTR  40
C                                                                       BDTR  50
C        PURPOSE                                                        BDTR  60
C           COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,     BDTR  70
C           DISTRIBUTED ACCORDING TO THE BETA DISTRIBUTION WITH         BDTR  80
C           PARAMETERS A AND B, IS LESS THAN OR EQUAL TO X.  F(A,B,X),  BDTR  90
C           THE ORDINATE OF THE BETA DENSITY AT X, IS ALSO COMPUTED.    BDTR 100
C                                                                       BDTR 110
C        USAGE                                                          BDTR 120
C           CALL BDTR(X,A,B,P,D,IER)                                    BDTR 130
C                                                                       BDTR 140
C        DESCRIPTION OF PARAMETERS                                      BDTR 150
C           X   - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.              BDTR 160
C           A   - BETA DISTRIBUTION PARAMETER (CONTINUOUS).             BDTR 170
C           B   - BETA DISTRIBUTION PARAMETER (CONTINUOUS).             BDTR 180
C           P   - OUTPUT PROBABILITY.                                   BDTR 190
C           D   - OUTPUT DENSITY.                                       BDTR 200
C           IER - RESULTANT ERROR CODE WHERE                            BDTR 210
C               IER= 0 --- NO ERROR                                     BDTR 220
C               IER=-1,+1  CDTR HAS BEEN CALLED AND AN ERROR HAS        BDTR 230
C                          OCCURRED.  SEE CDTR.                         BDTR 240
C               IER=-2 --- AN INPUT PARAMETER IS INVALID.  X IS LESS    BDTR 250
C                          THAN 0.0 OR GREATER THAN 1.0, OR EITHER A OR BDTR 260
C                          B IS LESS THAN 0.5 OR GREATER THAN 10**(+5). BDTR 270
C                          P AND D ARE SET TO -1.7E38.                   BDTR 280
C               IER=+2 --- INVALID OUTPUT.  P IS LESS THAN ZERO OR      BDTR 290
C                          GREATER THAN ONE.  P IS SET TO 1.7E38.        BDTR 300
C                                                                       BDTR 310
C        REMARKS                                                        BDTR 320
C           SEE MATHEMATICAL DESCRIPTION.                               BDTR 330
C                                                                       BDTR 340
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  BDTR 350
C           DLGAM                                                       BDTR 360
C           NDTR                                                        BDTR 370
C           CDTR                                                        BDTR 380
C                                                                       BDTR 390
C        METHOD                                                         BDTR 400
C           REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL          BDTR 410
C           DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,              BDTR 420
C           IBM RESEARCH REPORT RC-1094, 1963.                          BDTR 430
C                                                                       BDTR 440
C     ..................................................................BDTR 450
C                                                                       BDTR 460
      SUBROUTINE BDTR(X,A,B,P,D,IER)                                    BDTR 470
      DOUBLE PRECISION XX,DLXX,DL1X,AA,BB,G1,G2,G3,G4,DD,PP,XO,FF,FN,   BDTR 480
     1XI,SS,CC,RR,DLBETA                                                BDTR 490
C                                                                       BDTR 500
C        TEST FOR VALID INPUT DATA                                      BDTR 510
C                                                                       BDTR 520
      IF(A-(.5-1.E-5)) 640,10,10                                        BDTR 530
   10 IF(B-(.5-1.E-5)) 640,20,20                                        BDTR 540
   20 IF(A-1.E+5) 30,30,640                                             BDTR 550
   30 IF(B-1.E+5) 40,40,640                                             BDTR 560
   40 IF(X) 640,50,50                                                   BDTR 570
   50 IF(1.-X) 640,60,60                                                BDTR 580
C                                                                       BDTR 590
C        COMPUTE LOG(BETA(A,B))                                         BDTR 600
C                                                                       BDTR 610
   60 AA=DBLE(A)                                                        BDTR 620
      BB=DBLE(B)                                                        BDTR 630
      CALL DLGAM(AA,G1,IOK)                                             BDTR 640
      CALL DLGAM(BB,G2,IOK)                                             BDTR 650
      CALL DLGAM(AA+BB,G3,IOK)                                          BDTR 660
      DLBETA=G1+G2-G3                                                   BDTR 670
C                                                                       BDTR 680
C        TEST FOR X NEAR 0.0 OR 1.0                                     BDTR 690
C                                                                       BDTR 700
      IF(X-1.E-8) 80,80,70                                              BDTR 710
   70 IF((1.-X)-1.E-8) 130,130,140                                      BDTR 720
   80 P=0.0                                                             BDTR 730
      IF(A-1.) 90,100,120                                               BDTR 740
   90 D=1.7E38                                                          BDTR 750
      GO TO 660                                                         BDTR 760
  100 DD=-DLBETA                                                        BDTR 770
      IF(DD+1.68D02)  120,120,110                                       BDTR 780
  110 DD=DEXP(DD)                                                       BDTR 790
      D=SNGL(DD)                                                        BDTR 800
      GO TO 660                                                         BDTR 810
  120 D=0.0                                                             BDTR 820
      GO TO 660                                                         BDTR 830
  130 P=1.0                                                             BDTR 840
      IF(B-1.) 90,100,120                                               BDTR 850
C                                                                       BDTR 860
C        SET PROGRAM PARAMETERS                                         BDTR 870
C                                                                       BDTR 880
  140 XX=DBLE(X)                                                        BDTR 890
      DLXX=DLOG(XX)                                                     BDTR 900
      DL1X=DLOG(1.D0-XX)                                                BDTR 910
      XO=XX/(1.D0-XX)                                                   BDTR 920
      ID=0                                                              BDTR 930
C                                                                       BDTR 940
C        COMPUTE ORDINATE                                               BDTR 950
C                                                                       BDTR 960
      DD=(AA-1.D0)*DLXX+(BB-1.D0)*DL1X-DLBETA                           BDTR 970
      IF(DD-1.68D02) 150,150,160                                        BDTR 980
  150 IF(DD+1.68D02) 170,170,180                                        BDTR 990
  160 D=1.7E38                                                           BDTR1000
      GO TO 190                                                         BDTR1010
  170 D=0.0                                                             BDTR1020
      GO TO 190                                                         BDTR1030
  180 DD=DEXP(DD)                                                       BDTR1040
      D=SNGL(DD)                                                        BDTR1050
C                                                                       BDTR1060
C        A OR B OR BOTH WITHIN 1.E-8 OF 1.0                             BDTR1070
C                                                                       BDTR1080
  190 IF(ABS(A-1.)-1.E-8)  200,200,210                                  BDTR1090
  200 IF(ABS(B-1.)-1.E-8)  220,220,230                                  BDTR1100
  210 IF(ABS(B-1.)-1.E-8)  260,260,290                                  BDTR1110
  220 P=X                                                               BDTR1120
      GO TO 660                                                         BDTR1130
  230 PP=BB*DL1X                                                        BDTR1140
      IF(PP+1.68D02) 240,240,250                                        BDTR1150
  240 P=1.0                                                             BDTR1160
      GO TO 660                                                         BDTR1170
  250 PP=DEXP(PP)                                                       BDTR1180
      PP=1.D0-PP                                                        BDTR1190
      P=SNGL(PP)                                                        BDTR1200
      GO TO 600                                                         BDTR1210
  260 PP=AA*DLXX                                                        BDTR1220
      IF(PP+1.68D02) 270,270,280                                        BDTR1230
  270 P=0.0                                                             BDTR1240
      GO TO 660                                                         BDTR1250
  280 PP=DEXP(PP)                                                       BDTR1260
      P=SNGL(PP)                                                        BDTR1270
      GO TO 600                                                         BDTR1280
C                                                                       BDTR1290
C        TEST FOR A OR B GREATER THAN 1000.0                            BDTR1300
C                                                                       BDTR1310
  290 IF(A-1000.) 300,300,310                                           BDTR1320
  300 IF(B-1000.) 330,330,320                                           BDTR1330
  310 XX=2.D0*AA/XO                                                     BDTR1340
      XS=SNGL(XX)                                                       BDTR1350
      AA=2.D0*BB                                                        BDTR1360
      DF=SNGL(AA)                                                       BDTR1370
      CALL CDTR(XS,DF,P,DUMMY,IER)                                      BDTR1380
      P=1.0-P                                                           BDTR1390
      GO TO 670                                                         BDTR1400
  320 XX=2.D0*BB*XO                                                     BDTR1410
      XS=SNGL(XX)                                                       BDTR1420
      AA=2.D0*AA                                                        BDTR1430
      DF=SNGL(AA)                                                       BDTR1440
      CALL CDTR(XS,DF,P,DUMMY,IER)                                      BDTR1450
      GO TO 670                                                         BDTR1460
C                                                                       BDTR1470
C        SELECT PARAMETERS FOR CONTINUED FRACTION COMPUTATION           BDTR1480
C                                                                       BDTR1490
  330 IF(X-.5) 340,340,380                                              BDTR1500
  340 IF(AA-1.D0) 350,350,360                                           BDTR1510
  350 RR=AA+1.D0                                                        BDTR1520
      GO TO 370                                                         BDTR1530
  360 RR=AA                                                             BDTR1540
  370 DD=DLXX/5.D0                                                      BDTR1550
      DD=DEXP(DD)                                                       BDTR1560
      DD=(RR-1.D0)-(RR+BB-1.D0)*XX*DD +2.D0                             BDTR1570
      IF(DD) 420,420,430                                                BDTR1580
  380 IF(BB-1.D0) 390,390,400                                           BDTR1590
  390 RR=BB+1.D0                                                        BDTR1600
      GO TO 410                                                         BDTR1610
  400 RR=BB                                                             BDTR1620
  410 DD=DL1X/5.D0                                                      BDTR1630
      DD=DEXP(DD)                                                       BDTR1640
      DD=(RR-1.D0)-(AA+RR-1.D0)*(1.D0-XX)*DD +2.D0                      BDTR1650
      IF(DD) 430,430,420                                                BDTR1660
  420 ID=1                                                              BDTR1670
      FF=DL1X                                                           BDTR1680
      DL1X=DLXX                                                         BDTR1690
      DLXX=FF                                                           BDTR1700
      XO=1.D0/XO                                                        BDTR1710
      FF=AA                                                             BDTR1720
      AA=BB                                                             BDTR1730
      BB=FF                                                             BDTR1740
      G2=G1                                                             BDTR1750
C                                                                       BDTR1760
C        TEST FOR A LESS THAN 1.0                                       BDTR1770
C                                                                       BDTR1780
  430 FF=0.D0                                                           BDTR1790
      IF(AA-1.D0) 440,440,470                                           BDTR1800
  440 CALL DLGAM(AA+1.D0,G4,IOK)                                        BDTR1810
      DD=AA*DLXX+BB*DL1X+G3-G2-G4                                       BDTR1820
      IF(DD+1.68D02) 460,460,450                                        BDTR1830
  450 FF=FF+DEXP(DD)                                                    BDTR1840
  460 AA=AA+1.D0                                                        BDTR1850
C                                                                       BDTR1860
C        COMPUTE P USING CONTINUED FRACTION EXPANSION                   BDTR1870
C                                                                       BDTR1880
  470 FN=AA+BB-1.D0                                                     BDTR1890
      RR=AA-1.D0                                                        BDTR1900
      II=80                                                             BDTR1910
      XI=DFLOAT(II)                                                     BDTR1920
      SS=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))             BDTR1930
      SS=SS*XO                                                          BDTR1940
      DO 480 I=1,79                                                     BDTR1950
      II=80-I                                                           BDTR1960
      XI=DFLOAT(II)                                                     BDTR1970
      DD=(XI*(FN+XI))/((RR+2.D0*XI+1.D0)*(RR+2.D0*XI))                  BDTR1980
      DD=DD*XO                                                          BDTR1990
      CC=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))             BDTR2000
      CC=CC*XO                                                          BDTR2010
      SS=CC/(1.D0+DD/(1.D0-SS))                                         BDTR2020
  480 CONTINUE                                                          BDTR2030
      SS=1.D0/(1.D0-SS)                                                 BDTR2040
      IF(SS) 650,650,490                                                BDTR2050
  490 CALL DLGAM(AA+BB,G1,IOK)                                          BDTR2060
      CALL DLGAM(AA+1.D0,G4,IOK)                                        BDTR2070
      CC=G1-G2-G4+AA*DLXX+(BB-1.D0)*DL1X                                BDTR2080
      PP=CC+DLOG(SS)                                                    BDTR2090
      IF(PP+1.68D02) 500,500,510                                        BDTR2100
  500 PP=FF                                                             BDTR2110
      GO TO 520                                                         BDTR2120
  510 PP=DEXP(PP)+FF                                                    BDTR2130
  520 IF(ID) 540,540,530                                                BDTR2140
  530 PP=1.D0-PP                                                        BDTR2150
  540 P=SNGL(PP)                                                        BDTR2160
C                                                                       BDTR2170
C        SET ERROR INDICATOR                                            BDTR2180
C                                                                       BDTR2190
      IF(P) 550,570,570                                                 BDTR2200
  550 IF(ABS(P)-1.E-7) 560,560,650                                      BDTR2210
  560 P=0.0                                                             BDTR2220
      GO TO 660                                                         BDTR2230
  570 IF(1.-P) 580,600,600                                              BDTR2240
  580 IF(ABS(1.-P)-1.E-7) 590,590,650                                   BDTR2250
  590 P=1.0                                                             BDTR2260
      GO TO 660                                                         BDTR2270
  600 IF(P-1.E-8) 610,610,620                                           BDTR2280
  610 P=0.0                                                             BDTR2290
      GO TO 660                                                         BDTR2300
  620 IF((1.0-P)-1.E-8) 630,630,660                                     BDTR2310
  630 P=1.0                                                             BDTR2320
      GO TO 660                                                         BDTR2330
  640 IER=-2                                                            BDTR2340
      D=-1.7E38                                                          BDTR2350
      P=-1.7E38                                                          BDTR2360
      GO TO 670                                                         BDTR2370
  650 IER=+2                                                            BDTR2380
      P= 1.7E38                                                          BDTR2390
      GO TO 670                                                         BDTR2400
  660 IER=0                                                             BDTR2410
  670 RETURN                                                            BDTR2420
      END                                                               BDTR2430