Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmodo.bas
There are 2 other files named cmodo.bas in the archive. Click here to see a list.
00020REM ****************************************************************
00030REM CMODO CMODO CMODO CMODO CMODO CMODO
00040REM ****************************************************************
00050 FILES RFILE1,RFILE2,RFILE3
00090RESTORE#1
00091 INPUT# 1,I1,I2,I3
00100SCRATCH#1
00101 PRINT # 1,24,I2,I3
00110S6=0
00120PRINT L$
00130PRINT " EVALUATION OF A CHI-SQUARE DISTRIBUTION"
00140PRINT
00150PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE CHARACTERISTICS OF"
00160PRINT "A CHI-SQUARE 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=4 AND MAX=2000).";
00340GOSUB 9000
00350G=O1
00360IF G<4 THEN 460
00370IF G>2000 THEN 460
00380PRINT
00390PRINT "INPUT SCALE PARAMETER(MEAN DIVIDED BY DF).";
00400GOSUB 9000
00410IF O1 <= 0 THEN 440
00420S0=SQR(O1)
00430GOTO 180
00440PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00450GOTO 390
00460PRINT "REENTER. MINIMUM IS 4 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 CHI-SQUARE 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$
00630 REM
00640GOTO 300
00650PRINT"---------------------------------------------------------"
00660PRINT " CHI-SQUARE DISTRIBUTION"
00665M5=S0*S0*G
00667S5=SQR(2*S0*S0*S0*S0*G)
00670PRINT USING 680,G,S0*S0
00680:DEGREES OF FREEDOM =#####.## SCALE PARAMETER =#######.###
00685: MEAN =#####.## STAN. DEV. =#######.###
00687PRINT USING 685,M5,S5
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
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
01400P0=P9
01401GOSUB 5000
01540: ##.## PERCENTILE =#########.##
01550PRINT USING 1540,O1,X*S0*S0
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 GRAPHED."
02413GOSUB 6620
02414K0=J1*S0*S0
02415K1=J2*S0*S0
02416K2=(K1-K0)/19
02417M0=S0*S0*(G-2)
02418PRINT "WHEN YOU ARE READY FOR THE GRAPH TO BE DISPLAYED TYPE '1'.";
02419GOSUB 9000
02420PRINT L$
02421GOSUB 9400
02422PRINT L$
02423GOTO 490
05000REM***********************************************************
05005REM CHI-SQUARE PERCENTILE
05010REM INPUTS P2
05015REM OUTPUTS J2
05020J=P0
05025P9=P0
05030P2=P0
05035GOSUB 5165
05040X=J2
05045GOSUB 5500
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
05112IF ABS(E1-E2)<.001 THEN 5160
05113IF ABS(P9-P)<.0009 THEN 5160
05115GOSUB 5500
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 ********************
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)>.00001 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"
09419REM**************************************************
09420REM INPUT K0,K1, AND MO
09421REM DENSITY CALL 5345
09422REM******************************************************
09423: DF=######.## SCALE PARAMETER=#######.##
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*S0
09524GOSUB 5345
09526D6=D2
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
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
09690PRINT"---------------------------------------------------------"
09999 END