Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod3.bas
There are 2 other files named cmod3.bas in the archive. Click here to see a list.
00020REM****************************************************************
00030REM CMOD3 CMOD3 CMOD3 CMOD3 CMOD3 CMOD3
00040REM****************************************************************
00050REM
00060REM THIS MODULES COMBINES THE PRIOR DISTRIBUTION ON PI WITH THE
00070REM SAMPLE DATA IN ORDER TO GET THE POSTERIOR DISTRIBUTION ON PI.
00080REM
00090REM
00100REM*****************************************************************
00110REM
00120REM RUN THROUGH 1TOC1 R0=0 / R1=1 / R2=2 / R3=3 / R4=4 / R5=5
00130READ R0,R1,R2,R3,R4,R5
00140DATA 0,1,2,3,4,5
00150GOSUB 5160
00160 DIM S(1,5)
00170X=R0
00180V6=R0
00190Z0=R0
00200V7=R1
00210FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00250RESTORE#1
00251 INPUT# 1,I1,I2,I3
00260SCRATCH#1
00261 PRINT # 1,3,I2,I3
00270RESTORE#2
00271 INPUT# 2,A,B
00280PRINT L$
00290PRINT " POSTERIOR ANALYSIS BETA-BINOMIAL MODEL"
00300REM*******************************************************************
00310REM
00320REM IF A=0 THEN THE PARAMETERS OF PRIOR NOT PASSED
00330REM
00340REM*********************************************************************
00350IF A <> R0 THEN 780
00360PRINT
00370PRINT "THIS MODULE COMBINES THE PRIOR AND SAMPLE INFORMATION AND"
00380PRINT "PROVIDES THE POSTERIOR DISTRIBUTION ON PI."
00390PRINT
00400PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00410PRINT " 1. POSTERIOR ANALYSIS WITH INFORMATIVE PRIOR"
00420PRINT " 2. POSTERIOR ANALYSIS WITH NON-INFORMATIVE PRIOR"
00430PRINT " 3. ASSESSMENT OF YOUR PRIOR DISTRIBUTION ON PI"
00440PRINT " 4. EXIT"
00450GOSUB 9000
00460IF O1=R4 THEN 2160
00470IF O1=R3 THEN 2180
00480IF O1=R2 THEN 520
00490IF O1=R1 THEN 650
00500PRINT "REENTER. INPUT MUST BE NUMBER OF AVAILABLE OPTION."
00510GOTO 450
00520A1=R1
00530B1=R1
00540S1=R0
00550P5=R0
00560P7=R0
00570P8=R0
00580P6=R0
00590V5=R0
00600M1=R2
00610Z0=R1
00620V6=R1
00630PRINT L$
00640GOTO 790
00650PRINT L$
00660PRINT "INPUT PARAMETER A OF YOUR PRIOR DISTRIBUTION.";
00670GOSUB 9000
00680IF O1<1.1 THEN 760
00690A=O1
00700PRINT
00710PRINT "INPUT PARAMETER B";
00720GOSUB 9000
00730IF O1<1.1 THEN 760
00740B=O1
00750GOTO 780
00760PRINT "PARAMETERS MUST BE 1.1 OR GREATER."
00770GOTO 660
00780PRINT
00790PRINT "INPUT NUMBER OF SAMPLE OBSERVATIONS.";
00800GOSUB 9000
00810IF ABS(INT(O1)-O1)<.0001 THEN 840
00820PRINT "REENTER. SAMPLE SIZE MUST BE INTEGER VALUE."
00830GOTO 780
00840IF O1<R1 THEN 820
00850PRINT
00860N=O1
00870PRINT "INPUT NUMBER OF SUCCESSES.";
00880GOSUB 9000
00890IF ABS(INT(O1)-O1)<.0001 THEN 920
00900PRINT "REENTER. NUMBER OF SUCCESSES MUST BE INTEGER VALUE."
00910GOTO 880
00920IF O1 <= N THEN 950
00930PRINT "REENTER. MUST NOT BE GREATER THAN NUMBER OF OBSERVATIONS."
00940GOTO 880
00950IF O1 >= R0 THEN 980
00960PRINT "REENTER. MUST BE 0 OR GREATER."
00970GOTO 880
00980S0=O1
00990IF V6=R1 THEN 1090
01000A1=A
01010B1=B
01020M1=A+B
01030PRINT
01040PRINT "SOME OF THE CHARACTERISTICS OF THE POSTERIOR DISTRIBUTION ARE"
01050PRINT "BEING COMPUTED."
01060S1=(A*B)/(M1*M1*(M1+R1))
01070S0=O1
01080GOTO 1140
01090A=A1+S0
01100B=B1+N-S0
01110M=A+B
01120S2=(A*B)/(M*M*(M+R1))
01130IF V6=R1 THEN 2460
01140K5=R1
01150P1=.1
01160GOSUB 2190
01170S(K5,R1)=J2
01180P1=.25
01190GOSUB 2200
01200S(K5,R2)=J2
01210P1=.5
01220GOSUB 2200
01230S(K5,R3)=J2
01240P1=.75
01250GOSUB 2200
01260S(K5,R4)=J2
01270P1=.9
01280GOSUB 2200
01290S(K5,R5)=J2
01300IF V6=R1 THEN 1460
01310IF V7=R0 THEN 1380
01320P5=S(K5,R1)
01330P7=S(K5,R2)
01340P8=S(K5,R3)
01350P6=S(K5,R4)
01360V5=S(K5,R5)
01370GOTO 1610
01380PRINT L$
01390PRINT " SUMMARY OF BETA-BINOMIAL ANALYSIS"
01400PRINT
01410REM
01420PRINT" PRIOR BETA DISTRIBUTIONS POSTERIOR"
01430PRINT " ---------------------------------------------------------"
01440IF Z0=R0 THEN 1460
01450PRINT "NON-INFORMATIVE PRIOR"
01460M=A+B
01470:########.## HYPOTHETICAL SAMPLE SIZE#######.##
01480: ##.#### STANDARD DEVIATION ##.####
01490: ##.## 10TH PERCENTILE ##.##
01500: ##.## 25TH PERCENTILE ##.##
01510: ##.## 90TH PERCENTILE ##.##
01520: ##.## 50TH PERCENTILE ##.##
01530: ##.## 75TH PERCENTILE ##.##
01540: ##.## -##.## 50% HDR ##.## -##.##
01550: ##.## -##.## 75% HDR ##.## -##.##
01560: ##.## MODE ##.##
01570: ##.## -##.## 95% HDR ##.## -##.##
01580: ##.## MEAN ##.##
01590:########.## PARAMETER A #######.##
01600:########.## PARAMETER B #######.##
01610J5=.5
01620J8=1.1
01630GOSUB 7000
01640IF V7=R0 THEN 1680
01650H1=J1
01660H2=J2
01670GOTO 1790
01680 PRINT USING 1590,A1,A
01690 PRINT USING 1600,B1,B
01700 PRINT USING 1480,SQR(S1),SQR(S2)
01710 PRINT USING 1490,P5,S(K5,R1)
01720 PRINT USING 1500,P7,S(K5,R2)
01730 PRINT USING 1520,P8,S(K5,R3)
01740 PRINT USING 1530,P6,S(K5,R4)
01750 PRINT USING 1510,V5,S(K5,R5)
01760 PRINT USING 1580,A1/M1,A/M
01770 PRINT USING 1560,(A1-R1)/(M1-R2),(A-R1)/(M-R2)
01780 PRINT USING 1540,H1,H2,J1,J2
01790J5=.75
01800J8=J2/J1
01810GOSUB 7000
01820IF V7=R0 THEN 1860
01830H3=J1
01840H4=J2
01850GOTO 1870
01860 PRINT USING 1550,H3,H4,J1,J2
01870J5=.95
01880J8=J2/J1
01890GOSUB 7000
01900IF V7=R0 THEN 2000
01910H5=J1
01920IF H5>.01 THEN 1940
01930H5=.01
01940H6=J2
01950IF H6<.99 THEN 1970
01960H6=.99
01970IF V6=R1 THEN 2650
01980V7=R0
01990GOTO 1090
02000REM
02010IF J2<.99 THEN 2030
02020J2=.99
02030 PRINT USING 1570,H5,H6,J1,J2
02040PRINT " ----------------------------------------------------------"
02050PRINT "IF YOU WANT TO EVALUATE THE POSTERIOR TYPE '1',"
02060PRINT "IF YOU WANT TO EVALUATE THE PREDICTIVE DISTRIBUTION TYPE '2',"
02070PRINT "OTHERWISE TYPE '0'.";
02080SCRATCH#3
02081 PRINT # 3,A,B
02090GOSUB 9000
02100IF O1=R0 THEN 2160
02110IF O1=R1 THEN 2150
02120IF O1=R2 THEN 2170
02130PRINT "REENTER. MUST BE 0,1 OR 2."
02140GOTO 2090
02150CHAIN "CMODB"
02160CHAIN "RSTRT"
02170CHAIN "CMODX"
02180CHAIN "CMOD2"
02190E1=R0
02200A7=R0
02210GOSUB 5220
02220E2=R1
02230A8=.999999
02240A6=P1
02250S5=(A6-A7)/(A8-A7)
02260IF S5<.85 THEN 2280
02270S5=.85
02280IF S5>.15 THEN 2300
02290S5=.15
02300J2=E1+S5*(E2-E1)
02310J1=R0
02320GOSUB 2810
02330IF ABS(E1-E2)<.0001 THEN 2430
02340IF ABS(P-P1)<.005 THEN 2430
02350IF P>P1 THEN 2390
02360E1=J2
02370A7=P
02380GOTO 2250
02390A8=P
02400E2=J2
02410GOTO 2250
02420GOTO 2430
02430IF J2<.99 THEN 2450
02440J2=.99
02450RETURN
02460PRINT L$
02470PRINT " SUMMARY OF BETA-BINOMIAL ANALYSIS"
02480PRINT
02490PRINT " NON-INFORMATIVE PRIOR DISTRIBUTION A=1 B=1"
02500PRINT
02510PRINT " BETA POSTERIOR DISTRIBUTION"
02520PRINT "--------------------------------------------------------------"
02530: PARAMETER A ######.##
02540 PRINT USING 2530,A
02550: PARAMETER B ######.##
02560 PRINT USING 2550,B
02570: MODE ######.##
02580 PRINT USING 2570,(A-R1)/(A+B-R2)
02590: MEAN ######.##
02600 PRINT USING 2590,A/(A+B)
02610: STANDARD DEVIATION ######.####
02620 PRINT USING 2610,SQR(S2)
02630GOTO 1140
02640: 10TH PERCENTILE ######.##
02650 PRINT USING 2640,S(K5,R1)
02660: 25TH PERCENTILE ######.##
02670 PRINT USING 2660,S(K5,R2)
02680: 50TH PERCENTILE ######.##
02690 PRINT USING 2680,S(K5,R3)
02700: 75TH PERCENTILE ######.##
02710 PRINT USING 2700,S(K5,R4)
02720: 90TH PERCENTILE ######.##
02730 PRINT USING 2720,S(K5,R5)
02740: 50% HDR ##.## -##.##
02750 PRINT USING 2740,H1,H2
02760: 75% HDR ##.## -##.##
02770 PRINT USING 2760,H3,H4
02780: 95% HDR ##.## -##.##
02790 PRINT USING 2780,H5,H6
02800GOTO 2040
02810REM ***************************************************
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 <> R0 THEN 5035
05031IF A<R5 THEN 5035
05032IF B >= R5 THEN 5400
05035IF A+B>85 THEN 5280
05040P=R0
05045C6=R0
05050IF A <= B THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=R1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=R1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=R0 THEN 5115
05105IF D9=R1 THEN 5115
05107D9=LOG(D9)*(A-R1)+LOG(R1-D9)*(B-R1)
05108IF D9<-80 THEN 5115
05110P=P+W(I1)*EXP(D9)
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=R0 THEN 5155
05135A=C6
05140B=C7
05145P=R1-P
05150J2=R1-J2
05155RETURN
05160FOR I1=R1 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)**(R1/R3)
05285W2=(A*(R1-J2))**(R1/R3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)**(R1/R3)
05305W2=(A*(R1-J1))**(R1/R3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=R3*(W1*(R1-R1/9/B)-W2*(R1-R1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05340REM 2/16/76 CHANGED TO ALL LOG
05345D2=LOG(F0)+(A-R1)*LOG(J2)+(B-R1)*LOG(R1-J2)
05350RETURN
05390REM
05400REM *******************************************************
05403REM PEIZER PRATT APROXIMATION
05406D0=.02*(J2/B-(R1-J2)/A+(J2-.5)/(A+B))
05409D0=B-R1/R3-(A+B-R2/R3)*(R1-J2)+D0
05412Y=(B-.5)/((A+B-R1)*(R1-J2))
05415GOSUB 5436
05418D9=G6
05421Y=(A-.5)/((A+B-R1)*J2)
05424GOSUB 5436
05427Y3=D0*SQR((R1+J2*D9+(R1-J2)*G6)/((A+B-R5/6)*(R1-J2)*J2))
05430GOSUB 8000
05433RETURN
05436IF Y<.00001 THEN 5451
05439IF ABS(Y-R1)<.00001 THEN 5457
05442G6=(R1-Y*Y+R2*Y*LOG(Y))/(R1-Y)**R2
05445IF G6>-R1 THEN 5460
05446G6=-R1
05448GOTO 5460
05451G6=R1
05454GOTO 5460
05457G6=R0
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=R0
05878RETURN
05881IF G9<1.E+10 THEN 5890
05884G0=G9*(LOG(G9)-R1)
05887RETURN
05890G6=R1
05893IF 18<G5 THEN 5905
05896G6=G6*G5
05899G5=G5+R1
05902GOTO 5893
05905R8=R1/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+R1/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=R0
07010U0=R0
07012IF A<B THEN 7024
07014U0=R1
07016J7=A
07018A=B
07020J=R0
07022B=J7
07024J6=SQR((A*B)/(A+B)/(A+B)/(A+B+R1))
07026J=(A-R1)/(B-R1)
07028J0=(A-R1)/(A+B-R2)
07030GOSUB 5220
07032P9=R0
07034J9=J0
07036J6=1.6*J6
07038J8=A/(A+B)+J5*J6
07040IF J8<R1 THEN 7044
07042J8=.9999
07044IF A/(A+B)-J5*J6>R0 THEN 7050
07046J8=J8/1.E-08
07048GOTO 7052
07050J8=J8/(A/(A+B)-J5*J6)
07052GOSUB 7168
07054GOSUB 2810
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 2810
07074IF P<J5 THEN 7068
07076J6=J1
07078J7=J2
07080P0=P
07082GOTO 7110
07084J6=J1
07086J7=J2
07088P0=P
07090J8=(R1+J8)/R2
07092GOSUB 7168
07094GOSUB 2810
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+(R1-J1)*J9)/(J1*J6+(R1-J1)*J0)
07138GOSUB 7168
07140GOSUB 2810
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-R1)/(J8**(J+R1)-R1)
07170J2=J1*J8
07172IF J2<R1 THEN 7178
07174J8=J8*.95
07176GOTO 7172
07178RETURN
07180IF U0=R0 THEN 7194
07182J7=A
07184A=B
07186B=J7
07188J7=J2
07190J2=R1-J1
07192J1=R1-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=R1/(R1+.231642*Y4)
08025D=.398942*EXP(-X*X/R2)
08030C1=1.33027
08035C2=1.82126
08040C3=1.78148
08045C4=.356564
08050C5=.319382
08055P=R1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
08060IF X >= R0 THEN 8070
08065P=R1-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