Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmodn.bas
There are 2 other files named cmodn.bas in the archive. Click here to see a list.
00025REM ***************************************************************
00030REM
00035REM INVERSE CHI-SQUARE DISTRIBUTION
00040REM
00045REM ***************************************************************
00050REM CMODN CMODN CMODN CMODN CMODN CMODN
00055REM*************************************************************
00060 FILES RFILE1,RFILE2,RFILE3
00080RESTORE#1
00081 INPUT# 1,I1,I2,I3
00085RESTORE#3
00086 INPUT# 3,G,S0
00090SCRATCH#1
00091 PRINT # 1,23,I2,I3
00095S6=0
00100REM ************************************************************
00105PRINT L$
00110PRINT " EVALUATION OF INVERSE CHI-SQUARE DISTRIBUTION"
00115PRINT
00120PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE CHARACTERISTICS OF"
00125PRINT "AN INVERSE CHI-SQUARE DISTRIBUTION."
00130GOTO 215
00135PRINT L$
00140PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00145PRINT " 1. PERCENTILES"
00150PRINT " 2. HIGHEST DENSITY REGIONS"
00155PRINT " 3. PROBABILITY ABOVE AND BELOW SOME VALUE"
00160PRINT " 4. PROBABILITY BETWEEN TWO VALUES"
00165PRINT " 5. GRAPH OF THE DENSITY FUCNCTION"
00170PRINT " 6. EXIT MODULE"
00175GOSUB 9000
00180IF O1=6 THEN 1140
00185I=O1
00190GOTO 560
00195REM
00200REM S6=1 IF THE PARAMETERS OF THE DISTRIBUTION ARE ALREADY IN
00205REM S6=0 IF THE PARAMETERS ARE TO BE ENTERED
00210REM
00215IF S6=1 THEN 135
00220S6=1
00225REM *****************************************************************
00230REM
00235REM BEGINNING OF ROUTINE TO INPUT DISTRIBUTION PARAMETERS
00240REM
00245PRINT
00250PRINT "INPUT THE DEGREES OF FREEDOM.(MIN=6 AND MAX=2000)";
00255GOSUB 9000
00260G=O1
00265IF G>2000 THEN 325
00270IF G<6 THEN 315
00275PRINT
00280PRINT "INPUT THE SCALE PARAMETER(MEAN MULTIPLIED BY (DF-2)).";
00285GOSUB 9000
00290IF O1 <= 0 THEN 305
00295S0=SQR(O1)
00300GOTO 135
00305PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00310GOTO 275
00315PRINT "REENTER. DEGREES OF FREEDOM MUST BE AT LEAST 6."
00320GOTO 245
00325PRINT "REENTER. DEGREES OF FREEDOM CAN NOT BE GREATER THAN 2000."
00330GOTO 245
00335REM
00340REM END OF ROUTINE TO INPUT THE DISTRIBUTION PARAMETERS
00345REM
00350REM**************************************************************
00355REM
00360REM BEGINNINING OF ROUTINE THAT DETERMINES WHAT YOU WANT TO DO
00365REM AFTER YOU HAVE FINISHED A ROUTINE
00370REM
00375PRINT L$
00380PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00385PRINT " 1. FURTHER EVALUATE THIS DISTRIBUTION."
00390PRINT " 2. EVALUATE ANOTHER INVERSE CHI-SQUARE DISTRIBUTION"
00395PRINT " 3. EXIT MODULE"
00400GOSUB 9000
00405IF O1=3 THEN 1140
00410IF O1=1 THEN 440
00415IF O1=2 THEN 430
00420PRINT "REENTER. INPUT MUST BE 1,2 OR 3."
00425GOTO 400
00430S6=0
00435GOTO 445
00440S6=1
00445PRINT L$
00450 REM
00455GOTO 215
00460REM
00465REM END OF ROUTINE TO DETERMINE WHERE YOU ARE GOING
00470REM
00475REM****************************************************************
00480PRINT L$
00485PRINT"----------------------------------------------------------"
00490M5=S0*S0/(G-2)
00495S5=(2*(S0*S0) ** 2/((G-2) ** 2*(G-4))) ** .5
00500PRINT " INVERSE CHI-SQUARE DISTRIBUTION."
00505PRINT USING 510,G,S0*S0
00510:DEGREES OF FREEDOM =#####.## SCALE PARAMETER =#######.###
00515: MEAN =#####.## STAN. DEV. =#######.###
00520PRINT USING 515,M5,S5
00525 PRINT"-----------------------------------------------------------"
00530RETURN
00535REM ***************************************************************
00540REM
00545REM ROUTINE TO SELECT OPTION
00550REM
00555REM
00560IF I=3 THEN 745
00565IF I=1 THEN 860
00570IF I=4 THEN 975
00575IF I=2 THEN 620
00580IF I=5 THEN 1170
00585PRINT
00590PRINT "REENTER. YOU DID NOT SPECIFY 1,2,3,4 OR 5."
00595GOTO 175
00600REM
00605REM END OF ROUTINE TO SELECT OPTION
00610REM
00615REM **********************************************************
00620PRINT L$
00625REM
00630PRINT " OPTION 2: P% HIGHEST DENSITY REGIONS"
00635PRINT
00640PRINT "INPUT P% AS NUMBER FROM 5 THROUGH 99."
00645PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR P%."
00650GOSUB 485
00655PRINT "INPUT P%";
00660GOSUB 9000
00665IF O1=0 THEN 375
00670J5=O1/100
00675IF O1<5 THEN 715
00680IF O1>99 THEN 715
00685GOSUB 7800
00690J1=J1*S0*S0
00695J2=J2*S0*S0
00700PRINT USING 705,J5*100,J1,J2
00705: ###.#% HDR =(######.##,######.## )
00710GOTO 655
00715PRINT
00720PRINT "REENTER. INPUT MUST BE 0 OR P% FROM 1 THRU 99."
00725PRINT
00730GOTO 655
00735REM
00740REM ***********************************************************
00745PRINT L$
00750PRINT " OPTION 3: PROBABILITY ABOVE AND BELOW SOME VALUE"
00755PRINT
00760PRINT "ONLY VALUES GREATER THAN 0 ARE ACCEPTABLE."
00765PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR VALUE."
00770GOSUB 485
00775PRINT "INPUT VALUE";
00780GOSUB 9000
00785IF O1=0 THEN 375
00790X3=O1
00795IF X3>0 THEN 810
00800PRINT "REENTER. VALUE MUST BE GREATER THAN 0."
00805GOTO 780
00810X=S0*S0/X3
00815GOSUB 5500
00820X=O1
00825: PROB( X > ######.### ) =##.##
00830PRINT USING 825,O1,P
00835: PROB( X < ######.### ) =##.##
00840PRINT USING 835,O1,1-P
00845GOTO 775
00850REM
00855REM ********************************************************
00860PRINT L$
00865PRINT " OPTION 1: PERCENTILES"
00870PRINT
00875PRINT "INPUT PERCENTILE AS NUMBER FROM .5 THROUGH 99.5."
00880PRINT "WHEN YOU WANT TO EXIT ROUTINE INPUT '0' FOR PERCENTILE."
00885GOSUB 485
00890PRINT "INPUT PERCENTILE";
00895GOSUB 9000
00900IF O1=0 THEN 375
00905P9=O1/100
00910IF O1<.5 THEN 920
00915IF O1 <= 99.5 THEN 940
00920PRINT
00925PRINT "REENTER. INPUT MUST BE ZERO OR ACCEPTABLE PERCENTILE."
00930PRINT
00935GOTO 890
00940P0=1-P9
00945GOSUB 5000
00950: ##.## PERCENTILE =######.##
00955PRINT USING 950,O1,1/X*S0*S0
00960GOTO 890
00965REM
00970REM ******************************************************************
00975PRINT L$
00980PRINT " OPTION 4: PROBABILITY BETWEEN TWO VALUES"
00985PRINT
00990PRINT "BOTH VALUES MUST BE GREATER THAN 0."
00995PRINT "WHEN YOU WANT TO EXIT ROUTINE TYPE '0' FOR BOTH VALUES."
01000GOSUB 485
01005PRINT "INPUT SMALLER VALUE";
01010GOSUB 9000
01015X3=O1
01020IF X3=0 THEN 1050
01025IF X3>0 THEN 1050
01030PRINT
01035PRINT "REENTER. THE POINTS MUST BE GREATER THAN 0."
01040PRINT
01045GOTO 1005
01050PRINT "INPUT LARGER VALUE";
01055GOSUB 9000
01060IF O1 <> 0 THEN 1070
01065IF X3=0 THEN 375
01070IF X3=0 THEN 1035
01075X4=O1
01080IF X3 >= X4 THEN 1150
01085X=S0*S0/X3
01090GOSUB 5500
01095P4=P
01100X=S0*S0/X4
01105GOSUB 5500
01110P3=P
01115: PROB(######.### TO ######.###) =##.##
01120PRINT USING 1115,X3,X4,P4-P3
01125GOTO 1005
01130REM
01135REM ***************************************************************
01140CHAIN "RSTRT"
01145REM *************************************************************
01150PRINT
01155PRINT "REENTER. INPUT SMALLER VALUE FIRST."
01160PRINT
01165GOTO 1005
01170PRINT L$
01175PRINT "OPTION 5: GRAPH OF THE DENSITY FUNCTION OVER 99% HDR"
01180GOSUB 485
01185M0=S0*S0/(G+2)
01190J5=.99
01195PRINT "THESE ARE THE PARAMETERS OF THE DISTRIBUTION TO BE GRAPHED."
01200GOSUB 7800
01205K0=J1*S0*S0
01210K1=J2*S0*S0
01215PRINT "WHEN YOU ARE READY FOR THE GRAPH TO BE DISPLAYED TYPE '1'.";
01220GOSUB 9000
01225PRINT L$
01230GOSUB 9400
01235GOTO 375
01240F0=LOG(G0)
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 ********************
05345G9=.5*G
05348GOSUB 5850
05349F0=G0
05350F0=-F0+.5*G*LOG(.5*S0*S0)-.5*(G+2)*LOG(J2)-.5*S0*S0/J2
05355D2=F0
05360RETURN
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
05800Y3=X2
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 ****************************************************
07800REM *******************************************************
07802REM INVERSE CHI-SQUARE HIGHEST DENSITY REGION
07804REM INPUTS G J5
07806REM OUTPUTS J1 J2
07808REM
07814J8=5.1
07816GOSUB 7874
07818X=1/J2
07820GOSUB 5500
07822J3=P
07824X=1/J1
07826GOSUB 5500
07828J3=P-J3
07830IF ABS(J3-J5)>.00001 THEN 7834
07832RETURN
07834IF J3>J5 THEN 7840
07836J8=J8+5
07838GOTO 7816
07840J9=J8-5
07842J0=J8
07844J8=(J0+J9)/2
07846GOSUB 7874
07848X=1/J2
07850GOSUB 5500
07852J3=P
07854X=1/J1
07856GOSUB 5500
07858J3=P-J3
07859IF ABS(J1-J2)<.0001 THEN 7862
07860IF ABS(J3-J5)>.00001 THEN 7864
07862RETURN
07864IF J3>J5 THEN 7870
07866J9=J8
07868GOTO 7844
07870J0=J8
07872GOTO 7844
07874J=J8*(EXP(2*J8/(G+2))+1)/(EXP(2*J8/(G+2))-1)
07876J1=1/(J+J8)
07878J2=1/(J-J8)
07880RETURN
07882REM
07884REM END OF INVERSE CHI-SUARE HDR ROUTINE
07886REM***********************************************************
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: GRAPH OF INVERSE CHI SQUARE ##.# % HDR
09426GOTO 9520
09450 :#####.## I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL CONTINUE=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,J5*100
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'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL CONTINUE=1
09582PRINT "THE HIGHEST DENSITY REGION IS WITHIN TOO 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 ARE READY TO CONTINUE TYPE '1'.";
09587GOTO 9555
09999 END