Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmodf.bas
There are 2 other files named cmodf.bas in the archive. Click here to see a list.
00020REM
00030REM BEHRENS FISHER DISTRIBUTION
00040REM
00050REM**********************************************************
00060REM**** *
00070REM**** CMODF *
00080REM**** *
00090REM*********************************************************
00100 FILES RFILE1,RFILE2,RFILE3
00110X=0
00150RESTORE#1
00151 INPUT# 1,I1,I2,I3
00160SCRATCH#1
00161 PRINT # 1,15,I2,I3
00170RESTORE#2
00171 INPUT# 2,B1,B2,B3,B4,B5
00180S6=0
00190V8=0
00200V7=0
00210GOSUB 6060
00220IF B1=0 THEN 510
00230V8=1
00240V7=1
00250GOTO 510
00260PRINT L$
00270PRINT " BEHRENS-FISHER DISTRIBUTION"
00280PRINT
00290PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE CHARACTERISTICS OF A"
00300PRINT "BEHRENS-FISHER DISTRIBUTION. TYPE THE NUMBER OF THE OPTION"
00310PRINT "YOU WANT."
00320PRINT " 1. PERCENTILES"
00330PRINT " 2. HIGHEST DENSITY REGIONS"
00340PRINT " 3. PROBABILITY LESS THAN SOME VALUE"
00350PRINT " 4. PROBABILITY BETWEEN TWO VALUES"
00360PRINT " 5. GRAPH OF THE DENSITY FUNCTION"
00370PRINT " 6. EXIT"
00380GOSUB 9000
00390IF O1=6 THEN 470
00400IF O1=5 THEN 480
00410IF O1=1 THEN 480
00420IF O1=2 THEN 480
00430IF O1=3 THEN 480
00440IF O1=4 THEN 480
00450PRINT "REENTER. INPUT MUST BE NUMBER OF OPTION."
00460GOTO 380
00470CHAIN "RSTRT"
00480I1=O1
00490IF S6=1 THEN 1300
00500GOTO 1160
00510PRINT
00520IF V8=1 THEN 680
00530PRINT "INPUT THE PARAMETERS OF THE BEHRENS-FISHER DISTRIBUTION."
00540PRINT
00550PRINT "LET THE PARAMETERS OF THE TWO T DISTRIBUTIONS BE DENOTED BY"
00560PRINT " NU1 AND NU2 DEGREES OF FREEDOM"
00570PRINT " M1 AND M2 MEANS"
00580PRINT " K1 AND K2 SCALE PARAMETERS"
00590PRINT
00600PRINT "PSI=ARCTANGENT OF SQUARE ROOT (NU2 TIMES K1 DIVIDED BY"
00610PRINT " NU1 TIMES K2)";
00620GOSUB 9000
00630IF O1 >= 0 THEN 660
00640PRINT "REENTER. PSI MUST BE AT LEAST 0 AND NOT MORE THAN 90."
00650GOTO 600
00660IF O1>90 THEN 640
00670B3=O1
00680IF B3 <= 45 THEN 710
00690Z7=90-B3
00700GOTO 720
00710Z7=B3
00720IF V8=1 THEN 920
00730PRINT
00740PRINT "NU1=DEGREES OF FREEDOM OF 1ST T-DISTRIBUTION";
00750GOSUB 9000
00760IF O1 >= 6 THEN 810
00770PRINT
00780 REM
00790 PRINT"REENTER. DEGREES OF FREEDOM MUST BE AT LEAST 6."
00800GOTO 750
00810N0=O1
00820B1=O1
00830PRINT
00840PRINT "NU2=DEGREES OF FREEDOM OF THE 2ND DISTRIBUTION";
00850GOSUB 9000
00860IF O1 >= 6 THEN 890
00870PRINT "REENTER. MUST BE AT LEAST 6."
00880GOTO 830
00890B2=O1
00900N1=O1
00910IF B3 <= 45 THEN 940
00920N0=B2
00930N1=B1
00940Z7=Z7*1.74533E-02
00950IF V8=0 THEN 980
00960E0=B4
00970GOTO 1070
00980PRINT
00990PRINT "EPSILON=SQUARE ROOT OF K1 DIVIDED BY NU1 PLUS K2 DIVIDED BY NU2";
01000GOSUB 9000
01010IF O1>0 THEN 1050
01020PRINT
01030PRINT "REENTER. EPSILON MUST BE GREATER THAN 0."
01040GOTO 980
01050E0=O1
01060B4=E0
01070IF V8=0 THEN 1100
01080M0=B5
01090GOTO 1150
01100PRINT
01110PRINT "ZETA=MEAN OF 1ST MINUS MEAN OF 2ND";
01120GOSUB 9000
01130M0=O1
01140B5=M0
01150GOTO 260
01160K2=E0**2*N1/(((SIN(Z7)/COS(Z7))**2+1))
01170K1=(E0**2-(K2/N1))*N0
01180K1=K1/N0
01190K2=K2/N1
01200C1=K2/(K1+K2)
01210S1=1-C1
01220F1=(N1/(N1-2))*C1+(N0/(N0-2))*S1
01230F2=(N1**2/((N1-2)**2*(N1-4)))*C1**2+(N0**2/((N0-2)**2*(N0-4)))*S1**2
01240B=4+F1**2/F2
01250G=B
01260A=((B-2)/B)*F1
01270E0=(A*(K1+K2))**.5
01280L1=E0
01290S0=SQR(L1/(G-2))
01300PRINT L$
01310 ONI1 GOTO 1320,1610,1810,2010,2730
01320PRINT " OPTION 1: PERCENTILES"
01330PRINT
01340PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
01350PRINT "INPUT PERCENTILE AS NUMBER FROM 2.5 THROUGH 97.5."
01360GOSUB 2090
01370GOSUB 2380
01380GOSUB 2090
01390PRINT "INPUT PERCENTILE";
01400GOSUB 9000
01410IF O1=0 THEN 2120
01420IF O1<2.5 THEN 1590
01430IF O1>97.5 THEN 1590
01440IF O1=50 THEN 1480
01450IF O1>50 THEN 1510
01460P9=1-O1/100
01470GOTO 1520
01480B9=M0
01490P9=O1/100
01500GOTO 1560
01510P9=O1/100
01520GOSUB 2840
01530IF O1>50 THEN 1560
01540P9=O1/100
01550B9=2*M0-B9
01560: ##.## PERCENTILE =########.##
01570PRINT USING 1560,P9*100,B9
01580GOTO 1390
01590PRINT "REENTER. PERCENTILE MUST BE FROM 2.5 THRU 97.5."
01600GOTO 1390
01610PRINT " OPTION 2: HIGHEST DENSITY REGIONS"
01620PRINT
01630PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
01640PRINT "YOU CAN GET ANY P% HDR FROM 20 TO 98."
01650GOSUB 2090
01660GOSUB 2380
01670GOSUB 2090
01680PRINT "INPUT P%";
01690GOSUB 9000
01700IF O1=0 THEN 2120
01710IF O1<20 THEN 1790
01720IF O1>98 THEN 1790
01730P9=O1/200
01740P9=.5+P9
01750GOSUB 2840
01760: ##.#% HDR = (######.##,######.## )
01770PRINT USING 1760,O1,2*M0-B9,B9
01780GOTO 1680
01790PRINT "REENTER. MUST BE AT LEAST 20 AND NOT MORE THAN 98."
01800GOTO 1680
01810PRINT " OPTION 3: PROBABILITY LESS THAN SOME VALUE"
01820H9=1
01830PRINT "TO EXIT ROUTINE TYPE '-7777' WHEN ASKED FOR INPUT."
01840GOSUB 2090
01850GOSUB 2380
01860GOSUB 2090
01870PRINT "INPUT X";
01880GOSUB 9000
01890IF O1=-7777 THEN 2120
01900E2=(O1-M0)/E0
01910J2=E2
01920G=B
01930GOSUB 6000
01940: PROB ( BF <######.## ) =##.##
01950PRINT USING 1940,O1,P
01960: PROB ( BF >######.## ) =##.##
01970PRINT USING 1960,O1,1-P
01980GOTO 1870
01990GOTO 430
02000GOTO 1870
02010PRINT " OPTION 4: PROBABILITY BETWEEN TWO VALUES"
02020H9=1
02030PRINT
02040PRINT "TO EXIT ROUTINE TYPE '-7777'S WHEN ASKED FOR INPUT."
02050GOSUB 2090
02060GOSUB 2380
02070GOSUB 2090
02080GOTO 2510
02090 REM
02100PRINT"-----------------------------------------------------------"
02110RETURN
02120PRINT L$
02130S6=0
02140PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
02150PRINT " 1. FURTHER EVALUATE THIS DISTRIBUTION"
02160IF V7=0 THEN 2190
02170PRINT " 2. EXIT MODULE"
02180GOTO 2210
02190PRINT " 2. EVALUATE A DIFFERENT BEHRENS-FISHER DISTRIBUTION"
02200PRINT " 3. EXIT MODULE"
02210GOSUB 9000
02220IF V7=0 THEN 2270
02230IF O1=1 THEN 2330
02240IF O1=2 THEN 2320
02250PRINT "REENTER. INPUT MUST BE 1 OR 2."
02260GOTO 2210
02270IF O1=1 THEN 2330
02280IF O1=2 THEN 530
02290IF O1=3 THEN 2320
02300PRINT "REENTER. INPUT MUST BE 1,2 OR 3."
02310GOTO 2210
02320CHAIN "RSTRT"
02330S6=1
02340PRINT
02350V8=0
02360PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
02370GOTO 320
02380PRINT " BEHRENS-FISHER DISTRIBUTION"
02390PRINT "------------------------------------------------------"
02400:NU1=####.## NU2=####.## PSI=##.## DEGREES
02410PRINT USING 2400,B1,B2,B3
02420:EPSILON (SCALE) =#####.### ZETA (MEAN) =###.##
02430PRINT USING 2420,B4,B5
02440T1=(TAN(B3)**2*B4*B4)*B1/((TAN(B3)**2)+1)
02450T2=(B4**2-(1/(B1))*T1)/(1/B2)
02460T3=SQR((T1)/(B1-2)+T2/(B2-2))
02470:STANDARD DEVIATION=#####.###
02480PRINT USING 2470,T3
02490PRINT "------------------------------------------------------"
02500RETURN
02510PRINT "INPUT SMALLER VALUE";
02520GOSUB 9000
02530V4=O1
02540IF O1=-7777 THEN 2550
02550PRINT "INPUT LARGER VALUE";
02560GOSUB 9000
02570V5=O1
02580IF O1=-7777 THEN 2120
02590IF V4<V5 THEN 2620
02600PRINT "REENTER. ENTER SMALLER FIRST."
02610GOTO 2510
02620E2=(V4-M0)/E0
02630J2=E2
02640GOSUB 6000
02650F6=P
02660E2=(V5-M0)/E0
02670J2=E2
02680GOSUB 6000
02690F6=P-F6
02700: PROB (#####.## < X < #####.## ) =##.##
02710PRINT USING 2700,V4,V5,F6
02720GOTO 2510
02730REM OPTION 5 GRAPH
02740PRINT " OPTION 5: GRAPH OF DENSITY FUNCTION OVER 99% HDR"
02750GOSUB 2380
02760PRINT "THESE ARE THE PARAMETERS OF THE DISTRIBUTION TO BE GRAPHED."
02770J5=.99
02780GOSUB 7500
02790Q=SQR((G-2)/G)*S0*J2
02800K0=M0-Q
02810K1=M0+Q
02820GOSUB 9400
02830GOTO 2120
02840REM--*****ROUTINE TO CALC T PERCENTILE*****
02850E5=0
02860E6=2.5
02870E7=5
02880J2=E6
02890GOSUB 6000
02900IF ABS(P-P9) <= .00001 THEN 2980
02910IF P>P9 THEN 2950
02920E5=E6
02930E6=(E5+E7)/2
02940GOTO 2880
02950E7=E6
02960E6=(E5+E7)/2
02970GOTO 2880
02980B9=M0+E6*E0
02990RETURN
03000REM--*****END OF ROUTINE*****
05345G9=(1/2)*(G+1)
05350GOSUB 5850
05355D2=G0
05360G9=(1/2)*G
05365GOSUB 5850
05370D2=D2-G0
05375D2=D2+.5*G*LOG(L1)
05380D2=D2-.5*(G+1)*LOG(L1+(J2-M0)*(J2-M0))
05385RETURN
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 ****************************************************
06000REM****************************************************************
06002REM STUDENT'S T CDF ROUTINE
06004REM INPUT G J2
06006REM OUTPUT P
06008REM PRIOR GOSUB 6060
06010IF J2<6 THEN 6016
06012P=1
06014GOTO 6058
06016 DIM W(16),O(16)
06018Y3=J2
06020IF G <= 120 THEN 6026
06022GOSUB 8000
06024GOTO 6058
06026IF G=1 THEN 6098
06028J1=0
06030N=G
06032P=0
06034GOSUB 6084
06036D0=(J2-J1)*.5
06038D1=(J1+J2)*.5
06040FOR I1=1 TO 16
06042D9=D0*O(I1)+D1
06044IF D9=0 THEN 6050
06046IF D9=1 THEN 6050
06048P=P+W(I1)*(EXP(-(N+1)/2*LOG(1+D9*D9/N)))
06050NEXT I1
06052P=P*F0
06054P=P*D0
06056P=P+.5
06058RETURN
06060FOR I1=1 TO 16
06062READ W(I1),O(I1)
06064NEXT I1
06066DATA 2.71525E-02,-.989401
06068DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
06070DATA .124629,-.755404,.149596,-.617876
06072DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
06074DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
06076DATA .149596,.617876,.124629,.755404
06078DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
06080DATA .989401
06082RETURN
06084G9=(N+1)/2
06086GOSUB 5850
06088F0=G0
06090G9=N/2
06092GOSUB 5850
06094F0=EXP(F0-G0)/SQR(3.14159*N)
06096RETURN
06098REM FOLLOWING FOR NU=1
06100P=.5+1/3.14159*ATN(Y3)
06102RETURN
06104REM END OF STUDENT'S T CDF ROUTINE
06106REM*************************************************************
07500REM ************************************************************
07501REM STUDENT'S T DISTRIBUTION HIGHEST DENSITY REGIONS
07502REM INPUTS G J5
07503REM J2
07504REM
07505Z8=.5
07506N=G
07507X9=1
07508J1=0
07509J2=X9
07510GOSUB 6000
07511P=2*P-1
07512Z9=P
07513IF P>J5 THEN 7517
07514X9=X9+1
07515Z8=Z9
07516GOTO 7508
07517X0=X9-1
07518X2=X9
07519X9=X0+(J5-Z8)*(X2-X0)/(Z9-Z8)
07520J1=0
07521J2=X9
07522GOSUB 6000
07523P=2*P-1
07524IF ABS(P-J5)<.0001 THEN 7534
07525IF P<J5 THEN 7530
07526X2=X9
07527Z9=P
07528X9=X0+(J5-Z8)*(X2-X0)/(Z9-Z8)
07529GOTO 7520
07530X0=X9
07531Z8=P
07532X9=X0+(J5-Z8)*(X2-X0)/(Z9-Z8)
07533GOTO 7520
07534J2=X9
07535RETURN
07536REM
07537REM END OF STUDENT'S T HDR ROUTINE
07538REM *********************************************************
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)
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******************************************************
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
09471: GRAPH OF BEHRENS-FISHER ##.# % HDR
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 9471,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
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 WAN TO CONTINUE TYPE '1'.";
09587GOTO 9555
09999 END