Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmodd.bas
There are 2 other files named cmodd.bas in the archive. Click here to see a list.
00025 REM CMODD
00060 FILES RFILE1,RFILE2,RFILE3
00080RESTORE#1
00081 INPUT# 1,I1,I2,I3
00100REM
00105IF I3=8 THEN 1390
00110REM
00115REM ****************************************************************
00120L9=0
00125RESTORE#3
00126 INPUT# 3,G,S0
00130IF G=0 THEN 140
00135L9=1
00140REM ******************************************************************
00145REM
00150REM S6=1 IF WHEN ROUTINE IS EXITED YOU WANT TO ENTER
00155REM DIFFERENT INVERSE CHI DISTRIBUTION
00160REM
00165S6=0
00180SCRATCH#1
00181 PRINT # 1,13,I2,I3
00185REM
00190REM INVERSE-CHI DISTRIBUTION EVALUATION
00195REM
00200PRINT L$
00205PRINT " INVERSE CHI DISTRIBUTION EVALUATION"
00210PRINT
00215PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE CHARACTERISTICS OF"
00220PRINT "AN INVERSE CHI DISTRIBUTION."
00225GOTO 330
00235PRINT
00240PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00245PRINT " 1. PERCENTILES"
00250PRINT " 2. HIGHEST DENSITY REGIONS"
00255PRINT " 3. PROBABILITY ABOVE AND BELOW SOME VALUE"
00260PRINT " 4. PROBABILITY BETWEEN TWO VALUES"
00265PRINT " 5. GRAPH OF THE DENSITY FUNCTION"
00270PRINT " 6. EXIT MODULE"
00275GOSUB 9000
00280IF O1=6 THEN 1315
00285I=O1
00290GOTO 725
00320S6=0
00325L9=0
00330IF S6=1 THEN 235
00335IF L9=1 THEN 235
00340S6=1
00365PRINT
00370PRINT "INPUT THE DEGREES OF FREEDOM.(MIN=6 AND MAX=2000)";
00375GOSUB 9000
00380G=O1
00385IF G<6 THEN 440
00390IF G>2000 THEN 440
00395PRINT
00400PRINT "INPUT THE SCALE PARAMETER(SQUARE ROOT OF (MEAN "
00405PRINT "MULTIPLED TIMES THE DF-2)).";
00410GOSUB 9000
00415IF O1 <= 0 THEN 430
00420S0=O1
00425GOTO 235
00430PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00435GOTO 400
00440PRINT "REENTER. MINIMUM IS 6 AND MAXIMUM IS 2000."
00445GOTO 370
00450REM
00455REM END OF ROUTINE TO INPUT THE DISTRIBUTION PARAMETERS
00460REM
00465REM**************************************************************
00470REM
00475REM BEGINNINING OF ROUTINE THAT DETERMINES WHAT YOU WANT TO DO
00480REM AFTER YOU HAVE FINISHED A ROUTINE
00485REM
00490PRINT L$
00495PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00500PRINT " 1. FURTHER EVALUATE THIS DISTRIBUTION."
00505IF L9=1 THEN 550
00510PRINT " 2. EVALUATE ANOTHER INVERSE CHI DISTRIBUTION"
00515PRINT " 3. EXIT MODULE"
00520GOSUB 9000
00525IF O1=3 THEN 1315
00530IF O1=1 THEN 605
00535IF O1=2 THEN 595
00540PRINT "REENTER. INPUT MUST BE 1,2 OR 3."
00545GOTO 520
00550PRINT " 2. EVALUATE THE POSTERIOR DISTRIBUTION ON THE MEAN"
00555PRINT " 3. EXIT MODULE"
00560GOSUB 9000
00565IF O1=2 THEN 590
00570IF O1=3 THEN 1360
00575IF O1=1 THEN 605
00580PRINT "REENTER. INPUT MUST BE 1,2 OR 3."
00585GOTO 560
00590CHAIN "cmodc"
00595S6=0
00600GOTO 610
00605S6=1
00610PRINT L$
00620GOTO 330
00625REM
00630REM END OF ROUTINE TO DETERMINE WHERE YOU ARE GOING
00635REM
00640REM****************************************************************
00645PRINT L$
00650PRINT"----------------------------------------------------------"
00655PRINT " INVERSE-CHI DISTRIBUTION"
00660PRINT USING 665,G,S0
00665:DEGREES OF FREEDOM =#####.## SCALE PARAMETER =#######.###
00670M5=(S0*S0/(G-1.5)) ** .5
00675S5=S0*S0/(2*(G-2)*(G-5/3))
00680: MEAN =#####.## STAN. DEV. =#######.###
00685PRINT USING 680,M5,S5
00690 PRINT"----------------------------------------------------------"
00695RETURN
00700REM ***************************************************************
00705REM
00710REM ROUTINE TO SELECT OPTION
00715REM
00720REM
00725IF I=3 THEN 910
00730IF I=1 THEN 1025
00735IF I=4 THEN 1140
00740IF I=2 THEN 785
00745IF I=5 THEN 1650
00750PRINT
00755PRINT "REENTER. YOU DID NOT SPECIFY 1,2,3,4,5, OR 6."
00760GOTO 275
00765REM
00770REM END OF ROUTINE TO SELECT OPTION
00775REM
00780REM **********************************************************
00785PRINT L$
00790REM
00795PRINT " OPTION 2: HIGHEST DENSITY REGIONS"
00800PRINT
00805PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
00810PRINT "INPUT P% AS NUMBER FROM 1 THROUGH 99."
00815GOSUB 650
00820PRINT "INPUT P%";
00825GOSUB 9000
00830IF O1=0 THEN 490
00835J5=O1/100
00840IF O1<1 THEN 880
00845IF O1>99 THEN 880
00850GOSUB 7300
00855J1=J1*S0
00860J2=J2*S0
00865PRINT USING 870,J5*100,J1,J2
00870: ###.#% HDR =(######.##,######.## )
00875GOTO 820
00880PRINT
00885PRINT "REENTER. INPUT MUST BE 0 OR P% FROM 1 THRU 99."
00890PRINT
00895GOTO 820
00900REM
00905REM ***********************************************************
00910PRINT L$
00915PRINT " OPTION 3: PROBABILITY ABOVE AND BELOW SOME VALUE"
00920PRINT
00925PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
00930PRINT "INPUTTED VALUE MUST BE POSITIVE."
00935GOSUB 650
00940PRINT "INPUT VALUE";
00945GOSUB 9000
00950IF O1=0 THEN 490
00955X3=O1
00960IF X3>0 THEN 975
00965PRINT "REENTER. VALUE MUST BE GREATER THAN 0."
00970GOTO 945
00975X=S0*S0/X3/X3
00980GOSUB 5500
00985X=O1
00990: PROB( X > ######.### ) =##.###
00995PRINT USING 990,O1,P
01000: PROB( X < ######.### ) =##.###
01005PRINT USING 1000,O1,1-P
01010GOTO 940
01015REM
01020REM ********************************************************
01025PRINT L$
01030PRINT " OPTION 1: PERCENTILES"
01035PRINT
01040PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
01045PRINT "INPUT PERCENTILE AS NUMBER FROM .5 THROUGH 99.5."
01050GOSUB 650
01055PRINT "INPUT PERCENTILE";
01060GOSUB 9000
01065IF O1=0 THEN 490
01070P9=O1/100
01075IF O1<.5 THEN 1085
01080IF O1 <= 99.5 THEN 1105
01085PRINT
01090PRINT "REENTER. INPUT MUST BE 0 OR PERCENTILE FROM .5 THRU 99.5."
01095PRINT
01100GOTO 1055
01105P0=1-P9
01110GOSUB 5000
01115: ##.## PERCENTILE =######.##
01120PRINT USING 1115,O1,SQR(1/X)*S0
01125GOTO 1055
01130REM
01135REM ******************************************************************
01140PRINT L$
01145PRINT " OPTION 4: PROBABILITY X IS BETWEEN TWO VALUES"
01150PRINT
01155PRINT "TO EXIT ROUTINE TYPE '0' WHEN YOU ARE ASKED TO INPUT VALUES."
01160GOSUB 650
01165PRINT "INPUT SMALLER VALUE";
01170GOSUB 9000
01175X3=O1
01180IF X3=0 THEN 1210
01185IF X3>0 THEN 1210
01190PRINT
01195PRINT "THE POINTS MUST BE GREATER THAN 0. PLEASE RESPECIFY."
01200PRINT
01205GOTO 1165
01210PRINT "INPUT LARGER VALUE";
01215GOSUB 9000
01220IF O1 <> 0 THEN 1230
01225IF X3=0 THEN 490
01230IF X3=0 THEN 1195
01235X4=O1
01240IF X3<X4 THEN 1260
01245PRINT
01250PRINT "REENTER. INPUT SMALLER VALUE FIRST."
01255GOTO 1165
01260X=S0*S0/X3/X3
01265GOSUB 5500
01270P4=P
01275X=S0*S0/X4/X4
01280GOSUB 5500
01285P3=P
01290PRINT USING 1295,X3,X4,P4-P3
01295: PROB(######.### < X <######.### ) =##.###
01300GOTO 1165
01315IF L9=0 THEN 1360
01320PRINT "IF YOU WANT TO EVALUATE POSTERIOR ON THE MEAN TYPE '1'."
01325PRINT "IF YOU DO NOT TYPE '0'."
01330GOSUB 9000
01335IF O1=0 THEN 1360
01340IF O1=1 THEN 1355
01345PRINT "REENTER. INPUT MUST BE 0 OR 1."
01350GOTO 1330
01355CHAIN "CMODC"
01360CHAIN "RSTRT"
01365REM **************************************************************
01370REM
01375REM EXPECTED UTILITY COMPUTATION
01380REM
01385REM ***************************************************************
01390 DIM M(9),U(9),P(3)
01391 RESTORE #3
01392 FOR I=1 TO 9
01393 INPUT #3,M(I)
01394 NEXT I
01395 FOR I=1 TO 9
01396 INPUT #3,U(I)
01397 NEXT I
01398 FOR I=1 TO 3
01399 INPUT#3,P(I)
01400 NEXT I
01401 PRINTL$
01405S0=P(2)
01410G=P(1)
01415K8=0
01420P0=0
01425K5=2
01430IF M(1)=0 THEN 1460
01435X3=M(1)
01440X=S0*S0/X3/X3
01445GOSUB 5500
01450K8=1-P
01455P0=1-P
01460U0=0
01465U9=0
01470PRINT "INVERSE CHI DISTRIBUTION"
01475PRINT
01480PRINT "DEGREES OF FREEDOM =";G;" SCALE PARAMETER =";S0
01485PRINT
01490PRINT "THE EXPECTED UTILITY IS BEING COMPUTED."
01495FOR I7=M(1) TO M(9) STEP (M(9)-M(1))/110
01500I=I7
01505IF I=0 THEN 1560
01510IF I>M(K5) THEN 1640
01515X=S0*S0/I/I
01520GOSUB 5500
01525P=1-P
01530U1=(I-M(K5-1))*((U(K5)-U(K5-1))/(M(K5)-M(K5-1)))
01535U1=U(K5-1)+U1
01540U2=(U0+U1)/2
01545U9=U9+U2*(P-P0)
01550U0=U1
01555P0=P
01560NEXT I7
01565IF 1-P<1 THEN 1580
01570U9=1
01575GOTO 1600
01580K8=K8+1-P
01585U9=U9/(1-K8)
01590PRINT
01595:THE EXPECTED UTILITY IS ##.##.
01600PRINT USING 1595,U9
01605PRINT
01610PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01615GOSUB 9000
01620 DIM V(1)
01625 SCRATCH #3
01626 PRINT #3,1
01627 FOR I=1 TO 9
01628 PRINT #3,M(I)
01629 NEXT I
01630 FOR I=1 TO 9
01631 PRINT #3,U(I)
01632 NEXT I
01633FOR I=1 TO 3
01634 PRINT #3,P(I)
01635 NEXT I
01636 CHAIN"CMODP"
01640K5=K5+1
01645GOTO 1515
01650PRINT L$
01655PRINT "OPTION 5: GRAPH OF THE DENSITY FUNCTION OVER 99% HDR"
01660GOSUB 650
01665J5=.99
01670PRINT "THESE ARE THE PARAMETERS OF THE DISTRIBUTION TO BE GRAPHED."
01675GOSUB 7307
01680K0=J1*S0
01685K1=J2*S0
01690M0=(S0*S0/(G+1)) ** .5
01695PRINT "WHEN YOU ARE READY FOR THE GRAPH TO BE DISPLAYED TYPE '1'.";
01700GOSUB 9000
01705PRINT L$
01710GOSUB 9400
01715GOTO 490
05000REM***********************************************************
05005REM CHI-SQUARE PERCENTILE
05010REM INPUTS P2
05015REM OUTPUTS J2
05020J=P0
05025P9=P0
05030P2=P0
05035GOSUB 5165
05040X=J2
05045GOSUB 5500
05050IF ABS(P-P9)<.0001 THEN 5160
05055P4=P
05060E7=J2
05065P2=P2-.009
05070GOSUB 5165
05075E1=J2
05080E2=E7
05085X=E1
05090GOSUB 5500
05095P3=P
05100X3=(P9-P3)/(P4-P3)
05105X3=X3*(E2-E1)+E1
05110X=X3
05115GOSUB 5500
05120IF ABS(P9-P)<.0001 THEN 5160
05125IF P>P9 THEN 5145
05130E1=X3
05135P3=P
05140GOTO 5100
05145E2=X3
05150P4=P
05155GOTO 5100
05160RETURN
05165P3=P2
05170IF P2 <= .5 THEN 5180
05175P2=1-P2
05180A1=SQR(LOG(1/P2/P2))
05185A2=2.51552+.802853*A1+.010328*A1*A1
05190A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
05195J2=A1-A2
05200IF P3 <= .5 THEN 5215
05205P2=P3
05210GOTO 5220
05215J2=-J2
05220J2=G*(1-2/9/G+J2*SQR(2/9/G)) ** 3
05225RETURN
05230REM
05235REM********************* END OF CSQ PERCENTILE ********************
05335REM ********************************************************
05340REM SUBROUTINE FOR COMPUTING
05342REM DENSITY OF INVERSE CHI
05344REM ********************************************************
05345G9=.5*G
05350GOSUB 5860
05355F0=LOG(G0)
05360F0=LOG(2)-F0+.5*G*LOG(.5*S0*S0)-(G+1)*LOG(J2)-.5*S0*S0/J2 ** 2
05365D2=F0
05370RETURN
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
05805GOSUB 8000
05810RETURN
05812REM
05813REM END OF CHI-SQUARE CDF ROUTINE
05815REM *********************************************************
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 ****************************************************
07300REM *******************************************************
07301REM INVERSE CHI HIGHEST DENSITY REGION ROUTINE
07302REM INPUTS G J5
07303REM OUTPUTS J1 J2
07304REM
07307J8=1
07308GOSUB 7337
07309X=1/(J1*J1)
07310GOSUB 5500
07311J3=P
07312X=1/(J2*J2)
07313GOSUB 5500
07314J3=J3-P
07315IF ABS(J3-J5)>.0001 THEN 7317
07316RETURN
07317IF J3>J5 THEN 7320
07318J8=J8+1
07319GOTO 7308
07320J9=J8-1
07321J0=J8
07322J8=(J0+J9)/2
07323GOSUB 7337
07324X=1/(J1*J1)
07325GOSUB 5500
07326J3=P
07327X=1/(J2*J2)
07328GOSUB 5500
07329J3=J3-P
07330IF ABS(J3-J5)>.0001 THEN 7332
07331RETURN
07332IF J3>J5 THEN 7335
07333J9=J8
07334GOTO 7322
07335J0=J8
07336GOTO 7322
07337J=J8*(1+EXP(-2*J8/(G+1)))/(1-EXP(-2*J8/(G+1)))
07338J1=1/SQR(J+J8)
07339J2=1/SQR(J-J8)
07340RETURN
07341REM
07342REM END OF ROUTINE FOR INVERSE CHI HDR'S
07343REM**************************************************
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)
08021 IF X*X/2<80 THEN 8025
08022 D=0
08023 GOTO 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"
09419REM**************************************************
09420REM INPUT K0,K1, AND MO
09421REM DENSITY CALL 5345
09422REM******************************************************
09423: DEGREES OF FREEDOM =#######.## SCALE =######.##
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
09522J2=M0
09523PRINT USING 9423,G,S0
09524GOSUB 5345
09526D6=D2
09528IF K1-K2-K0<.01 THEN 9582
09530FOR J2=K0 TO K1-K2 STEP K2
09532GOSUB 5345
09533K7=EXP(D2-D6+LOG(50))
09534 IF K7>=1 THEN 9536
09535 K7=1
09536GOSUB 9458
09538NEXT J2
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"CAN NOT CALCULATE HDR. TYPE '1' TO CONTINUE"
09587GOTO 9555
09593PRINT USING 9423,J5
09999 END