10 X$=SYS(CHR$(7)):R$=LEFT(X$,7):L$=RIGHT(X$,8) 20 REM**************************************************************** 30 REM CMOD3 CMOD3 CMOD3 CMOD3 CMOD3 CMOD3 40 REM**************************************************************** 50 REM 60 REM THIS MODULES COMBINES THE PRIOR DISTRIBUTION ON PI WITH THE 70 REM SAMPLE DATA IN ORDER TO GET THE POSTERIOR DISTRIBUTION ON PI. 80 REM 90 REM IT IS ACCESSED BY: COMPONENT 1 100 REM 110 REM***************************************************************** 120 REM 130 REM NEEDED FOR NDTR 135 DIM S(5,5) 140 X=0 150 V7=1 160 DIM #1,A1(111),A2(111),A3(111) 170 OPEN R$ AS FILE 1 210 A1(0)=3 220 A=A2(0):B=A2(1) 230 PRINT L$ 240 PRINT " POSTERIOR ANALYSIS BETA-BINOMIAL MODEL" 250 REM******************************************************************* 260 REM 270 REM IF A=0 THEN THE PARAMETERS OF PRIOR NOT PASSED 280 REM 290 REM********************************************************************* 300 IF A <> 0 THEN 470 310 PRINT 320 PRINT "THIS MODULE COMBINES THE PRIOR AND SAMPLE INFORMATION AND" 330 PRINT "PROVIDES THE POSTERIOR DISTRIBUTION ON PI." 340 PRINT 350 PRINT "INPUT PARAMETER A OF YOUR PRIOR DISTRIBUTION."; 360 GOSUB 9000 370 IF O1<1.1 THEN 450 380 A=O1 390 PRINT 400 PRINT "INPUT PARAMETER B"; 410 GOSUB 9000 420 IF O1<1.1 THEN 450 430 B=O1 440 GOTO 470 450 PRINT "PARAMETERS MUST BE 1.1 OR GREATER." 460 GOTO 350 470 GOSUB 5160 480 PRINT 490 PRINT "INPUT NUMBER OF SAMPLE OBSERVATIONS."; 500 GOSUB 9000 510 IF ABS(INT(O1)-O1)<.0001 THEN 540 520 PRINT "REENTER. SAMPLE SIZE MUST BE INTEGER VALUE." 530 GOTO 480 540 IF O1<1 THEN 520 550 PRINT 560 N=O1 570 PRINT "INPUT NUMBER OF SUCCESSES."; 580 GOSUB 9000 590 IF ABS(INT(O1)-O1)<.0001 THEN 620 600 PRINT "REENTER. NUMBER OF SUCCESSES MUST BE INTEGER VALUE." 610 GOTO 580 620 IF O1 <= N THEN 650 630 PRINT "REENTER. MUST NOT BE GREATER THAN NUMBER OF OBSERVATIONS." 640 GOTO 580 650 IF O1 >= 0 THEN 680 660 PRINT "REENTER. MUST BE 0 OR GREATER." 670 GOTO 580 680 A1=A 690 B1=B 700 M1=A+B 710 PRINT 720 PRINT "SOME OF THE CHARACTERISTICS OF THE POSTERIOR DISTRIBUTION ARE" 730 PRINT "BEING COMPUTED." 740 S1=(A*B)/(M1*M1*(M1+1)) 750 S0=O1 760 GOTO 810 770 A=A1+S0 780 B=B1+N-S0 790 M=A+B 800 S2=(A*B)/(M*M*(M+1)) 810 K5=1 820 P1=.1 830 GOSUB 1750 840 S(K5,1)=J2 850 P1=.25 860 GOSUB 1760 870 S(K5,2)=J2 880 P1=.5 890 GOSUB 1760 900 S(K5,3)=J2 910 P1=.75 920 GOSUB 1760 930 S(K5,4)=J2 940 P1=.9 950 GOSUB 1760 960 S(K5,5)=J2 970 IF V7=0 THEN 1040 980 P5=S(K5,1) 990 P7=S(K5,2) 1000 P8=S(K5,3) 1010 P6=S(K5,4) 1020 V5=S(K5,5) 1030 GOTO 1260 1040 PRINT L$ 1050 PRINT " SUMMARY OF BETA-BINOMIAL ANALYSIS" 1060 PRINT 1080PRINT" PRIOR BETA DISTRIBUTIONS POSTERIOR" 1100PRINT" --------------------------------------------------------" 1110 M=A+B 1260 J5=.5 1270 J8=1.1 1280 GOSUB 7000 1290 IF V7=0 THEN 1330 1300 H1=J1 1310 H2=J2 1320 GOTO 1440 1330A$="########.## PARAMETER A #######.##" 1331 PRINT USING A$ ,A1,A 1340A$="########.## PARAMETER B #######.##" 1341 PRINT USING A$ ,B1,B 1350A$=" .#### STANDARD DEVIATION .####" 1351 PRINT USING A$ ,S1,S2 1360A$=" .## 10TH PERCENTILE .##" 1361 PRINT USING A$ ,P5,S(K5,1) 1370A$=" .## 25TH PERCENTILE .##" 1371 PRINT USING A$ ,P7,S(K5,2) 1380A$=" .## 50TH PERCENTILE .##" 1381 PRINT USING A$ ,P8,S(K5,3) 1390A$=" .## 75TH PERCENTILE .##" 1391 PRINT USING A$ ,P6,S(K5,4) 1400A$=" .## 90TH PERCENTILE .##" 1401 PRINT USING A$ ,V5,S(K5,5) 1410A$=" .## MEAN .##" 1411 PRINT USING A$ ,A1/M1,A/M 1420A$=" .## MODE .##" 1421 PRINT USING A$ ,(A1-1)/(M1-2),(A-1)/(M-1) 1430A$=" .## - .## 50% HDR .## - .##" 1431 PRINT USING A$ ,H1,H2,J1,J2 1440 J5=.75 1450 J8=J2/J1 1460 GOSUB 7000 1470 IF V7=0 THEN 1510 1480 H3=J1 1490 H4=J2 1500 GOTO 1520 1510A$=" .## - .## 75% HDR .## - .##" 1511 PRINT USING A$ ,H3,H4,J1,J2 1520 J5=.95 1530 J8=J2/J1 1540 GOSUB 7000 1550 IF V7=0 THEN 1600 1560 H5=J1 1565 IF H5>.01 THEN 1570 1568 H5=.01 1570 H6=J2 1575 IF H6<.99 THEN 1580 1577 H6=.99 1580 V7=0 1590 GOTO 770 1600 REM 1608 IF J2<.99 THEN 1620 1610 J2=.99 1620A$=" .## - .## 95% HDR .## - .##" 1621 PRINT USING A$ ,H5,H6,J1,J2 1630PRINT" ---------------------------------------------------------" 1650 PRINT "IF YOU WANT TO EVALUATE POSTERIOR DISTRIBUTION TYPE '1'." 1660 PRINT "IF NOT, TYPE '0'." 1670 GOSUB 9000 1680 IF O1=0 THEN 1740 1690 IF O1=1 THEN 1720 1700 PRINT "REENTER. MUST BE 0 OR 1." 1710 GOTO 1670 1720 A3(0)=A:A3(1)=B 1730 CLOSE 1: CHAIN "CMODB" 1740 CLOSE 1: CHAIN "RSTRT" 1750 E1=0 1760 A7=0 1770 GOSUB 5220 1780 E2=1 1790 A8=.999999 1800 A6=P1 1810 S5=(A6-A7)/(A8-A7) 1820 IF S5<.85 THEN 1840 1830 S5=.85 1840 IF S5>.15 THEN 1860 1850 S5=.15 1860 J2=E1+S5*(E2-E1) 1870 J1=0 1880 GOSUB 5000 1882 IF ABS(E1-E2)<.0001 THEN 1990 1890 IF ABS(P-P1)<.005 THEN 1990 1910 IF P>P1 THEN 1950 1920 E1=J2 1930 A7=P 1940 GOTO 1810 1950 A8=P 1960 E2=J2 1970 GOTO 1810 1980 GOTO 1990 1990 IF J2<.99 THEN 2010 2000 J2=.99 2010 RETURN 5000 REM *************************************************** 5005 REM BETA CDF ROUTINE 5010 REM INPUT A B J2 5015 REM OUTPUT P 5020 REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220 5025 DIM W(16),O(16) 5030 GOSUB 5340 5035 IF A+B>85 THEN 5280 5040 P=0 5045 C6=0 5050 IF A>1 THEN 5080 5055 C6=A 5060 C7=B 5065 A=C7 5070 B=C6 5075 J2=1-J2 5080 D0=(J2-J1)*.5 5085 D1=(J1+J2)*.5 5090 FOR I1=1 TO 16 5095 D9=D0*O(I1)+D1 5100 IF D9=0 THEN 5115 5105 IF D9=1 THEN 5115 5107 D9=LOG(D9)*(A-1)+LOG(1-D9)*(B-1) 5108 IF D9<-80 THEN 5115 5110 P=P+W(I1)*EXP(D9) 5115 NEXT I1 5120 P=P*F0 5125 P=P*D0 5130 IF C6=0 THEN 5155 5135 A=C6 5140 B=C7 5145 P=1-P 5150 J2=1-J2 5155 RETURN 5160 FOR I1=1 TO 16 5165 READ W(I1),O(I1) 5170 NEXT I1 5175 DATA 2.71525E-02,-.989401 5180 DATA 6.22535E-02,-.944575,9.51585E-02,-.865631 5185 DATA .124629,-.755404,.149596,-.617876 5190 DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02 5195 DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017 5200 DATA .149596,.617876,.124629,.755404 5205 DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02 5210 DATA .989401 5215 RETURN 5220 G9=A+B 5225 GOSUB 5850 5230 F0=G0 5235 G9=A 5240 GOSUB 5850 5245 F0=F0-G0 5250 G9=B 5255 GOSUB 5850 5260 F0=F0-G0 5265 IF A+B>85 THEN 5275 5270 F0=EXP(F0) 5275 RETURN 5280 W1=(B*J2)^(1/3) 5285 W2=(A*(1-J2))^(1/3) 5290 GOSUB 5325 5295 I1=P 5300 W1=(B*J1)^(1/3) 5305 W2=(A*(1-J1))^(1/3) 5310 GOSUB 5325 5315 P=I1-P 5320 RETURN 5325 Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A) 5330 GOSUB 8000 5335 RETURN 5340 REM 2/16/76 CHANGED TO ALL LOG 5345 D2=F0+(A-1)*LOG(J2)+(B-1)*LOG(1-J2) 5350 IF D2<-80 THEN 5370 5355 IF D2>85 THEN 5380 5360 D2=EXP(D2) 5365 RETURN 5370 D2=1.E-37 5375 RETURN 5380 D2=1.E+37 5385 RETURN 5390 REM 5395 REM END OF BETA CDF ROUTINE 5400 REM ******************************************************* 5850 REM **************************************************** 5852 REM LOG GAMMA ROUTINE 5853 REM INPUT G9 5854 REM OUTPUT G0 5860 G5=G9 5863 IF G9 <= 1.E+30 THEN 5872 5866 G0=1.E+38 5869 RETURN 5872 IF G9>1.E-09 THEN 5881 5875 G0=0 5878 RETURN 5881 IF G9<1.E+10 THEN 5890 5884 G0=G9*(LOG(G9)-1) 5887 RETURN 5890 G6=1 5893 IF 18.001 THEN 7042 7040 GOTO 7108 7042 IF P0>J5 THEN 7052 7044 J7=J8 7046 J8=J8*2 7048 P9=P0 7050 GOTO 7028 7052 J9=J7 7054 J0=J8 7056 J8=(J5-P9)/(P0-P9) 7058 IF J8>.15 THEN 7064 7060 J8=.15 7062 GOTO 7068 7064 IF J8<.85 THEN 7068 7066 J8=.85 7068 J8=J8*(J0-J9)+J9 7070 GOSUB 7096 7072 IF J8<-99 THEN 7080 7074 GOSUB 5000 7076 J3=P 7077 IF ABS(J1-J2)<.0001 THEN 7108 7078 IF ABS(J3-J5)>.001 THEN 7082 7080 GOTO 7108 7082 IF J3>J5 THEN 7090 7084 J9=J8 7086 P9=J3 7088 GOTO 7056 7090 J0=J8 7092 P0=J3 7094 GOTO 7056 7096 J1=(J8^J-1)/(J8^(J+1)-1) 7098 J2=J8*J1 7100 IF J2<1 THEN 7106 7102 J8=J8*.95 7104 GOTO 7100 7106 RETURN 7108 IF U0=0 THEN 7122 7110 J7=A 7112 A=B 7114 B=J7 7116 J7=J2 7118 J2=1-J1 7120 J1=1-J7 7122 RETURN 7124 REM END OF BETA HDR ROUTINE 7126 REM************************************************* 8000 REM ********************************************************** 8001 REM ROUTINE CALCULATES THE CDF FOR NORMAL DISTRIBUTION 8002 REM INPUT Y3 8003 REM OUTPUT P 8004 REM 8005 Y4=ABS(Y3) 8010 X1=X 8015 X=Y3 8020 T=1/(1+.231642*Y4) 8021 IF X*X/2<80 THEN 8025 8022 D=0 8023 GOTO 8030 8025 D=.398942*EXP(-X*X/2) 8030 C1=1.33027 8035 C2=1.82126 8040 C3=1.78148 8045 C4=.356564 8050 C5=.319382 8055 P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5) 8060 IF X >= 0 THEN 8070 8065 P=1-P 8070 X=X1 8075 RETURN 8076 REM 8077 REM END OF NORMAL CDF ROUTINE 8078 REM ********************************************************** 9000 REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED. 9005 INPUT O1 9015 IF O1=-9999 THEN 9025 9020 RETURN 9025 CLOSE 1: CHAIN "RSTRT" 9035 REM*************END ROUTINE 9999 END