Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod2.bas
There are 2 other files named cmod2.bas in the archive. Click here to see a list.
00030REM CMOD2
00120FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00160RESTORE#1
00161 INPUT# 1,I1,I2,I3
00170SCRATCH#1
00171 PRINT # 1,2,I2,I3
00180S9=0
00190V7=0
00200V2=1
00210S8=0
00220GOSUB 5160
00230 DIM S(4,5)
00240PRINT L$
00250PRINT " PRIOR DISTRIBUTION - BETA-BINOMIAL MODEL"
00260PRINT
00270PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A BETA DISTRIBUTION"
00280PRINT "TO YOUR PRIOR BELIEFS ABOUT PI. WE BEGIN BY ASKING YOU TO"
00290PRINT "SPECIFY THE 25TH, 50TH AND 75TH PERCENTILES OF YOUR PRIOR"
00300PRINT "DISTRIBUTION."
00310PRINT
00320PRINT "SPECIFY 25TH. YOUR BETTING ODDS ARE 3 TO 1 THAT PI IS"
00330PRINT "GREATER THAN THIS VALUE.";
00340GOSUB 9000
00350IF O1 <= 0 THEN 370
00360IF O1<.97 THEN 390
00370PRINT "REENTER. 25TH MUST BE GREATER THAN 0 AND LESS THAN .97."
00380GOTO 340
00390Q1=O1
00400PRINT
00410PRINT "SPECIFY 50TH. YOUR BETTING ODDS ARE EVEN THAT PI IS"
00420PRINT "GREATER THAN THIS VALUE.";
00430GOSUB 9000
00440IF O1>Q1 THEN 470
00450PRINT "REENTER. 50TH MUST BE GREATER THAN 25TH."
00460GOTO 430
00470IF O1 >= .03 THEN 500
00480PRINT "REENTER. 50TH MUST BE AT LEAST .03."
00490GOTO 430
00500IF O1 <= .97 THEN 530
00510PRINT "REENTER. 50TH MUST NOT BE GREATER THAN .97."
00520GOTO 430
00530Q2=O1
00540PRINT
00550PRINT "SPECIFY 75TH. YOUR BETTING ODDS ARE 1 TO 3 THAT PI IS"
00560PRINT "GREATER THAN THIS VALUE.";
00570GOSUB 9000
00580IF O1>Q2 THEN 610
00590PRINT "REENTER. 75TH MUST BE GREATER THAN 50TH AND LESS THAN 1."
00600GOTO 570
00610IF O1 >= 1 THEN 590
00620Q3=O1
00650REM Q1 25TH PERCENTILE
00660REM Q2 50TH PERCENTILE
00670REM Q3 75TH PERCENTILE
00680PRINT
00690PRINT "POSSIBLE APPROXIMATE DISTRIBUTIONS ARE BEING COMPUTED."
00700X=SQR(Q2*(1-Q1))-SQR(Q1*(1-Q2))
00710X=X*X
00720C=.114/X
00730X7=1/X
00740A=C*Q2+1/3
00750B=C*(1-Q2)+1/3
00760K5=1
00770GOSUB 1660
00810IF ABS((S(1,2)-Q1))<.01 THEN 830
00820GOTO 900
00830IF ABS((S(1,3)-Q2))<.01 THEN 850
00840GOTO 900
00850IF ABS((S(1,4)-Q3)) >= .01 THEN 900
00860GOTO 2590
00900PRINT L$
00910PRINT "HERE ARE SOME OF THE PERCENTILES OF FOUR BETA DISTRIBUTIONS"
00920PRINT "THAT HAVE BEEN FITTED TO YOUR PERCENTILE SPECIFICATIONS."
00930PRINT
00940PRINT " 10TH 25TH 50TH 75TH 90TH"
00950GOSUB 2560
00960X=SQR(Q2*(1-Q3))-SQR(Q3*(1-Q2))
00970X=X*X
00980C=.114/X
00990A=C*Q2+1/3
01000B=C*(1-Q2)+1/3
01010K5=2
01020GOSUB 1660
01030C2=.057*(1/X+X7)
01040A=C2*Q2+1/3
01050B=C2*(1-Q2)+1/3
01060K5=3
01070GOSUB 1660
01080A=5*P(1,1)+3*P(2,1)+P(3,1)
01090A=A/9
01100B=5*P(1,2)+3*P(2,2)+P(3,2)
01110B=B/9
01120K5=4
01130GOSUB 1660
01140PRINT
01150GOSUB 4780
01160GOSUB 9000
01170IF O1=0 THEN 1190
01180GOTO 1500
01190PRINT L$
01200PRINT "RESPECIFY PERCENTILES."
01210PRINT
01220GOTO 320
01230GOSUB 9000
01240IF O1>0 THEN 1270
01250PRINT "REENTER. 25TH MUST BE GREATER THAN 0 AND LESS THAN 1."
01260GOTO 1230
01270IF O1 >= 1 THEN 1250
01280Q1=O1
01290PRINT "50TH";
01300GOSUB 9000
01310IF O1>Q1 THEN 1340
01320PRINT "REENTER. 50TH MUST BE GREATER THAN 25TH AND LESS THAN 1."
01330GOTO 1300
01340IF O1 >= 1 THEN 1320
01350Q2=O1
01360PRINT "75TH";
01370GOSUB 9000
01380IF O1>Q2 THEN 1410
01390PRINT "REENTER. 75TH MUST BE GREATER THAN 50TH AND LESS THAN 1."
01400GOTO 1360
01410IF O1 >= 1 THEN 1390
01420Q3=O1
01430GOTO 700
01440REM*************************************************************
01450REM
01460REM USER DOES NOT WANT TO RESPECIFY PERCENTILES
01470REM PICKS ONE OF FOUR APPROXIMATE DISTRIBUTIONS AS
01480REM TENTATIVE FIT
01490REM
01500IF O1=1 THEN 1570
01510IF O1=2 THEN 1570
01520IF O1=3 THEN 1570
01530IF O1=4 THEN 1570
01540PRINT "REENTER. MUST BE 0 OR NUMBER OF A DISTRIBUTION."
01550GOSUB 9000
01560GOTO 1170
01570K5=O1
01580A=P(K5,1)
01590B=P(K5,2)
01600GOTO 4150
01610REM**********************************************************
01620REM
01630REM IF EITHER PARAMETER IS NOT GREATER THAN 1.15
01640REM ADJUSTMENTS ARE MADE TO BOTH SO THAT BOTH ARE
01650REM
01660IF A>B THEN 2230
01670IF A>1.15 THEN 2310
01680IF K5=4 THEN 1880
01690IF K5=3 THEN 1800
01700REM A IS LESS THAN 1.15
01710REM
01720REM SET A=1.15*(K5-1)*.4 AND LET Q2 BE AN ESTIMATE OF MEAN
01730REM
01740B=B*(1.15+(K5-1)*.4)/A
01750A=1.15+(K5-1)*.4
01760GOTO 2310
01770REM
01780REM SET A=1.15 AND LET Q2 BE ESTIMATE OF MODE
01790REM
01800B=.15/(Q2)+.85
01810A=1.15
01820GOTO 2310
01830REM
01840REM SET M0=A+B WHERE B IS ESTIMATED TWICE-ONCE ASSUMING
01850REM Q2=MEAN AND AGAIN ASSUMING Q2=MODE KEEP UPPING M0
01860REM UNTIL BOTH PARAMETERS ARE GREATER THAN 1.15
01870REM
01880A=1.15
01890M0=A/Q2
01900M0=.5*M0+.5*A/Q2+1
01910B=M0-A
01920IF B>1.15 THEN 2310
01930A=A+.1
01940GOTO 1890
01950REM******************************************************
01960REM
01970REM B IS LESS THAN 1.15
01980REM
01990IF K5=3 THEN 2100
02000IF K5=4 THEN 2200
02010REM
02020REM SET B=1.15 AND LET Q2 BE ESTIMATE OF MEAN
02030A=A*(1.15+(K5-1)*.4)/B
02040B=1.15+(K5-1)*.4
02050GOTO 2310
02060REM
02070REM SET M0=A+B WHERE IT IS ESTIMATED TWICE-ONCE UNDER
02080REM ASSUMPTION Q2=MEAN AND AGAIN ASSUMING Q2=MODE AVERAGE
02090REM
02100B=1.15
02110M0=B/(1-Q2)
02120M0=.5*M0+.5*((2*Q2-B-1)/(Q2-1))
02130A=M0-B
02140IF A>1.15 THEN 2310
02150B=B+.1
02160GOTO 2100
02170REM
02180REM SET B=1.15 AND LET Q2=MODE
02190REM
02200B=1.15
02210A=(1-.85*Q2)/(1-Q2)
02220GOTO 2310
02230IF B<1.15 THEN 1990
02240IF A+B<2000 THEN 2270
02250A=Q2*2000
02260B=2000-A
02270REM *********************************************************
02280REM
02290REM COMPUTE PERCENTILES
02300REM
02310P1=.1
02320P(K5,1)=A
02330P(K5,2)=B
02340GOSUB 4500
02350S(K5,1)=J2
02360P1=.9
02370GOSUB 4530
02380S(K5,5)=J2
02390P1=.75
02400E1=S(K5,1)
02410A7=.1
02420GOSUB 4550
02430S(K5,4)=J2
02440P1=.25
02450E1=S(K5,1)
02460A7=.1
02470GOSUB 4550
02480S(K5,2)=J2
02490P1=.5
02500E2=S(K5,4)
02510A8=.75
02520GOSUB 4550
02530S(K5,3)=J2
02540IF S9=1 THEN 2620
02550IF K5=1 THEN 2580
02560PRINT USING 2570,K5,S(K5,1),S(K5,2),S(K5,3),S(K5,4),S(K5,5)
02570: ##. ##.## ##.## ##.## ##.## ##.##
02580RETURN
02590PRINT L$
02600PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE BETA"
02610PRINT "DISTRIBUTION FITTED TO YOUR PERCENTILE SPECIFICATIONS."
02620PRINT
02630M=A+B
02640: HYPOTHETICAL SAMPLE SIZE (M)######.##
02650: STANDARD DEVIATION ##.##
02660: 10TH PERCENTILE ##.##
02670: 25TH PERCENTILE ##.##
02680: 50TH (MEDIAN) ##.##
02690: 75TH PERCENTILE ##.##
02700: 90TH PERCENTILE ##.##
02710: 50% HDR ##.## -##.##
02720: 75% HDR ##.## -##.##
02730: MODE ##.##
02740: 95% HDR ##.## -##.##
02750: MEAN ##.##
02760: PARAMETER A ######.##
02770: PARAMETER B ######.##
02780J5=.5
02790J8=S(K5,4)/S(K5,2)
02800GOSUB 7000
02810IF V7=1 THEN 2840
02820 PRINT USING 2640,A+B
02830IF V7=0 THEN 2850
02840 PRINT USING 2730,(A-1)/(M-2)
02850 PRINT USING 2660,S(K5,1)
02860 PRINT USING 2670,S(K5,2)
02870 PRINT USING 2680,S(K5,3)
02880 PRINT USING 2690,S(K5,4)
02890 PRINT USING 2700,S(K5,5)
02900IF V7=1 THEN 2940
02910IF J2<.99 THEN 2930
02920J2=.99
02930GOTO 2960
02940J1=H(1)
02950J2=H(2)
02960 PRINT USING 2710,J1,J2
02970IF V7=1 THEN 3050
02980H(1)=J1
02990H(2)=J2
03000J5=.75
03010GOSUB 7000
03020IF J2<.99 THEN 3040
03030J2=.99
03040GOTO 3070
03050J1=H(3)
03060J2=H(4)
03070 PRINT USING 2720,J1,J2
03080IF V7=1 THEN 3160
03090H(3)=J1
03100H(4)=J2
03110J5=.95
03120GOSUB 7000
03130IF J2<.99 THEN 3150
03140J2=.99
03150GOTO 3180
03160J1=H(5)
03170J2=H(6)
03180 PRINT USING 2740,J1,J2
03190IF V7=1 THEN 4200
03200H(5)=J1
03210H(6)=J2
03220IF S9=1 THEN 3290
03230PRINT "IF YOU DO NOT FEEL THAT THE HYPOTHETICAL SAMPLE SIZE (M)"
03240PRINT "REFLECTS YOUR PRIOR INFORMATION ABOUT PI YOU CAN SPECIFY A"
03250PRINT "DIFFERENT VALUE FOR M. THIS WILL NOT AFFECT THE MEDIAN BUT"
03260PRINT "WILL CHANGE THE HDRS AND OTHER PERCENTILES. A LARGER M WILL"
03270PRINT "RESULT IN SHORTER INTERVALS, AND A SMALLER M IN LONGER ONES."
03280PRINT
03290IF V2=0 THEN 3410
03300V2=0
03310M2=3
03320M7=3*S(K5,3)*M2*(M2-2)+M2
03330M7=M7/(3*M2-4)
03340IF M7<1.5 THEN 3360
03350IF M2-M7>1.5 THEN 3380
03360M2=M2+1
03370GOTO 3320
03380M7=INT(M2+1)
03390IF M7<M THEN 3420
03400M7=M
03410:IF YOU WANT TO CHANGE M TYPE NEW VALUE (MIN =###.####).
03420 PRINT USING 3410,M7
03430PRINT "IF YOU DO NOT WANT TO CHANGE M TYPE '0'";
03440GOSUB 9000
03450IF O1=0 THEN 3610
03460IF O1>9999 THEN 3500
03470IF O1 >= M7 THEN 3520
03480PRINT "REENTER. INPUT MUST BE 0 OR AN ACCEPTABLE M VALUE."
03490GOTO 3440
03500PRINT "REENTER. M MUST NOT BE GREATER THAN 9999."
03510GOTO 3440
03520M=O1
03530IF ABS(S(K5,3)*M-S(K5,3)*(M-2)-1)>.019 THEN 3590
03540A=3*S(K5,3)*M*(M-2)+M
03550A=A/(3*M-4)
03560B=M-A
03570S9=1
03580GOTO 4050
03590O1=S(K5,3)
03600GOTO 3780
03610PRINT
03620IF S8=0 THEN 3640
03630GOTO 3670
03640PRINT "TO CHANGE THE CENTERING OF THE DISTRIBUTION, SPECIFY"
03650S8=1
03660PRINT "A DIFFERENT MEDIAN. THIS WILL NOT AFFECT THE VALUE OF M."
03670REM
03680PRINT "IF YOU WANT TO CHANGE MEDIAN TYPE NEW VALUE ELSE '0'.";
03690GOSUB 9000
03700IF O1=0 THEN 4415
03710IF O1<.03 THEN 3730
03720IF O1 <= .97 THEN 3750
03730PRINT "REENTER. MEDIAN MUST BE AT LEAST .03 AND NOT MORE THAN .97."
03740GOTO 3690
03750V2=1
03760IF O1<1 THEN 3780
03770GOTO 3730
03780L4=O1*M
03790L3=O1*(M-2)+1
03800IF L4>L3 THEN 3840
03810A=L3
03820L3=L4
03830L4=A
03840A=.5*(L3+L4)
03850B=M-A
03860GOTO 3930
03870L4=A
03880L3=A0
03890GOTO 3930
03900L4=A0
03910L3=A
03920A=.5*(L3+L4)
03930S9=1
03940P1=.5
03950B=M-A
03960GOSUB 4500
03970IF ABS(J2-O1)<.004 THEN 4050
03980IF J2>O1 THEN 4020
03990L3=A
04000A=.5*(L4+A)
04010GOTO 3950
04020L4=A
04030A=.5*(L3+A)
04040GOTO 3950
04050IF A<1.1 THEN 4070
04060IF B>1.1 THEN 4150
04070PRINT L$
04080PRINT "THE MEDIAN YOU SPECIFIED IS NOT CONSISTENT WITH THE M VALUE"
04090PRINT "YOU ARE USING. EXTREME VALUES (HIGH OR LOW) REQUIRE LARGER"
04100PRINT "M VALUES. YOU SHOULD EITHER MODERATE YOUR MEDIAN ESTIMATE"
04110PRINT "OR INCREASE THE M. WHEN YOU WANT TO CONTINUE TYPE '1'."
04120GOSUB 9000
04130V2=0
04140GOTO 1580
04150PRINT L$
04160PRINT "HERE ARE SOME CHARACTERISTICS OF THE BETA DISTRIBUTION YOU"
04170PRINT "ARE NOW CONSIDERING."
04180IF S9=0 THEN 2620
04190GOTO 2310
04200REM
04210PRINT "THIS COMPLETES THE FITTING OF A PRIOR DISTRIBUTION ON PI."
04220PRINT "IF YOU ARE NOT GOING TO CONTINUE WITH THE ANALYSIS AT THIS"
04230PRINT "TIME YOU SHOULD RECORD THE PARAMETERS OF YOUR PRIOR."
04240PRINT
04250PRINT "IF YOU WANT TO CONTINUE THE ANALYSIS TYPE '1', ELSE '0'.";
04260GOSUB 9000
04270IF O1=1 THEN 4330
04280IF O1=0 THEN 4320
04290PRINT
04300PRINT "REENTER. INPUT MUST BE 0 OR 1."
04310GOTO 4260
04320CHAIN "RSTRT"
04330SCRATCH#2
04331 PRINT # 2,A,B
04340PRINT L$
04350PRINT "IF YOU WANT TO DO A PREPOSTERIOR ANALYSIS TYPE '1', ELSE '0'.";
04360GOSUB 9000
04370IF O1=1 THEN 4410
04380IF O1=0 THEN 4400
04390PRINT "REENTER. INPUT MUST BE 0 OR 1."
04400CHAIN "CMOD3"
04410CHAIN "CMODAB"
04415J5=.99
04417GOSUB 7000
04418PRINT L$
04420M0=(A-1)/(A+B-2)
04421K0=J1
04422K1=J2
04423GOSUB 9400
04424PRINT "IF YOU WISH TO CHANGE YOUR PRIOR TYPE '1', ELSE '0'.";
04425GOSUB 9000
04427IF O1=1 THEN 4150
04429V7=1
04430PRINT L$
04440PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE PRIOR"
04450PRINT "DISTRIBUTION FITTED TO YOUR PRIOR BELIEFS ABOUT PI."
04460PRINT
04470 PRINT USING 2760,A
04480 PRINT USING 2770,B
04490GOTO 2640
04500E1=0
04510A7=0
04520GOSUB 5220
04530E2=1
04540A8=.999999
04550A6=P1
04560S5=(A6-A7)/(A8-A7)
04570IF S5<.85 THEN 4590
04580S5=.85
04590IF S5>.15 THEN 4610
04600S5=.15
04610J2=E1+S5*(E2-E1)
04620J1=0
04630GOSUB 5000
04640IF ABS(E1-E2)<.0001 THEN 4750
04650IF ABS(P-P1)<.005 THEN 4750
04660IF J2>.99 THEN 4760
04670IF P>P1 THEN 4710
04680E1=J2
04690A7=P
04700GOTO 4560
04710A8=P
04720E2=J2
04730GOTO 4560
04740GOTO 4750
04750IF J2<.99 THEN 4770
04760J2=.99
04770RETURN
04780PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE"
04790PRINT "WHICH MOST CLOSELY CORRESPONDS TO YOUR PRIOR BELIEFS."
04800PRINT "YOU CAN EITHER TENTATIVELY ACCEPT THIS DISTRIBUTION OR"
04810PRINT "RESPECIFY THE PERCENTILES."
04820PRINT
04830PRINT "IF YOU WANT ONE OF THESE DISTRIBUTIONS TYPE ITS NUMBER."
04840PRINT "IF YOU WANT TO RESPECIFY THE PERCENTILES TYPE '0'."
04850RETURN
05000REM BETA CDF
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
05345D2=LOG(F0)+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350RETURN
05390REM
05400REM PEIZER-PRATT APPROX
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
05600REM BETA PER'TILE
05615J1=0
05618P2=P1/100+.02
05619IF P2<1 THEN 5621
05620P2=.999
05621GOSUB 5762
05624GOSUB 5000
05627IF ABS(P1/100-P)<.0001 THEN 5750
05630IF P1/100>P THEN 5639
05633E4=J2
05636P8=P
05637A7=P2-.02
05638GOTO 5669
05639E1=J2
05642P7=P
05643A8=P2+.02
05645IF A8<1 THEN 5651
05648A8=.9999
05651P0=A8
05654P2=P0
05657GOSUB 5762
05658GOSUB 5000
05659IF P1/100 <= P THEN 5665
05660J2=J2+.03
05661IF J2<.999 THEN 5658
05662J2=.99991
05663GOTO 5658
05665E4=J2
05666P8=P
05667GOTO 5693
05669IF A7>0 THEN 5675
05672A7=.00001
05675P0=A7
05678P2=P0
05681GOSUB 5762
05682GOSUB 5000
05683IF P <= P1/100 THEN 5690
05684J2=J2-.03
05685IF J2>.001 THEN 5682
05686J2=.0001
05687GOTO 5682
05690P7=P
05691E1=J2
05693E2=E4
05696J2=(P1/100-P7)/(P8-P7)
05699IF J2>.05 THEN 5708
05705J2=.05
05706GOTO 5711
05708IF J2<.85 THEN 5711
05709J2=.85
05711J2=E1+J2*(E2-E1)
05714J1=0
05717GOSUB 5000
05720IF ABS(E1-E2)<.0001 THEN 5750
05726IF ABS(P-P1/100)<.0001 THEN 5750
05729IF P>P1/100 THEN 5741
05732E1=J2
05735P7=P
05738GOTO 5696
05741E2=J2
05744P8=P
05747GOTO 5696
05750RETURN
05762REM-INVERSE BETA
05777P3=P2
05780IF P2 <= .5 THEN 5786
05783P2=1-P2
05786A1=SQR(LOG(1/P2/P2))
05789A2=2.51552+.802853*A1+.010328*A1*A1
05792A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
05795J2=A1-A2
05798IF P3 <= .5 THEN 5810
05801P2=P3
05804J2=-J2
05807GOTO 5813
05810REM
05813L=J2*J2-3
05816L=L/6
05819H=2/(1/(2*A-1)+1/(2*B-1))
05822D=1/(2*B-1)-2/(2*A-1)
05825W=J2*(H+L) ** .5/H-D*(L+5/6-2/H/3)
05828X4=A/(A+B*EXP(2*W))
05831J2=X4
05837RETURN
05840REM********************* END INVERSE BETA **********************
05850REM LOG GAMMA
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
07000REM**************************************************************
07002REM BETA HDR ROUTINE
07004REM INPUT A B J5
07006REM OUTPUT J1 J2
07008J7=0
07010U0=0
07012IF A<B THEN 7024
07014U0=1
07016J7=A
07018A=B
07020J=0
07022B=J7
07024J6=SQR((A*B)/(A+B)/(A+B)/(A+B+1))
07026J=(A-1)/(B-1)
07028J0=(A-1)/(A+B-2)
07030GOSUB 5220
07032P9=0
07034J9=J0
07036J6=1.6*J6
07038J8=A/(A+B)+J5*J6
07040IF J8<1 THEN 7044
07042J8=.9999
07044IF A/(A+B)-J5*J6>0 THEN 7050
07046J8=J8/1.E-08
07048GOTO 7052
07050J8=J8/(A/(A+B)-J5*J6)
07052GOSUB 7168
07054GOSUB 5000
07056P0=P
07058IF ABS(P0-J5)<.0001 THEN 7180
07060IF P0>J5 THEN 7084
07062J0=J1
07064J9=J2
07066P9=P
07068J8=J8*1.5
07070GOSUB 7168
07072GOSUB 5000
07074IF P<J5 THEN 7068
07076J6=J1
07078J7=J2
07080P0=P
07082GOTO 7110
07084J6=J1
07086J7=J2
07088P0=P
07090J8=(1+J8)/2
07092GOSUB 7168
07094GOSUB 5000
07096IF P>J5 THEN 7090
07098J0=J1
07100J9=J2
07102P9=P
07104GOTO 7110
07106J6=J1
07108J7=J2
07110J8=(J5-P9)/(P0-P9)
07112J2=J6
07114GOSUB 5345
07116J1=EXP(D2)
07118J2=J0
07120GOSUB 5345
07122C=J8
07123IF J1<1.E-14 THEN 7126
07124J1=(-2*J1+SQR(4*J1*J1+4*C*(J2+J1)*(J2-J1)))/(2*(J2-J1))
07125GOTO 7127
07126J1=(-2*J1+SQR(4*C*(J2+J1)*(J2-J1)))/(2*(J2-J1))
07127IF J1<.9 THEN 7132
07128J1=.9
07130GOTO 7136
07132IF J1>.1 THEN 7136
07134J1=.1
07136J8=(J1*J7+(1-J1)*J9)/(J1*J6+(1-J1)*J0)
07138GOSUB 7168
07140GOSUB 5000
07142J3=P
07144IF ABS(J1-J2)<.0001 THEN 7180
07146IF ABS(J3-J5)>.0001 THEN 7150
07148GOTO 7180
07150IF J3>J5 THEN 7160
07152J0=J1
07154J9=J2
07156P9=J3
07158GOTO 7110
07160J6=J1
07162J7=J2
07164P0=J3
07166GOTO 7110
07168J1=(J8**J-1)/(J8**(J+1)-1)
07170J2=J1*J8
07172IF J2<1 THEN 7178
07174J8=J8*.95
07176GOTO 7172
07178RETURN
07180IF U0=0 THEN 7194
07182J7=A
07184A=B
07186B=J7
07188J7=J2
07190J2=1-J1
07192J1=1-J7
07194RETURN
07196REM END OF BETA HDR ROUTINE
07198REM*************************************************
08000REM NORMAL CDF
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
09000REM 1 INPUT
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09272PRINTUSING9466,J2,MID$(U$,1,K7)
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 YOUR PRIOR BETA ##.#% HDR
09426GOTO 9520
09450:#####.## I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
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
09462PRINT USING 9466,J2,MID$(T$,1,K7)
09464GOTO 9474
09466:#####.## I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
09468PRINTUSING 9466,J2,MID$(S$,1,K7)
09470GOTO 9474
09472PRINTUSING9466,J2,MID$(U$,1,K7)
09474RETURN
09510REM ***************************************************************
09520K2=(K1-K0)/19
09522J2=M0
09523 PRINT USING 9423,J5*100
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
09548PRINTUSING9450,K1,MID$(U$,1,K7)
09556RETURN
09560IF J2>M0+.4*K2 THEN 9572
09561IF J2<M0-.4*K2 THEN 9566
09562PRINTUSING9566,J2,MID$(T$,1,K7)
09564RETURN
09566:##.#^^^^ I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
09568PRINT USING 9566,J2,MID$(S$,1,K7)
09570RETURN
09572PRINT USING 9566,J2,MID$(U$,1,K7)
09574RETURN
09575PRINTUSING 9580,K1,MID$(U$,1,K7)
09577GOTO 9556
09580:##.#^^^^ I'LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
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 WANT TO CONTINUE TYPE '1'.";
09587GOTO 9556
09999END