Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmodab.bas
There are 2 other files named cmodab.bas in the archive. Click here to see a list.
00020REM**********************************************************
00030REM CMODAB CMODAB CMODAB CMODAB CMODAB CMODAB
00040REM**********************************************************
00050REM
00060REM THIS MODULE WILL ASSIST YOU IN CARRYING OUT A PREPOSTERIOR
00070REM ANALYSIS FOR THE BETA BINOMIAL MODEL
00080REM
00090REM****************************************************
00100 FILES RFILE1,RFILE2,RFILE3
00140RESTORE#1
00141 INPUT# 1,I1,I2,I3
00150SCRATCH#1
00151 PRINT # 1,37,I2,I3
00160 REM
00170X=0
00180PRINT L$
00190 DIM B(500)
00200GOSUB 5160
00210PRINT L$
00220K5=0
00230P=0
00240P4=0
00250PRINT " BETA BINOMIAL PREPOSTERIOR ANALYSIS"
00260PRINT
00270PRINT "THIS MODULE WILL ASSIST YOU IN CARRYING OUT A PREPOSTERIOR"
00280PRINT "ANALYSIS USING YOUR PRIOR DISTRIBUTION AND AN ADVERSARY"
00290PRINT "PRIOR DISTRIBUTION."
00300PRINT
00310RESTORE#2
00311 INPUT# 2,A,B
00320IF A <> 0 THEN 470
00330PRINT "INPUT THE PARAMETERS OF YOUR PRIOR DISTIBUTION."
00340PRINT
00350PRINT "INPUT A.";
00360GOSUB 9000
00370IF O1 >= 1.1 THEN 400
00380PRINT "REENTER. PARAMETERS MUST BE AT LEAST 1.1."
00390GOTO 360
00400A=O1
00410PRINT "INPUT B.";
00420GOSUB 9000
00430IF O1 >= 1.1 THEN 460
00440PRINT "REENTER. PARAMETERS MUST BE AT LEAST 1.1."
00450GOTO 420
00460B=O1
00465PRINT L$
00470GOSUB 5220
00490PRINT "THERE ARE TWO STEPS TO THE ANALYSIS. THE PURPOSE OF THE"
00500PRINT "FIRST STEP IS TO GIVE YOU A ROUGH IDEA OF THE EFFECT OF"
00510PRINT "DIFFERENT SAMPLE SIZES ON THE EXPECTED MEANS OF YOUR"
00520PRINT "ADVERSARY POSTERIOR DISTRIBUTION. THESE ARE THE MEANS"
00530PRINT "YOU WOULD EXPECT ACCORDING TO YOUR PRIOR DISTRIBUTION."
00540PRINT
00550PRINT "ONCE YOU HAVE A ROUGH IDEA OF THE SAMPLE SIZE YOU WANT"
00560PRINT "YOU CAN PROCEED TO THE SECOND STEP AND LOOK IN MORE DETAIL"
00570PRINT "AT YOUR EXPECTED ADVERSARY POSTERIOR DISTRIBUTIONS. YOU"
00580PRINT "WILL BE ABLE TO GET THE PROBABILITY THAT PI IS GREATER "
00590PRINT "THAN CERTAIN VALUES TO BE SPECIFIED BY YOU."
00600PRINT
00610GOSUB 2530
00620PRINT "INPUT THE PARAMETERS OF THE ADVERSARY PRIOR DISTRIBUTION."
00630PRINT
00640PRINT "INPUT A.";
00650GOSUB 9000
00660A1=O1
00662IF O1 >= 1.1 THEN 670
00664PRINT "PARAMETERS MUST BE AT LEAST 1.1. PLEASE RESPECIFY."
00666GOTO 640
00670PRINT "INPUT B.";
00680GOSUB 9000
00690B1=O1
00692IF O1 >= 1.1 THEN 700
00694PRINT "PARAMETER MUST BE AT LEAST 1.1. RESPECIFY."
00696GOTO 670
00700PRINT
00710L1=.41
00720U1=.5
00730PRINT L$
00740PRINT "INPUT THE NUMBER OF DIFFERENT SAMPLE SIZES YOU WANT TO"
00750PRINT "CONSIDER. (EXIT=0 MAX=10)";
00760GOSUB 9000
00770IF O1=0 THEN 1450
00780IF O1 >= 1 THEN 810
00790PRINT "REENTER. NUMBER MUST BE AT LEAST 1 BUT NOT GREATER THAN 10."
00800GOTO 760
00810IF O1>10 THEN 790
00820N5=O1
00830GOSUB 850
00840GOTO 940
00850GOTO 930
00860FOR J2=L1 TO U1 STEP .01
00870GOSUB 1350
00880NEXT J2
00890L1=U1+.01
00900U1=L1+.1
00910IF U1 <= .99 THEN 930
00920U1=.99
00930RETURN
00940PRINT
00950PRINT "INPUT THE SAMPLE SIZES."
00960PRINT
00970FOR K9=1 TO N5
00980 REM
00990 PRINT "SAMPLE SIZE ";K9;
01000GOSUB 9000
01010IF O1 >= 1 THEN 1040
01020PRINT "REENTER. SAMPLE SIZE MUST BE AT LEAST 1."
01030GOTO 1000
01040V(K9)=O1
01050NEXT K9
01060M=A/(A+B)
01070T1=A1+B1
01080M1=A1/(A1+B1)
01090PRINT L$
01100PRINT "HERE ARE THE MEANS OF THE PRIOR DISTRIBUTIONS AND THE EXPECTED"
01110PRINT "MEANS OF THE ADVERSARY POSTERIOR DISTRIBUTION."
01120PRINT
01130PRINT "INVESTIGATOR PRIOR DISTRIBUTIONS ADVERSARY"
01140PRINT "------------------------------------------------------------"
01150:####.## PARAMETER A #####.##
01160PRINT USING 1150,A,A1
01170:####.## PARAMETER B #####.##
01180PRINT USING 1170,B,B1
01190: ##.## MEAN ##.##
01200PRINT USING 1190,M,M1
01210PRINT "============================================================"
01220PRINT " SAMPLE EXPECTED MEAN OF ADVERSARY"
01230PRINT " SIZE POSTERIOR DISTRIBUTION"
01240FOR K9=1 TO N5
01250N=V(K9)
01260: ### ##.##
01270PRINT USING 1260,N,(N*M+A1)/(N+T1)
01280NEXT K9
01290PRINT "--------------------------------------------------------"
01300GOSUB 850
01310PRINT "IF YOU WANT TO TRY MORE N VALUES TYPE '1' ELSE '0'.";
01320GOSUB 9000
01330IF O1=0 THEN 1450
01340GOTO 730
01350J1=0
01360IF P>.9999 THEN 1430
01370GOSUB 5000
01380IF P<.0001 THEN 1430
01390K5=K5+1
01400I(K5)=P-P4
01410P4=P
01420E(K5)=J2-.005
01430RETURN
01440REM::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
01450PRINT L$
01460PRINT "THIS IS THE BEGINNING OF THE SECOND STEP IN THE ANALYSIS."
01470PRINT
01480PRINT "THE MODULE WILL COMPUTE AND PRINT THE PROBABILITY THAT PI"
01490PRINT "IS LESS THAN PI' FOR THE EXPECTED ADVERSARY POSTERIOR"
01500PRINT "DISTRIBUTION. YOU ARE TO SPECIFY PI' AND THE SAMPLE SIZE."
01510PRINT
01520PRINT "THE MINIMUM SAMPLE SIZE IS 5 AND THE MAXIMUM IS 200."
01530PRINT
01540PRINT "INPUT THE SAMPLE SIZE YOU WANT TO CONSIDER. (NONE=0)";
01550GOSUB 9000
01560IF O1=0 THEN 2230
01570IF O1>200 THEN 1600
01580IF O1 >= 5 THEN 1620
01600PRINT "REENTER. MINIMUM=5. MAXIMUM=200"
01610GOTO 1550
01620PRINT
01630N5=1
01640V(1)=O1
01650N=O1
01660PRINT "INPUT THE NUMBER OF PI' VALUES YOU WANT TO SPECIFY (MAX=4).";
01670GOSUB 9000
01680IF O1 >= 1 THEN 1710
01690PRINT "REEENTER. NUMBER MUST BE AT LEAST 1 BUT NOT GREATER THAN 4."
01700GOTO 1670
01710IF O1>4 THEN 1690
01720N9=O1
01730PRINT
01740FOR K4=1 TO N9
01750 REM
01760 PRINT"INPUT PI' ";
01770GOSUB 9000
01780IF O1 >= .01 THEN 1810
01790PRINT "REENTER. PI' MUST BE AT LEAST .01 BUT NOT GREATER THAN .99."
01800GOTO 1770
01810IF O1>.99 THEN 1790
01820Q(K4)=O1
01830NEXT K4
01840T7=0
01850T8=0
01860A5=A
01870B5=B
01880N1=A+B
01890R1=A
01900PRINT L$
01910PRINT "HERE ARE THE PRIOR AND POSTERIOR PROBABILITIES THAT PI IS"
01920PRINT "LESS THAN PI' FOR A SAMPLE OF SIZE N."
01930PRINT
01940: N =####
01950PRINT USING 1940,V(1)
01960PRINT
01970PRINT " YOU ADVERSARY"
01980PRINT " PI' PRIOR/POSTERIOR PRIOR POSTERIOR"
01990PRINT "----------------------------------------------------------"
02000FOR Q7=1 TO N9
02010A=A5
02020B=B5
02030GOSUB 5220
02040J2=Q(Q7)
02050J1=0
02060GOSUB 5000
02070P2=P
02080REM
02090A=A1
02100B=B1
02110J1=0
02120GOSUB 5220
02130GOSUB 5000
02140P3=P
02150A6=A
02160S8=A+B
02170GOSUB 2410
02180: ##.## ##.## ##.## ##.##
02190PRINT USING 2180,J2,P2,P3,T7
02200NEXT Q7
02210PRINT "----------------------------------------------------------"
02220GOTO 2240
02230PRINT L$
02240PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
02250PRINT " 1. TRY A DIFFERENT SET OF N AND/OR PI' VALUES."
02260PRINT " 2. DO ANOTHER PREPOSTERIOR ANALYSIS WITH A DIFFERENT"
02270PRINT " ADVERSARY."
02280PRINT " 3. EXIT MODULE"
02290GOSUB 9000
02300A=A5
02310B=B5
02320IF O1=1 THEN 1510
02330IF O1=2 THEN 2402
02340IF O1=3 THEN 2350
02350PRINT
02360PRINT "IF YOU WANT TO DO THE POSTERIOR ANALYSIS TYPE '1' ELSE '0'.";
02370GOSUB 9000
02380IF O1=0 THEN 2400
02390CHAIN "CMOD3"
02400CHAIN "RSTRT"
02402PRINT L$
02403GOTO 620
02410GOSUB 2580
02420K7=1
02430T7=0
02440FOR A=A6 TO A6+N
02450B=S8+N-A
02460GOSUB 5220
02470J1=0
02480GOSUB 5000
02490T7=T7+P*B(K7)
02500K7=K7+1
02510NEXT A
02520RETURN
02530PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
02540GOSUB 9000
02550PRINT L$
02560RETURN
02570REM************************** BETA-BINOMIAL *******************
02580 DIM N(4),D(5)
02590R=INT(R1*N/N1)
02600N(1)=R+R1-1
02610N(2)=N+N1-R-R1-1
02620N(3)=N
02630N(4)=N1-1
02640D(1)=R
02650D(2)=R1-1
02660D(3)=N-R
02670D(4)=N1-R1-1
02680D(5)=N+N1-1
02690F0=0
02700FOR K5=1 TO 4
02710G9=N(K5)
02720G9=G9+1
02730GOSUB 5850
02740F0=F0+G0
02750G9=D(K5)
02760G9=G9+1
02770GOSUB 5850
02780F0=F0-G0
02790NEXT K5
02800G9=D(5)
02810G9=G9+1
02820GOSUB 5850
02830F0=F0-G0
02840N8=N
02850F8=EXP(F0)
02860B(R+1)=F8
02870T8=F8
02880F6=F8
02890FOR R8=R+1 TO N
02900GOSUB 3020
02910F8=F8*F7
02920T8=T8+F8
02930B(R8+1)=F8
02940NEXT R8
02950F8=F6
02960FOR R8=R TO 1 STEP -1
02970GOSUB 3020
02980F8=F8/F7
02990B(R8)=F8
03000NEXT R8
03010RETURN
03020F7=(R1-1+R8)*(N8+1-R8)/(R8*(N8+N1-R1-R8))
03030RETURN
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=-1
05448GOTO 5460
05451G6=1
05454GOTO 5460
05457G6=0
05460RETURN
05465REM END BETA CDF ROUTINE
05470REM***********************************************************
05600REM************************************************************
05606REM BETA PERCENTILE ROUTINE
05609REM INPUTS P1
05612REM OUTPUTS J2
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
05753REM END OF PERCENTILE ROUTINE
05756REM
05759REM ***************************************************
05762REM*********************************************************
05768REM INVERSE BETA FUNCTION
05771REM INPUTS P2
05774REM OUTPUTS J2
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 ****************************************************
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 ****************************************************
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)
08021IF X*X/2<80 THEN 8025
08022D=0
08023GOTO 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
09999END