Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod8.bas
There are 2 other files named cmod8.bas in the archive. Click here to see a list.
00020REM ***********************************************************************
00030REM CMOD8 CMOD8 CMOD8 CMOD8 CMOD8 CMOD8
00040REM **********************************************************************
00050REM *********************************************************************
00060REM
00070REM SPECIFYING TYPICAL PRIOR M-GROUP PROPORTIONS
00080REM
00090REM **********************************************************************
00100 DIM P(5,5)
00110FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00150GOSUB 5160
00160RESTORE#1
00161 INPUT# 1,I1,I2,I3
00170SCRATCH#1
00171 PRINT # 1,8,I2,I3
00180PRINT L$
00190X=0
00200S9=0
00210S8=0
00220V7=0
00230V2=0
00240 DIM S(4,5)
00250PRINT L$
00260PRINT " SPECIFICATION OF A PRIOR FOR SIMULTANEOUS PROPORTIONS"
00270PRINT
00280PRINT "THIS MODULE WILL ASSIST YOU IN SPECIFYING A PRIOR DISTRIBUTION"
00290PRINT "FOR YOUR ANALYSIS. THE ASSUMPTION IS MADE THAT YOUR INFORMATION"
00300PRINT "ABOUT THE GROUPS YOU PLAN TO OBSERVE (OR HAVE OBSERVED) - AND "
00310PRINT "ABOUT ANY OTHER GROUPS YOU MIGHT OBSERVE - IS EXCHANGEABLE."
00320PRINT "THUS IT IS ASSUMED YOUR INFORMATION ABOUT A GIVEN GROUP IS "
00330PRINT "EQUIVALENT TO YOUR INFORMATION ABOUT ANY OTHER GROUP AS CONCERNS"
00340PRINT "THE PROBABILITY OF SUCCESS IN EACH OF THESE GROUPS."
00350PRINT
00360PRINT "CONSIDER A PARTICULAR GROUP UNDER STUDY. WE NEED TO ASSESS"
00370PRINT "YOUR KNOWLEDGE CONCERNING THE PROBABILITY OF SUCCESS (PI) FOR"
00380PRINT "THIS GROUP."
00390PRINT
00400PRINT "WE BEGIN BY ASKING YOU TO SPECIFY THE 25TH, 50TH AND 75TH"
00410PRINT "PERCENTILES OF YOUR PRIOR DISTRIBUTION FOR PI."
00420PRINT
00430PRINT
00440PRINT "SPECIFY 25TH. YOUR BETTING ODDS ARE 3 TO 1 THAT PI IS GREATER"
00450PRINT "THAN THIS VALUE.";
00460GOSUB 9000
00470IF O1>0 THEN 500
00480PRINT "REENTER. 25TH MUST BE GREATER THAN 0."
00490GOTO 460
00500IF O1<.97 THEN 530
00510PRINT "REENTER. 25TH MUST BE LESS THAN .97."
00520GOTO 460
00530Q1=O1
00540PRINT
00550PRINT "SPECIFY 50TH. YOUR BETTING ODDS ARE EVEN THAT PI IS GREATER"
00560PRINT "THAN THIS VALUE.";
00570GOSUB 9000
00580IF O1>Q1 THEN 610
00590PRINT "REENTER. 50TH MUST BE GREATER THAN 25TH."
00600GOTO 570
00610IF O1 >= .03 THEN 640
00620PRINT "REENTER. 50TH MUST NOT BE LESS THAN .03."
00630GOTO 570
00640IF O1 <= .97 THEN 670
00650PRINT "REENTER. 50TH MUST NOT BE GREATER THAN .97."
00660GOTO 570
00670Q2=O1
00680PRINT
00690PRINT "SPECIFY 75TH. YOUR BETTING ODDS ARE 3 TO 1 THAT PI IS LESS THAN"
00700PRINT "THIS VALUE.";
00710GOSUB 9000
00720IF O1>Q2 THEN 750
00730PRINT "REENTER. 75TH MUST BE GREATER THAN 50TH."
00740GOTO 710
00750IF O1<1 THEN 780
00760PRINT "REENTER. 75TH MUST BE LESS THAN 1."
00770GOTO 710
00780Q3=O1
00790PRINT
00800REM ******************************************************************
00810REM
00820PRINT "POSSIBLE APPROXIMATE DISTRIBUTIONS ARE BEING COMPUTED."
00830REM
00840REM ******************************************************************
00850X=SQR(Q2*(1-Q1))-SQR(Q1*(1-Q2))
00860X=X*X
00870C=.114/X
00880X7=1/X
00890A=C*Q2+1/3
00900B=C*(1-Q2)+1/3
00910K5=1
00920GOSUB 1660
00930IF ABS((S(1,2)-Q1))<.01 THEN 950
00940GOTO 990
00950IF ABS((S(1,3)-Q2))<.01 THEN 970
00960GOTO 990
00970IF ABS((S(1,4)-Q3)) >= .01 THEN 990
00980GOTO 2060
00990PRINT L$
01000PRINT "HERE ARE THE PERCENTILES OF FOUR BETA DISTRIBUTIONS WHICH HAVE"
01010PRINT "BEEN FITTED TO YOUR PERCENTILE SPECIFICATIONS."
01020PRINT
01030PRINT" 10TH 25TH 50TH 75TH 90TH"
01040GOSUB 2010
01050X=SQR(Q2*(1-Q3))-SQR(Q3*(1-Q2))
01060X=X*X
01070C=.114/X
01080A=C*Q2+1/3
01090B=C*(1-Q2)+1/3
01100K5=2
01110GOSUB 1660
01120C2=.057*(1/X+X7)
01130A=C2*Q2+1/3
01140B=C2*(1-Q2)+1/3
01150K5=3
01160GOSUB 1660
01170A=P(1,1)+P(2,1)+P(3,1)
01180A=A/3
01190B=P(1,2)+P(2,2)+P(3,2)
01200B=B/3
01210K5=4
01220GOSUB 1660
01230REM ********************************************************
01240REM
01250REM END OF APPROXIMATION ROUTINE
01260REM
01270REM *********************************************************
01280PRINT
01290GOSUB 4700
01300GOSUB 9000
01310IF O1=0 THEN 1380
01320GOTO 1500
01330REM *************************************************************
01340REM
01350REM ROUTINE FOR RESPECIFYING PERCENTILES
01360REM
01370REM *************************************************************
01380PRINT L$
01390PRINT "RESPECIFY PERCENTILES."
01400PRINT
01410GOTO 440
01420REM
01430REM
01440REM END OF SPECIFYING PERCENTILES
01450REM
01460REM ***************************************************************
01470REM
01480REM DETERMINES WHICH DISTRIBUTION HE WANTS
01490REM
01500IF O1=1 THEN 1620
01510IF O1=2 THEN 1620
01520IF O1=3 THEN 1620
01530IF O1=4 THEN 1620
01540PRINT "REENTER. MUST BE 0 OR NUMBER OF A DISTRIBUTION."
01550GOSUB 9000
01560GOTO 1310
01570REM ********************************************************
01580REM
01590REM PRINT PERCENTILES
01600REM
01610REM
01620K5=O1
01630A=P(K5,1)
01640B=P(K5,2)
01650GOTO 3620
01660IF A>1.15 THEN 1780
01670M2=2
01680M7=M2/Q2
01690M8=(M2-1)/Q2+2
01700IF M7>M8 THEN 1720
01710M7=M8
01720IF M7-M2>1.15*K5 THEN 1750
01730M2=M2+1
01740GOTO 1680
01750A=M2
01760B=M7-M2
01770GOTO 1790
01780IF B<1.15 THEN 1670
01790IF A+B<8000 THEN 1820
01800A=8000*Q2
01810B=8000-A
01820P1=.1
01830P(K5,1)=A
01840P(K5,2)=B
01850GOSUB 4400
01860S(K5,1)=J2
01870P1=.25
01880GOSUB 4420
01890S(K5,2)=J2
01900P1=.5
01910GOSUB 4420
01920S(K5,3)=J2
01930P1=.75
01940GOSUB 4420
01950S(K5,4)=J2
01960P1=.9
01970GOSUB 4420
01980S(K5,5)=J2
01990IF S9=1 THEN 2090
02000IF K5=1 THEN 2030
02010PRINT USING 2020,K5,S(K5,1),S(K5,2),S(K5,3),S(K5,4),S(K5,5)
02020: ## ##.### ##.### ##.### ##.### ##.##
02030RETURN
02040REM
02050REM
02060PRINT L$
02070PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE BETA"
02080PRINT "DISTRIBUTION FITTED TO YOUR PERCENTILE SPECIFICATIONS."
02090PRINT
02100M=A+B
02110: HYPOTHETICAL SAMPLE SIZE (T)########.##
02120: STANDARD DEVIATION ##.##
02130: 10TH PERCENTILE ##.##
02140: 25TH PERCENTILE ##.##
02150: 50TH (MEDIAN) ##.##
02160: 75TH PERCENTILE ##.##
02170: 90TH PERCENTILE ##.##
02180: 50% HDR ##.## - ##.##
02190: 75% HDR ##.## - ##.##
02200: MODE ##.##
02210: 95% HDR ##.## - ##.##
02220: MEAN ##.##
02230: PARAMETER A ########.##
02240: PARAMETER B ########.##
02250J5=.5
02260J8=S(K5,4)/S(K5,2)
02270GOSUB 7000
02280PRINT USING 2110,A+B
02290IF V7=0 THEN 2310
02300PRINT USING 2200,(A-1)/(M-2)
02310PRINT USING 2130,S(K5,1)
02320PRINT USING 2140,S(K5,2)
02330PRINT USING 2150,S(K5,3)
02340PRINT USING 2160,S(K5,4)
02350PRINT USING 2170,S(K5,5)
02360IF V7=1 THEN 2400
02370IF J2<.99 THEN 2390
02380J2=.99
02390GOTO 2420
02400J1=H(1)
02410J2=H(2)
02420PRINT USING 2180,J1,J2
02430IF V7=1 THEN 2520
02440H(1)=J1
02450H(2)=J2
02460J5=.75
02470J8=J1/J2
02480GOSUB 7000
02490IF J2<.99 THEN 2510
02500J2=.99
02510GOTO 2540
02520J1=H(3)
02530J2=H(4)
02540PRINT USING 2190,J1,J2
02550IF V7=1 THEN 2640
02560H(3)=J1
02570H(4)=J2
02580J5=.95
02590J8=J2/J1
02600GOSUB 7000
02610IF J2<.99 THEN 2630
02620J2=.99
02630GOTO 2660
02640J1=H(5)
02650J2=H(6)
02660PRINT USING 2210,J1,J2
02670PRINT
02680IF V7=1 THEN 3670
02690H(5)=J1
02700H(6)=J2
02710IF S9=1 THEN 2850
02720PRINT "IF YOU DO NOT FEEL THAT THE HYPOTHETICAL SAMPLE SIZE (T) "
02730PRINT "REFLECTS YOUR PRIOR INFORMATION ABOUT THE PROPORTION YOU CAN"
02740PRINT "SPECIFY A DIFFERENT T. THIS WILL NOT AFFECT THE MEDIAN BUT"
02750PRINT "WILL CHANGE OTHER PERCENTILES AND THE HDRS. A LARGER T WILL"
02760PRINT "RESULT IN SHORTER INTERVALS, AND A SMALLER T IN LONGER ONES."
02770PRINT
02780REM ***********************************************************
02790REM
02800REM ROUTINE TO INSURE THAT THE PARAMETERS ARE OK AND TO
02810REM MINIMAL VALUE THAT M CAN ASSUME
02820REM
02830REM
02840IF V2=1 THEN 2970
02850M2=3
02860V2=1
02870M7=3*S(K5,3)*M2*(M2-2)+M2
02880M7=M7/(3*M2-4)
02890IF M7<1.15 THEN 2910
02900IF M2-M7>1.15 THEN 2930
02910M2=M2+1
02920GOTO 2870
02930M7=INT(M2+1)
02940IF M7<M THEN 2960
02950M7=M
02960:IF YOU WANT TO CHANGE T, TYPE THE NEW VALUE (MIN=#####.###).
02970PRINT USING 2960,M7
02980PRINT "IF YOU DO NOT WANT TO CHANGE T, TYPE '0'."
02990GOSUB 9000
03000IF O1=0 THEN 3100
03010IF O1>9999 THEN 3050
03020IF O1 >= M7 THEN 3070
03030PRINT "REENTER. INPUT MUST BE 0 OR AN ACCEPTABLE VALUE."
03040GOTO 2990
03050PRINT "REENTER. T MUST NOT BE GREATER THAN 9999."
03060GOTO 2990
03070M=O1
03080O1=S(K5,3)
03090GOTO 3310
03100PRINT
03110IF S8=0 THEN 3150
03120PRINT "IF YOU WANT TO CHANGE THE MEDIAN TYPE THE NEW VALUE ELSE '0'.";
03130GOSUB 9000
03140GOTO 3220
03150PRINT "YOU CAN CHANGE THE CENTERING OF THE DISTRIBUTION BY SPECIFYING"
03160S8=1
03170PRINT "A DIFFERENT MEDIAN. THIS WILL NOT AFFECT THE VALUE OF T."
03180PRINT
03190PRINT "IF YOU WANT TO CHANGE THE MEDIAN TYPE THE NEW VALUE."
03200PRINT "IF YOU DO NOT TYPE '0'."
03210GOSUB 9000
03220IF O1=0 THEN 4320
03230IF O1 >= .03 THEN 3260
03240PRINT "REENTER. MEDIAN MUST BE AT LEAST .03."
03250GOTO 3210
03260IF O1 <= .97 THEN 3290
03270PRINT "REENTER. MEDIAN MUST NOT BE GREATER THAN .97."
03280GOTO 3210
03290A=O1*M
03300V2=0
03310S9=1
03320L4=O1*M
03330L3=O1*(M-2)+1
03340IF L4>L3 THEN 3380
03350A=L3
03360L3=L4
03370L4=A
03380A=.5*(L3+L4)
03390P1=.5
03400B=M-A
03410GOSUB 4400
03420IF ABS(J2-O1)<.004 THEN 3500
03430IF J2>O1 THEN 3470
03440L3=A
03450A=.5*(L4+A)
03460GOTO 3380
03470L4=A
03480A=.5*(L3+A)
03490GOTO 3380
03500IF A<1.1 THEN 3520
03510IF B>1.1 THEN 3620
03520PRINT L$
03530PRINT "THE MEDIAN YOU SPECIFIED IS NOT CONSISTENT WITH THE T"
03540V2=1
03550PRINT "VALUE YOU ARE USING. EXTREME MEDIAN VALUES REQUIRE LARGE"
03560PRINT "T VALUES. YOU SHOULD EITHER MODERATE YOUR MEDIAN ESTIMATE"
03570PRINT "OR INCREASE THE T VALUE."
03580PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
03590GOSUB 9000
03600A=P(K5,1)
03610B=P(K5,2)
03620PRINT L$
03630PRINT "HERE ARE SOME CHARACTERISTICS OF THE BETA DISTRIBUTION YOU"
03640PRINT "ARE NOW CONSIDERING."
03650IF S9=0 THEN 2090
03660GOTO 1820
03670PRINT "THIS COMPLETES THE SPECIFICATION OF THE PRIOR FOR PI. IF YOU"
03680PRINT "ARE NOT GOING DIRECTLY TO THE POSTERIOR ANALYSES YOU SHOULD "
03690PRINT "RECORD THE VALUES OF A AND B."
03700PRINT
03710PRINT "WHEN YOU ARE READY TO CONTINUE, TYPE '1'."
03720GOSUB 9000
03730PRINT L$
03740PRINT "YOU MUST NOW SPECIFY THE IMPORTANCE OF THE INFORMATION YOU HAVE"
03750PRINT "JUST GIVEN RELATIVE TO THE DATA YOU PLAN TO OBSERVE (OR HAVE"
03760PRINT "OBSERVED)."
03770PRINT
03780PRINT "THIS IMPORTANCE IS ASSESSED BY COMPARING YOUR KNOWLEDGE TO THAT"
03790PRINT "OBTAINED FROM EXACT INFORMATION CONCERNING THE PROBABILITIES OF"
03800PRINT "SUCCESS FOR SOME NUMBER OF GROUPS. THESE NUMBERS WILL BE USED,"
03810PRINT "TOGETHER WITH THE ACTUAL NUMBER OF GROUPS YOU OBSERVE, TO"
03820PRINT "DETERMINE THE RELATIVE IMPORTANCE OF PRIOR AND SAMPLE INFOR-"
03830PRINT "MATION."
03840PRINT
03850PRINT "FIRST CONSIDER THE OVERALL MEAN PROBABILITY OF SUCCESS TAKEN"
03860PRINT "ACROSS ALL GROUPS YOU MIGHT CONCEIVABLY OBSERVE. HOW MANY "
03870PRINT "GROUPS (T1) IS YOUR KNOWLEDGE CONCERNING THIS MEAN WORTH? (IN"
03880PRINT "MOST CASES A SMALL VALUE - BETWEEN 1 AND 5 - IS A GOOD CHOICE.)"
03890PRINT "T1 = ";
03900GOSUB 9000
03910IFO1<1THEN3920
03911IFO1<100THEN3940
03920PRINT "REENTER. INPUT MUST BE BETWEEN 1 AND 100."
03930GOTO 3890
03940T1=O1
03950PRINT
03960PRINT "NEXT CONSIDER THE BETWEEN GROUPS VARIANCE FOR PROBABILITY OF"
03970PRINT "SUCCESS, AGAIN TAKEN ACROSS ALL GROUPS. HOW MANY GROUPS (T2) IS"
03980PRINT "YOUR KNOWLEDGE CONCERNING THIS VARIANCE WORTH? (IN MOST CASES"
03990PRINT "A SMALL VALUE - BETWEEN 5 AND 10 - IS A GOOD CHOICE.)"
04000PRINT "T2 = ";
04010GOSUB 9000
04020IF O1<5 THEN 4030
04022IF O1 <= 100 THEN 4050
04030PRINT "REENTER. INPUT MUST BE BETWEEN 5 AND 100."
04040GOTO 4000
04050T2=O1
04060PRINT
04070PRINT "IF YOU ARE NOT GOING DIRECTLY TO THE POSTERIOR ANALYSES, YOU"
04080PRINT "SHOULD RECORD THE VALUES OF T1 AND T2."
04090PRINT
04100PRINT "IF YOU WANT TO DO THE POSTERIOR ANALYSIS TYPE '1' ELSE '0'.";
04110GOSUB 9000
04120IF O1=1 THEN 4180
04130IF O1=0 THEN 4170
04140PRINT
04150PRINT "REENTER. INPUT MUST BE 0 OR 1."
04160GOTO 4110
04170CHAIN "RSTRT"
04180SCRATCH#2
04181 PRINT # 2,A,B,T1,T2
04190REM THE 4 VALUES BEING PASSED ARE THE PARAMETERS OF THE FITTED BETA
04200REM (A AND B) AND THE PRIOR SAMPLE SIZES FOR THE MEAN AND VARIANCE
04210REM OF GAMMA (T1 AND T2). THESE LAST 2 ARE ETA AND NU + 1 IN THE
04220REM NOTATION OF LEWIS,WANG & NOVICK (1975).
04230PRINT
04240PRINT "IF ALL YOUR GROUP SAMPLE SIZES ARE THE SAME TYPE '1' ELSE '0'."
04250GOSUB 9000
04260IF O1=0 THEN 4300
04270IF O1=1 THEN 4310
04280PRINT "REENTER. INPUT MUST BE 0 OR 1."
04290GOTO 4250
04300CHAIN "CMOD9"
04310CHAIN "CMODA"
04320V7=1
04330PRINT L$
04340PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE PRIOR DISTRIBUTION"
04350PRINT "FITTED TO YOUR PRIOR BELIEFS ABOUT PI."
04360PRINT
04370PRINT USING 2230,A
04380PRINT USING 2240,B
04390GOTO 2110
04400E1=0
04410A7=0
04420E2=1
04430A8=.999999
04440A6=P1
04450GOSUB 5220
04460S5=(A6-A7)/(A8-A7)
04470IF S5<.85 THEN 4500
04480S5=.85
04490GOTO 4520
04500IF S5>.15 THEN 4520
04510S5=.15
04520J2=E1+S5*(E2-E1)
04530IF J2 >= 1 THEN 4650
04540J1=0
04550GOSUB 5000
04560IF ABS(E1-E2)<.0001 THEN 4670
04570IF ABS(P-P1)<.002 THEN 4670
04580IF P>P1 THEN 4630
04590E1=J2
04600A7=P
04610GOTO 4460
04620A8=P
04630E2=J2
04640GOTO 4460
04650J2=1
04660GOTO 4670
04670IF J2<.99 THEN 4690
04680J2=.99
04690RETURN
04700PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE"
04710PRINT "WHICH MOST CLOSELY CORRESPONDS TO YOUR PRIOR BELIEFS."
04720PRINT "YOU CAN EITHER TENTATIVELY ACCEPT THIS DISTRIBUTION OR"
04730PRINT "RESPECIFY THE PERCENTILES."
04740PRINT
04750PRINT "TYPE THE NUMBER OF THE DISTRIBUTION YOU WANT OR"
04760PRINT "TYPE '0' IF YOU WANT TO RESPECIFY THE PERCENTILES."
04770RETURN
04780REM **************************************************************
04790REM
04800REM APPENDED GOSUBS FOLLOW
04810REM
04820REM *****************************************************************
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=10**37
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 ****************************************************
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
07124J1=(-2*J1+SQR(4*J1*J1+4*C*(J2+J1)*(J2-J1)))/(2*(J2-J1))
07126IF 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 **********************************************************
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
09999END