Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmodac.bas
There are 2 other files named cmodac.bas in the archive. Click here to see a list.
00020REM****************************************************************
00030REM CMODAC CMODAC CMODAC CMODAC CMODAC CMODAC
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
00100REM
00110REM*****************************************************************
00120REM
00130REM NEEDED FOR NDTR
00140X=0
00150V7=1
00160 FILES RFILE1,RFILE2,RFILE3
00200RESTORE#1
00201 INPUT# 1,I1,I2,I3
00210SCRATCH#1
00211 PRINT # 1,38,I2,I3
00220RESTORE#2
00221 INPUT# 2,A,B
00230PRINT L$
00240PRINT " POSTERIOR ANALYSIS BETA PASCAL MODEL"
00250REM*******************************************************************
00260REM
00270REM IF A=0 THEN THE PARAMETERS OF PRIOR NOT PASSED
00280REM
00290REM*********************************************************************
00300IF A <> 0 THEN 470
00310PRINT
00320PRINT "THIS MODULE COMBINES THE PRIOR AND SAMPLE INFORMATION AND"
00330PRINT "PROVIDES THE POSTERIOR DISTRIBUTION ON PI."
00340PRINT
00350PRINT "INPUT PARAMETER A OF YOUR PRIOR DISTRIBUTION.";
00360GOSUB 9000
00370IF O1<1.1 THEN 450
00380A=O1
00390PRINT
00400PRINT "INPUT PARAMETER B";
00410GOSUB 9000
00420IF O1<1.1 THEN 450
00430B=O1
00440GOTO 470
00450PRINT "PARAMETERS MUST BE 1.1 OR GREATER."
00460GOTO 350
00470GOSUB 5160
00480PRINT
00490PRINT "INPUT NUMBER OF SAMPLE OBSERVATIONS.";
00500GOSUB 9000
00510IF ABS(INT(O1)-O1)<.0001 THEN 540
00520PRINT "REENTER. MUST BE INTERGER 1 OR GREATER."
00530GOTO 480
00540IF O1<1 THEN 510
00550PRINT
00560N=O1
00570PRINT "INPUT NUMBER OF SUCCESSES.";
00580GOSUB 9000
00590IF ABS(INT(O1)-O1)<.0001 THEN 620
00600PRINT "REENTER. MUST INTERGER 0 OR GREATER."
00610GOTO 580
00620IF O1 <= N THEN 650
00630PRINT "REENTER. MUST NOT BE GREATER THAN NUMBER OF OBSERVATIONS."
00640GOTO 550
00650IF O1 >= 0 THEN 680
00660PRINT "REENTER. MUST BE 0 OR GREATER."
00670GOTO 550
00680A1=A
00690B1=B
00700M1=A+B
00710PRINT
00720PRINT "SOME OF THE CHARACTERISTICS OF THE POSTERIOR DISTRIBUTION ARE"
00730PRINT "BEING COMPUTED."
00740S1=(A*B)/(M1*M1*(M1+1))
00750S0=O1
00760GOTO 810
00770A=A1+S0
00780B=B1+N-S0
00790M=A+B
00800S2=(A*B)/(M*M*(M+1))
00810K5=1
00820P1=.1
00830GOSUB 1750
00840S(K5,1)=J2
00850P1=.25
00860GOSUB 1760
00870S(K5,2)=J2
00880P1=.5
00890GOSUB 1760
00900S(K5,3)=J2
00910P1=.75
00920GOSUB 1760
00930S(K5,4)=J2
00940P1=.9
00950GOSUB 1760
00960S(K5,5)=J2
00970IF V7=0 THEN 1040
00980P5=S(K5,1)
00990P7=S(K5,2)
01000P8=S(K5,3)
01010P6=S(K5,4)
01020V5=S(K5,5)
01030GOTO 1260
01040PRINT L$
01050PRINT " SUMMARY OF BETA PASCAL ANALYSIS"
01060PRINT
01090PRINT" PRIOR BETA DISTRIBUTIONS POSTERIOR"
01100PRINT"---------------------------------------------------------"
01110M=A+B
01120:########.## HYPOTHETICAL SAMPLE SIZE#######.##
01130: ##.#### STANDARD DEVIATION ##.####
01140: ##.## 10TH PERCENTILE ##.##
01150: ##.## 25TH PERCENTILE ##.##
01160: ##.## 90TH PERCENTILE ##.##
01170: ##.## 50TH PERCENTILE ##.##
01180: ##.## 75TH PERCENTILE ##.##
01190: ##.## -##.## 50% HDR ##.## -##.##
01200: ##.## -##.## 75% HDR ##.## -##.##
01210: ##.## MODE ##.##
01220: ##.## -##.## 95% HDR ##.## -##.##
01230: ##.## MEAN ##.##
01240:########.## PARAMETER A #######.##
01250:########.## PARAMETER B #######.##
01260J5=.5
01270J8=1.1
01280GOSUB 7000
01290IF V7=0 THEN 1330
01300H1=J1
01310H2=J2
01320GOTO 1440
01330PRINT USING 1240,A1,A
01340PRINT USING 1250,B1,B
01350PRINT USING 1130,SQR(S1),SQR(S2)
01360PRINT USING 1140,P5,S(K5,1)
01370PRINT USING 1150,P7,S(K5,2)
01380PRINT USING 1170,P8,S(K5,3)
01390PRINT USING 1180,P6,S(K5,4)
01400PRINT USING 1160,V5,S(K5,5)
01410PRINT USING 1230,A1/M1,A/M
01420PRINT USING 1210,(A1-1)/(M1-2),(A-1)/(M-1)
01430PRINT USING 1190,H1,H2,J1,J2
01440J5=.75
01450J8=J2/J1
01460GOSUB 7000
01470IF V7=0 THEN 1510
01480H3=J1
01490H4=J2
01500GOTO 1520
01510PRINT USING 1200,H3,H4,J1,J2
01520J5=.95
01530J8=J2/J1
01540GOSUB 7000
01550IF V7=0 THEN 1600
01560H5=J1
01562IF H5>.01 THEN 1570
01564H5=.01
01570H6=J2
01572IF H6<.99 THEN 1580
01574H6=.99
01580V7=0
01590GOTO 770
01600REM
01610 IF J2<.99 THEN 1620
01615 J2=.99
01620PRINT USING 1220,H5,H6,J1,J2
01630PRINT"---------------------------------------------------------"
01640REM
01650PRINT "IF YOU WANT TO EVALUATE POSTERIOR DISTRIBUTION TYPE '1'."
01655PRINT "IF YOU WANT TO EVALUATE PREDICTIVE DISTRIBUTION TYPE '2'."
01660PRINT "IF NOT, TYPE '0'."
01670GOSUB 9000
01680IF O1=0 THEN 1740
01690IF O1=1 THEN 1720
01695IF O1=2 THEN 1742
01700PRINT "REENTER. MUST BE 0,1 OR 2."
01710GOTO 1670
01720SCRATCH#3
01721 PRINT # 3,A,B
01730CHAIN "CMODB"
01740CHAIN "RSTRT"
01742SCRATCH#3
01743 PRINT # 3,A,B
01745CHAIN "CMODY"
01750E1=0
01760E2=1
01770GOSUB 5220
01780J2=(E1+E2)/2
01790J1=0
01800GOSUB 5000
01810IF ABS(P-P1)<.009 THEN 1870
01820IF P>P1 THEN 1850
01830E1=J2
01840GOTO 1780
01850E2=J2
01860GOTO 1780
01870RETURN
05000REM ***************************************************
05005REM BETA CDF ROUTINE
05010REM INPUT A B J2
05015REM OUTPUT P
05020REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220
05025 DIM W(16),O(16)
05030GOSUB 5340
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A>1 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=F0+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350IF D2<-80 THEN 5370
05355IF D2>85 THEN 5380
05360D2=EXP(D2)
05365RETURN
05370D2=1.E-37
05375RETURN
05380D2=1.E+37
05385RETURN
05390REM
05395REM END OF BETA CDF ROUTINE
05400REM *******************************************************
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 ****************************************************
07000REM**************************************************************
07002REM BETA HDR ROUTINE
07004REM INPUT A B J5
07006REM OUTPUT J1 J2
07008J8=2
07010GOSUB 5220
07012REM
07014J7=1
07016GOSUB 7082
07018IF J8<-99 THEN 7028
07020P9=0
07022GOSUB 5000
07024P0=P
07026IF ABS(P0-J5)>.001 THEN 7030
07028RETURN
07030IF P0>J5 THEN 7038
07031J7=J8
07032J8=J8*J8
07034P9=P0
07036GOTO 7016
07038J9=J7
07040J0=J8
07042J8=(J5-P9)/(P0-P9)
07044IF J8>.15 THEN 7050
07046J8=.15
07048GOTO 7054
07050IF J8<.85 THEN 7054
07052J8=.85
07054J8=J8*(J0-J9)+J9
07056GOSUB 7082
07058IF J8<-99 THEN 7066
07060GOSUB 5000
07062J3=P
07064IF ABS(J3-J5)>.001 THEN 7068
07066RETURN
07068IF J3>J5 THEN 7076
07070J9=J8
07072P9=J3
07074GOTO 7042
07076J0=J8
07078P0=J3
07080GOTO 7042
07082J=(A-1)/(B-1)
07084J1=(J8**J-1)/(J8**(J+1)-1)
07086J2=J8*J1
07088IF J2<1 THEN 7094
07090J8=J8*.95
07092GOTO 7084
07094RETURN
07096REM END OF BETA HDR ROUTINE
07098REM*************************************************
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)
08022IF X*X<80 THEN 8025
08023D=0
08024GOTO 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