Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-06 - 43,50422/cmod57.bas
There are 2 other files named cmod57.bas in the archive. Click here to see a list.
00011REM********************************************************************
00012REM
00013REM     CMOD57     CMOD57     CMOD57     CMOD57     CMOD57
00014REM
00020REM ****************************************************************
00030REM    GAMMMA     GAMMMA     GAMMMA     GAMMA     GAMMA     GAMMA
00040REM ****************************************************************
00050  FILES RFILE1,RFILE2,RFILE3
00090RESTORE#1
00091  INPUT#  1,I1,I2,I3
00100SCRATCH#1
00101  PRINT #  1,57,I2,I3
00110S6=0
00120PRINT L$
00130PRINT "         EVALUATION OF A GAMMA DISTRIBUTION"
00140PRINT
00150PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE CHARACTERISTICS OF"
00160PRINT "A GAMMA DISTRIBUTION."
00170GOTO 300
00180PRINT L$
00190PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00200PRINT "     1. PERCENTILES"
00210PRINT "     2. HIGHEST DENSITY REGIONS"
00220PRINT "     3. PROBABILITY LESS THAN SOME VALUE"
00230PRINT "     4. PROBABILITY BETWEEN TWO VALUES"
00240PRINT "     5. GRAPH OF DENSITY FUNCTION"
00245PRINT "     6. EXIT MODULE"
00250GOSUB 9000
00260IF O1=6 THEN 1880
00270I=O1
00280PRINT
00290GOTO 710
00300IF S6=1 THEN 180
00310S6=1
00320PRINT
00330PRINT "INPUT DEGREES OF FREEDOM. (MIN=6 AND MAX=2000).";
00340GOSUB 9000
00341A=O1
00350G=O1
00360IF G<6 THEN 460
00370IF G>2000 THEN 460
00380PRINT
00390PRINT "INPUT THE SCALE PARAMETER.(DF/MEAN)";
00400GOSUB 9000
00410IF O1 <= 0 THEN 440
00420S1=O1
00421G=2*G
00422S0=1/(SQR(2*S1))
00430GOTO 180
00440PRINT "REENTER.  SCALE PARAMETER MUST BE POSITIVE."
00450GOTO 390
00460PRINT "REENTER.  MINIMUM IS 6 AND MAXIMUM IS 2000."
00470GOTO 330
00480PRINT L$
00490PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00500PRINT "     1. FURTHER EVALUATE THIS DISTRIBUTION."
00510PRINT "     2. EVALUATE ANOTHER GAMMA DISTRIBUTION."
00520PRINT "     3. EXIT MODULE"
00530GOSUB 9000
00540IF O1=3 THEN 1880
00550IF O1=1 THEN 610
00560IF O1=2 THEN 590
00570PRINT "REENTER.  INPUT MUST BE 1,2 OR 3."
00580GOTO 530
00590S6=0
00600GOTO 620
00610S6=1
00620PRINT L$
00630REM
00640GOTO 300
00650PRINT"----------------------------------------------------------"
00660PRINT "                 GAMMA DISTRIBUTION"
00665M5=G*S0
00667S5=(2*S0*S0*G) ** .5
00670PRINT  USING 680,G/2,S1
00680:DEGREES OF FREEDOM =#####.##      SCALE PARAMETER =#######.###
00685:              MEAN =#####.##           STAN. DEV. =#######.###
00687PRINT  USING 685,(G)*S0*S0,(2*(G)*S0*S0)**.5
00690PRINT"----------------------------------------------------------" 
00700RETURN
00710REM
00720REM  SELECT OPTION
00730IF I=2 THEN 800
00740IF I=1 THEN 1230
00750IF I=4 THEN 1570
00760IF I=3 THEN 1030
00764IF I=5 THEN 2400
00770PRINT
00780PRINT "REENTER.  YOU DID NOT SPECIFY 1,2,3,4,5, OR 6."
00790GOTO 250
00800PRINT L$
00810REM
00820PRINT "            OPTION 2: P% HIGHEST DENSITY REGIONS"
00830PRINT
00840PRINT "INPUT P% AS NUMBER FROM 5 TO 99."
00850PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR P%."
00860GOSUB 650
00870PRINT "INPUT P%";
00880GOSUB 9000
00890IF O1=0 THEN 480
00900J5=O1/100
00910IF O1<5 THEN 990
00920IF O1>99 THEN 990
00930GOSUB 6600
00940J1=J1*S0*S0
00950J2=J2*S0*S0
00960PRINT  USING 970,J5*100,J1,J2
00970:                           ###.#% HDR =(######.##,######.## )
00980GOTO 870
00990PRINT
01000PRINT "REENTER.  INPUT MUST BE AT LEAST 5 AND NOT GREATER THAN 99."
01010PRINT
01020GOTO 870
01030PRINT L$
01040PRINT "         OPTION 3: PROBABILITY LESS THAN SOME VALUE"
01050PRINT
01060PRINT "INPUTTED VALUE MUST BE POSITIVE."
01070PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR VALUE."
01080GOSUB 650
01090PRINT "INPUT VALUE";
01100GOSUB 9000
01110IF O1=0 THEN 480
01120X3=O1
01130IF X3>0 THEN 1160
01140PRINT "THE POINT MUST BE GREATER THAN 0. TRY AGAIN."
01150GOTO 1100
01160X=X3/S0/S0
01170GOSUB 5500
01180:                              PROB( X < ######.## ) =##.##
01190PRINT  USING 1180,O1,P
01200:                              PROB( X > ######.## ) =##.##
01210PRINT  USING 1200,O1,1-P
01220GOTO 1090
01230PRINT L$
01240PRINT "                 OPTION 1: PERCENTILES"
01250PRINT
01260PRINT "INPUT PERCENTILE AS NUMBER FROM .5 THRU 99.5."
01270PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR PERCENTILE."
01280GOSUB 650
01290PRINT "INPUT PERCENTILE";
01300GOSUB 9000
01310IF O1=0 THEN 480
01320P9=O1/100
01330P9=1-P9
01340IF O1<.5 THEN 1360
01350IF O1 <= 99.5 THEN 1400
01360PRINT
01370PRINT "REENTER.  INPUT MUST BE 0 OR PERCENTILE FROM .5 THRU 99.5."
01380PRINT
01390GOTO 1290
01400E2=S0*2
01410E1=0
01420X3=(E1+E2)/2
01430X=S0*S0/X3/X3
01440GOSUB 5500
01450P=1-P
01460IF E2-E1<.0001*S0 THEN 1520
01470IF P>P9 THEN 1500
01480E1=X3
01490GOTO 1420
01500E2=X3
01510GOTO 1420
01520X3=S0**4/X3/X3
01540:                              ##.## PERCENTILE =#########.##
01550PRINT  USING 1540,O1,X3
01560GOTO 1290
01570PRINT L$
01580PRINT "         OPTION 4: PROBABILITY BETWEEN TWO VALUES"
01590PRINT
01600PRINT "BOTH VALUES MUST BE GREATER THAN 0."
01610PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR BOTH VALUES."
01620GOSUB 650
01630PRINT "INPUT SMALLER VALUE";
01640GOSUB 9000
01650X3=O1
01660IF X3=0 THEN 1720
01670IF X3>0 THEN 1720
01680PRINT
01690PRINT "REENTER. THE POINTS MUST BE GREATER THAN 0."
01700PRINT
01710GOTO 1630
01720PRINT "INPUT LARGER VALUE";
01730GOSUB 9000
01740IF O1 <> 0 THEN 1760
01750IF X3=0 THEN 480
01760IF X3=0 THEN 1690
01770X4=O1
01780IF X3 >= X4 THEN 1890
01790X=X3/S0/S0
01800GOSUB 5500
01810P3=P
01820X=X4/S0/S0
01830GOSUB 5500
01840P4=P
01850PRINT  USING 1860,X3,X4,P4-P3
01860:                         PROB(#####.##  TO #####.## ) =##.##
01870GOTO 1630
01880CHAIN "RSTRT"
01890PRINT
01900PRINT "REENTER.  INPUT SMALLER VALUE FIRST."
01910PRINT
01920GOTO 1630
02400PRINT L$
02405PRINT "OPTION 5:  GRAPH OF THE DENSITY FUNCTION OVER 99% HDR"
02410GOSUB 650
02411J5=.99
02412PRINT "THESE ARE THE PARAMETERS OF THE DISTRIBUTION TO BE GRAPHED."
02413GOSUB 6620
02414K0=J1*S0*S0
02415K1=J2*S0*S0
02416K2=(K1-K0)/19
02417M0=S0*(G-2)*S0
02418GOSUB 9400
02420GOTO 490
05345REM   ********************************************************
05350REM       CHI-SQUARE DENSITY FUNCTION
05355REM   ********************************************************
05360G9=.5*G
05365GOSUB 5860
05370F0=G0
05375F0=.5*(G-2)*LOG(J2)-.5*J2/S0/S0
05380D2=F0
05385RETURN
05500REM ********************************************************
05510REM          CHI-SQUARE CDF ROUTINE-LOWER TAIL
05520REM                INPUT     G        X
05530REM                OUTPUT    P
05540REM*************************
05550REM   PEIZER-PRATT-APPROXIMATION
05560REM   PROB(CHISQ<X) WITH DF=G
05570REM*************************
05580X2=X/2
05590T=X2-G*.5+1/3-.04/G
05600Y=(G*.5-.5)/X2
05610IF Y<.00001 THEN 5680
05620IF ABS(Y-1)<.00001 THEN 5700
05630G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05640IF G6>-1 THEN 5710
05650Y3=T
05660GOTO 5720
05670GOTO 5710
05680G6=1
05690GOTO 5710
05700G6=0
05710Y3=T*SQR((1+G6)/X2)
05720GOSUB 8000
05730P2=P
05740RETURN
05850REM ****************************************************
05852REM        LOG GAMMA ROUTINE
05853REM           INPUT G9
05854REM           OUTPUT G0
05860G5=G9
05863IF G9 <= 1.E+30 THEN 5872
05866G0=1.E+38
05869RETURN
05872IF G9>1.E-09 THEN 5881
05875G0=0
05878RETURN
05881IF G9<1.E+10 THEN 5890
05884G0=G9*(LOG(G9)-1)
05887RETURN
05890G6=1
05893IF 18<G5 THEN 5905
05896G6=G6*G5
05899G5=G5+1
05902GOTO 5893
05905R8=1/G5**2
05908G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
05911C1=8.33333E-02
05914C2=2.77778E-03
05917C3=7.93651E-04
05920C4=5.95238E-04
05923G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
05926RETURN
05927REM          END OF LOG GAMMA ROUTINE
05928REM ****************************************************
06600REM     HDR ROUTINE FOR CHI-SQUARE
06605REM            INPUTS G,J5
06610REM              OUTPUTS J1,J2
06615REM*************************************************************
06620J8=1
06625GOSUB 6770
06630X=J1
06635GOSUB 5500
06640J3=P
06645X=J2
06650GOSUB 5500
06655J3=P-J3
06660IF ABS(J3-J5)>.00001 THEN 6670
06665RETURN
06670IF J3>J5 THEN 6685
06675J8=J8+5
06680GOTO 6625
06685J9=J8-5
06690J0=J8
06695J8=(J0+J9)/2
06700GOSUB 6770
06705X=J1
06710GOSUB 5500
06715J3=P
06720X=J2
06725GOSUB 5500
06730J3=P-J3
06732IF ABS(J0-J9)<.000001 THEN 6740
06735IF ABS(J3-J5)>.0001 THEN 6745
06740RETURN
06745IF J3>J5 THEN 6760
06750J9=J8
06755GOTO 6695
06760J0=J8
06765GOTO 6695
06770J=J8*(EXP(2*J8/(G-2))+1)/(EXP(2*J8/(G-2))-1)
06775J1=J-J8
06780J2=J+J8
06785RETURN
06790REM
06795REM       END   OF    CHID-SQUARE    HDR    ROUTINE
06800REM************************************************************
08000REM **********************************************************
08001REM      ROUTINE CALCULATES THE CDF FOR NORMAL DISTRIBUTION
08002REM               INPUT       Y3
08003REM               OUTPUT      P
08004REM
08005Y4=ABS(Y3)
08010X1=X
08015X=Y3
08020T=1/(1+.231642*Y4)
08021IF X*X/2<80 THEN 8025
08022D=0
08023GOTO 8030
08025D=.398942*EXP(-X*X/2)
08030C1=1.33027
08035C2=1.82126
08040C3=1.78148
08045C4=.356564
08050C5=.319382
08055P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
08060IF X >= 0 THEN 8070
08065P=1-P
08070X=X1
08075RETURN
08076REM
08077REM        END OF NORMAL CDF ROUTINE
08078REM **********************************************************
09000REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09035REM*************END ROUTINE
09400 REM
09412  T$=">>>>>>>>>1>>>>>>>>>2>>>>>>>>>3>>>>>>>>>4>>>>>>>>>5>>>>>>"
09414  S$="\\\\\\\\\I\\\\\\\\\I\\\\\\\\\I\\\\\\\\\I\\\\\\\\\I"
09416  U$="/////////I/////////I/////////I/////////I/////////I/////////I"
09419REM**************************************************
09420REM INPUT K0,K1, AND MO
09421REM DENSITY CALL 5345
09422REM******************************************************
09426GOTO 9520
09450  :#####.## I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLCONTINUE=1
09458IF ABS(J2)>9999 THEN 9560
09459IF ABS(J2)<.01 THEN 9560
09460IF J2>M0+.4*K2 THEN 9472
09461IF J2<M0-.4*K2 THEN 9466
09462  PRINT  USING 9466,J2,MID$(T$,1,K7-(1)+1)
09464GOTO 9474
09466:#####.## I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
09468  PRINT  USING 9466,J2,MID$(S$,1,K7-(1)+1)
09470GOTO 9474
09472  PRINT  USING 9466,J2,MID$(U$,1,K7-(1)+1)
09474RETURN
09510REM ***************************************************************
09520K2=(K1-K0)/19
09521PRINT "WHEN YOU ARE READY FOR THE GRAPH TO BE DISPLAYED TYPE '1'.";
09522J2=M0
09523GOSUB 9000
09524GOSUB 5345
09525PRINT L$
09526D6=D2
09527PRINT  USING 9539,J5*100
09528IF K1-K2-K0<.01 THEN 9582
09530FOR J2=K0 TO K1-.9*K2 STEP K2
09532GOSUB 5345
09533K7=EXP(D2-D6+LOG(50)) 
09534 IF K7>=1 THEN 9536
09535 K7=1
09536GOSUB 9458
09538NEXT J2
09539:         GRAPH OF GAMMA ##.# % HDR
09540J2=K1
09542GOSUB 5345
09543K7=EXP(D2-D6+LOG(50)) 
09544 IF K7>=1 THEN 9546
09545 K7=1
09546IF J2<.01 THEN 9575
09547IF J2>9999.99 THEN 9575
09548  PRINT  USING 9450,K1,MID$(U$,1,K7-(1)+1)
09555GOSUB 9000
09556RETURN
09560IF J2>M0+.4*K2 THEN 9572
09561IF J2<M0-.4*K2 THEN 9566
09562  PRINT  USING 9566,J2,MID$(T$,1,K7-(1)+1)
09564RETURN
09566:##.##^^^^ I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
09568  PRINT  USING 9566,J2,MID$(S$,1,K7-(1)+1)
09570RETURN
09572  PRINT  USING 9566,J2,MID$(U$,1,K7-(1)+1)
09574RETURN
09575  PRINT  USING 9580,K1,MID$(U$,1,K7-(1)+1)
09577GOTO 9555
09580  :##.##^^^^ I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLCONTINUE=1
09582PRINT "THE HIGHEST DENSITY REGION IS WITHIN TO SMALL A RANGE TO BE"
09583PRINT "ACCURATELY DISPLAYED.  THE UPPER AND LOWER BOUNDS ARE"
09584PRINT "WITHIN .00X AND .00Y OF EACH OTHER.  THIS INTERVAL CAN BE"
09585PRINT "OBTAINED FROM THE HIGHEST DENSITY REGION MODULE."
09586PRINT "WHEN YOU WANT TO CONTINUE TYPE '1'";
09587GOTO 9555
09999 END