Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmodu.bas
There are 2 other files named cmodu.bas in the archive. Click here to see a list.
00020REM******************************************************************
00030REM CMODU CMODU CMODU CMODU CMODU CMODU
00040REM******************************************************************
00050 FILES RFILE1,RFILE2,RFILE3
00060 DIM L(10)
00100RESTORE#1
00101 INPUT# 1,I1,I2,I3
00110SCRATCH#1
00111 PRINT # 1,30,I2,I3
00120X=0
00130GOSUB 5160
00140PRINT L$
00150PRINT " COMPARISON OF TWO STANDARD DEVIATIONS"
00160PRINT
00170PRINT "THIS MODULE COMPUTES THE PROBABILITY THAT THE RATIO OF TWO"
00180PRINT "STANDARD DEVIATIONS DISTRIBUTED AS INDEPENDENT INVERSE"
00190PRINT "CHIS IS GREATER THAN K, WHERE K CAN TAKE ON ANY POSITIVE"
00200PRINT "VALUE."
00210PRINT
00220PRINT
00230REM*****************************************************************
00240PRINT "INPUT THE PARAMETERS OF THE INVERSE CHI DISTRIBUTION ON THE"
00250PRINT "FIRST STANDARD DEVIATION."
00260REM*****************************************************************
00270PRINT
00280PRINT "INPUT THE DEGREES OF FREEDOM.";
00290GOSUB 9000
00300IF O1>2 THEN 330
00310PRINT "REENTER. THE DEGREES OF FREEDOM MUST BE GREATER THAN 2."
00320GOTO 290
00330V1=O1
00340PRINT "INPUT THE SCALE PARAMETER.";
00350GOSUB 9000
00360IF O1>0 THEN 390
00370PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00380GOTO 350
00390L1=O1
00400PRINT
00410REM*****************************************************************
00420PRINT "INPUT THE PARAMETERS OF THE SECOND INVERSE CHI DISTRIBUTION."
00430REM******************************************************************
00440PRINT
00450PRINT "INPUT THE DEGREES OF FREEDOM";
00460GOSUB 9000
00470IF O1>2 THEN 500
00480PRINT "REENTER. DEGREES OF FREEEDOM MUST BE GREATER THAN 2."
00490GOTO 460
00500V2=O1
00510PRINT "INPUT THE SCALE PARAMETER.";
00520GOSUB 9000
00530IF O1>0 THEN 560
00540PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00550GOTO 520
00560L2=O1
00570PRINT L$
00580PRINT "THIS MODULE COMPUTES AND PRINTS THE PROBABILITY THAT THE RATIO"
00590IF L1/SQR(V1+1)<L2/SQR(V2+1) THEN 650
00600B2=1
00610PRINT " -----------------"
00620PRINT " SIGMA-1 / SIGMA-2"
00630PRINT " -----------------"
00640GOTO 690
00650B2=2
00660PRINT " -----------------"
00670PRINT " SIGMA-2 / SIGMA-1"
00680PRINT " -----------------"
00690PRINT "IS GREATER THAN K, WHERE K CAN TAKE ON ANY POSITIVE VALUE."
00700PRINT "THIS RATIO WAS CHOSEN BECAUSE FOR THE MODAL ESTIMATES OF"
00710PRINT "THE STANDARD DEVIATIONS THIS RATIO IS GREATER THAN 1."
00720PRINT
00730PRINT "YOU CAN SPECIFY UP TO 6 DIFFERENT VALUES AT A TIME. AFTER"
00740PRINT "THE PROBABILITIES FOR THESE K VALUES HAVE BEEN COMPUTED AND "
00750PRINT "PRINTED YOU ARE GIVEN THE OPPORTUNITY TO SPECIFY MORE."
00760PRINT
00770PRINT "HOW MANY K VALUES DO YOU WANT TO SPECIFY";
00780GOSUB 9000
00790IF O1 >= 1 THEN 820
00800PRINT "REENTER. N MUST BE BETWEEN 1 AND 6 (INCLUSIVE)."
00810GOTO 780
00820IF O1>6 THEN 800
00830N=O1
00840FOR K5=1 TO N
00850REM
00860PRINT"INPUT VALUE ";K5;" ";
00870GOSUB 9000
00880IF O1>0 THEN 910
00890PRINT "REENTER. K MUST BE POSITIVE."
00900GOTO 870
00910L(K5)=O1
00920NEXT K5
00930PRINT L$
00940PRINT " COMPARISON OF TWO STANDARD DEVIATIONS"
00950PRINT " SIGMA-1 INVERSE CHI DISTRIBUTIONS SIGMA-2"
00960PRINT "============================================================"
00970:#####.## DEGREES OF FREEDOM ######.##
00980PRINT USING 970,V1,V2
00990:#####.## SCALE PARAMETER ######.##
01000PRINT USING 990,L1,L2
01010:#####.## MEAN ######.##
01020PRINT USING 1010,L1/SQR(V1-1.5),L2/SQR(V2-1.5)
01030:#####.## MODE ######.##
01040PRINT USING 1030,L1/SQR(V1+1),L2/SQR(V2+1)
01050PRINT "============================================================"
01060IF B2=2 THEN 1110
01070F7=(V1/(L1*L1))/(V2/(L2*L2))
01080A=V2/2
01090B=V1/2
01100GOTO 1120
01110F7=(V2/(L2*L2))/(V1/(L1*L1))
01120A=V1/2
01130B=V2/2
01140GOSUB 5220
01150FOR K5=1 TO N
01160K=L(K5)**2
01170J1=0
01180J2=(K*F7*A)/(B+K*F7*A)
01190GOSUB 5000
01200P=1-P
01210IF B2=1 THEN 1250
01220: PROB( SIGMA-2/SIGMA-1 > ###.## ) =##.##
01230PRINT USING 1220,L(K5),P
01240GOTO 1270
01250: PROB( SIGMA-1/SIGMA-2 > ###.## ) =##.##
01260PRINT USING 1250,L(K5),P
01270NEXT K5
01280PRINT "============================================================"
01290PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01300GOSUB 9000
01310PRINT L$
01320PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
01330PRINT " 1. FURTHER EVALUATE THIS RATIO"
01340PRINT " 2. EVALUATE A DIFFERENT RATIO"
01350PRINT " 3. EXIT MODULE"
01360GOSUB 9000
01370IF O1=3 THEN 1420
01380IF O1=2 THEN 1430
01390IF O1=1 THEN 1460
01400PRINT "REENTER. INPUT MUST BE THE NUMBER OF AN OPTION."
01410GOTO 1360
01420CHAIN "RSTRT"
01430B2=0
01440PRINT L$
01450GOTO 240
01460PRINT L$
01470GOTO 760
05000REM ***************************************************
05005REM BETA CDF ROUTINE
05010REM INPUT A B J1 J2
05015REM OUTPUT P
05020REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220
05025 DIM W(16),O(16)
05030IF J1 <> 0 THEN 5035
05031IF A<5 THEN 5035
05032IF B >= 5 THEN 5400
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A <= B THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=0 THEN 5115
05105IF D9=1 THEN 5115
05107D9=LOG(D9)*(A-1)+LOG(1-D9)*(B-1)
05108IF D9<-80 THEN 5115
05110P=P+W(I1)*EXP(D9)
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=0 THEN 5155
05135A=C6
05140B=C7
05145P=1-P
05150J2=1-J2
05155RETURN
05160FOR I1=1 TO 16
05165READ W(I1),O(I1)
05170NEXT I1
05175DATA 2.71525E-02,-.989401
05180DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
05185DATA .124629,-.755404,.149596,-.617876
05190DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
05195DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
05200DATA .149596,.617876,.124629,.755404
05205DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
05210DATA .989401
05215RETURN
05220G9=A+B
05225GOSUB 5850
05230F0=G0
05235G9=A
05240GOSUB 5850
05245F0=F0-G0
05250G9=B
05255GOSUB 5850
05260F0=F0-G0
05265IF A+B>85 THEN 5275
05270F0=EXP(F0)
05275RETURN
05280W1=(B*J2)**(1/3)
05285W2=(A*(1-J2))**(1/3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)**(1/3)
05305W2=(A*(1-J1))**(1/3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05340REM 2/16/76 CHANGED TO ALL LOG
05345D2=LOG(F0)+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350RETURN
05390REM
05400REM *******************************************************
05403REM PEIZER PRATT APROXIMATION
05406D0=.02*(J2/B-(1-J2)/A+(J2-.5)/(A+B))
05409D0=B-1/3-(A+B-2/3)*(1-J2)+D0
05412Y=(B-.5)/((A+B-1)*(1-J2))
05415GOSUB 5436
05418D9=G6
05421Y=(A-.5)/((A+B-1)*J2)
05424GOSUB 5436
05427Y3=D0*SQR((1+J2*D9+(1-J2)*G6)/((A+B-5/6)*(1-J2)*J2))
05430GOSUB 8000
05433RETURN
05436IF Y<.00001 THEN 5451
05439IF ABS(Y-1)<.00001 THEN 5457
05442G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05445IF G6>-1 THEN 5460
05446G6=-1
05448GOTO 5460
05451G6=1
05454GOTO 5460
05457G6=0
05460RETURN
05465REM END BETA CDF ROUTINE
05470REM***********************************************************
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/G5
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 ****************************************************
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
09999 END