10X$=SYS(CHR$(7)):R$=LEFT(X$,7):L$=RIGHT(X$,8) 20 REM************************************************************* 30 REM CMODG CMODG CMODG CMODG CMODG CMODG 40 REM*********************************************************** 50 REM 60 REM PRIOR FOR MULTINOMIAL DIRCHLET 70 REM 80 REM******************************************************* 90 DIM #1,A1(111),A2(111),A3(111) 100 OPEN R$ AS FILE 1 140A1(0)=16 150 DIM Q(10) 160 DIM R(10) 170 FOR K5=1 TO 10 180 R(K5)=0 190 NEXT K5 200 V9=0 210 K3=0 220 GOSUB 5160 230 X=0 240 V8=0 250 V7=0 260 S8=0 270 S7=0 280 K6=0 290 DIM P(10) 300 I2=1 310 PRINT L$ 320 PRINT " PRIOR DISTRIBUTION - MULTINOMIAL MODEL" 330 PRINT 340 PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A PRIOR DISTRIBUTION TO" 350 PRINT "YOUR BELIEFS ABOUT THE PROPORTIONS IN A MULTINOMIAL ANALYSIS." 360 PRINT "THE MODULE ATTEMPTS TO FIT A DIRICHLET DISTRIBUTION TO YOUR" 370 PRINT "BELIEFS ABOUT THE PROPORTIONS CONSIDERED JOINTLY, AND BETA" 380 PRINT "DISTRIBUTIONS TO YOUR BELIEFS ABOUT THE PROPORTIONS CONSIDERED" 390 PRINT "SEPARATELY. THERE CAN BE AS MANY AS 10 CATEGORIES." 400 PRINT 410 PRINT "HOW MANY CATEGORIES ARE THERE IN YOUR ANALYSIS"; 420 GOSUB 9000 430 S0=O1 440 PRINT 450 IF ABS(S0-INT(S0))<.0001 THEN 480 460 PRINT "REENTER. NUMBER OF CATEGORIES MUST BE AN INTEGER." 470 GOTO 420 480 IF S0<3 THEN 500 490 IF S0<11 THEN 630 500 PRINT "THE NUMBER OF CATEGORIES MUST BE GREATER THAN 2 BUT NOT GREATER" 510 PRINT "THAN 10. IF THERE ARE ONLY TWO CATEGORIES THEN YOU SHOULD SELECT" 520 PRINT "A BINARY MODEL FOR YOUR ANALYSIS." 530 PRINT 540 PRINT "IF YOU WANT TO RESPECIFY THE NUMBER OF CATEGORIES TYPE '1'." 550 PRINT "IF YOU WANT TO EXIT THE MODULE TYPE '0'." 560 GOSUB 9000 570 IF O1=1 THEN 410 580 IF O1=0 THEN 600 590 GOTO 410 600 CLOSE 1: CHAIN "RSTRT" 610 PRINT "REENTER. INPUT MUST 0 OR 1." 620 GOTO 560 630 GOSUB 700 640 GOTO 980 650 REM **************************************************************** 660 REM 670 REM ROUTINE FOR INPUTTING JOINT MODE ESTIMATE 680 REM 690 REM **************************************************************** 700 PRINT "INPUT YOUR ESTIMATES OF THE PROPORTION OF THE POPULATION IN" 710 PRINT "EACH CATEGORY. THESE ESTIMATES SHOULD SUM TO 1.0" 720 S1=0 730 FOR I1=1 TO S0 740 PRINT 750A$="CATEGORY ####, " 751PRINTUSING A$,I1; 770 GOSUB 9000 780 P(I1)=O1 790 IF P(I1) <= 0 THEN 820 800 IF P(I1) >= 1 THEN 820 810 GOTO 840 820 PRINT "REENTER. ESTIMATE MUST BE BETWEEN 0 AND 1." 830 GOTO 750 840 S1=S1+P(I1) 850 IF I1=1 THEN 870 852 IF I1=S0 THEN 870 860A$=" ACCUMULATED PROPORTION = #.###" 861PRINT USING A$,S1; 870 NEXT I1 890 IF ABS(S1-1)>10^(-10) THEN 910 900 RETURN 910 PRINT 920 PRINT "PROPORTIONS DO NOT SUM TO 1.0. RESPECIFY." 930 PRINT 940 GOTO 720 950 REM ************* END OF ROUTINE - INPUT JOINT MODE *************** 960 REM 970 REM 980 PRINT L$ 990 PRINT " SPECIFY THE 25TH AND 75TH PERCENTILES OF YOUR PRIOR" 1000 PRINT "DISTRIBUTION ON THE PROPORTION IN EACH OF THE CATEGORIES." 1010 PRINT "WE ARE ASSUMING YOUR ESTIMATE WAS A MEASURE OF THE CENTRAL" 1020 PRINT "TENDENCY OF THE PRIOR DISTRIBUTION." 1030 PRINT 1040 PRINT "CATEGORY ESTIMATE" 1050 FOR K5=1 TO S0 1060A$=" ## .##" 1061 PRINT USING A$ ,K5,P(K5) 1080 PRINT " 25TH"; 1090 GOSUB 9000 1100 Q1=O1 1110 IF O1 <= 0 THEN 1120 1115 IF O1= 1 THEN 1190 1185 IF O1>P(K5) THEN 1210 1190 PRINT "75TH MUST BE GREATER THAN ESTIMATE AND LESS THAN 1. RESPECIFY." 1200 GOTO 1160 1210 Q2=P(K5) 1220 GOSUB 1320 1230 K3=K3+W8 1240 NEXT K5 1250 W8=K3/S0 1252 IF W8<8000THEN 1260 1254 W8=8000 1260 GOTO 1510 1270 REM ***************************************************************** 1280 REM 1290 REM ESTIMATE A AND B PARAMETERS FROM PERCENTILES 1300 REM 1310 REM ******************************************************************* 1320 D1=SQR(Q2*(1-Q1))-SQR(Q1*(1-Q2)) 1330 D1=D1*D1 1340 D3=SQR(Q2*(1-Q3))-SQR(Q3*(1-Q2)) 1350 D3=D3*D3 1360 C=.056*(1/D1+1/D3) 1370 A=C*Q2+1/3 1380 IF A>1.15 THEN 1470 1390 A=2 1400 W8=A/P(K5) 1410 W7=(A-1)/P(K5)+S0 1420 IF W8>W7 THEN 1440 1430 W8=W7 1440 IF W8-A>1.15 THEN 1500 1450 A=A+1 1460 GOTO 1400 1470 B=C*(1-Q2)+1/3 1480 IF B <= 1.15 THEN 1390 1490 W8=A+B 1500 RETURN 1510 PRINT L$ 1520 PRINT " HERE ARE THE PERCENTILES OF THE MARGINAL BETA DISTRIBUTIONS" 1530 PRINT "FITTED TO YOUR SPECIFICATIONS." 1540 PRINT 1550A$=" HYPOTHETICAL SAMPLE SIZE (A)=######.##" 1551 PRINT USING A$ ,W8 1570 IF S0>8 THEN 1590 1580 PRINT 1600PRINT"CATEGORY JOINT ESTIMATE 25TH 50TH 75TH" 1610 REM ****************************************************************** 1620 FOR K5=1 TO S0 1630 R(K5)=.5*(P(K5)*W8+(W8-S0)*P(K5)+1) 1640 A=R(K5) 1650 B=W8-A 1660 GOSUB 5220 1670 P1=25 1680 GOSUB 3110 1690 Q(1)=J2 1700 P1=50 1710 GOSUB 3130 1720 Q(2)=J2 1730 P1=75 1740 GOSUB 3130 1750 IF J2<.99 THEN 1770 1760 J2=.99 1770 Q(3)=J2 1780A$=" ## .## .## .## .##" 1781 PRINT USING A$ ,K5,P(K5),Q(1),Q(2),Q(3) 1800 NEXT K5 1810 PRINT 1820 IF V8=1 THEN 1900 1830 PRINT "IF YOU DO NOT FEEL THAT THE HYPOTHETICAL SAMPLE SIZE (A) REFLECTS" 1840 PRINT "YOUR PRIOR INFORMATION ABOUT THE PROPORTIONS YOU CAN SPECIFY A" 1850 PRINT "DIFFERENT (A). THIS WILL NOT AFFECT THE JOINT ESTIMATE BUT WILL" 1860 PRINT "CHANGE THE MARGINAL PERCENTILES. A LARGER (A) WILL RESULT IN " 1870 PRINT "SMALLER INTERPERCENTILE DIFFERENCES AND A SMALLER (A) IN LARGER" 1880 PRINT "ONES." 1890 PRINT 1900 PRINT"IF YOU WANT TO CHANGE (A) TYPE NEW VALUE, ELSE '0'"; 1920 GOSUB 9000 1925 IF O1>9999 THEN 2210 1930 IF O1=0THEN 2240 1940 V8=1 1950 K7=1 1960 FOR K5=2TO S0 1970 IF P(K7)M7 THEN 2080 2070 M8=M7 2080 IF M8-A>1.15 THEN 2110 2090 A=A+1 2100 GOTO 2040 2110 IF M850 THEN 3220 3140 E2=A/(A+B) 3150 IF E2>(A-1)/(A+B-2) THEN 3170 3160 E2=(A-1)/(A+B-2) 3170 J1=0 3180 J2=E2 3190 GOSUB 5000 3200 A8=P 3210 GOTO 3240 3220 A8=.999999 3230 E2=1 3240 A6=P1/100 3250 S9=.85 3260 IF S9<(A6-A7)/(A8-A7) THEN 3280 3270 S9=(A6-A7)/(A8-A7) 3280 IF S9>.15 THEN 3300 3290 S9=.15 3300 J2=E1+S9*(E2-E1) 3310 J1=0 3320 GOSUB 5000 3322 IF ABS(E1-E2)<.0001THEN3410 3330 IF ABS(P-A6)<.001 THEN 3410 3340 IF P>A6 THEN 3380 3350 E1=J2 3360 A7=P 3370 GOTO 3250 3380 E2=J2 3390 A8=P 3400 GOTO 3250 3410 RETURN 3420 REM************************************************************* 3430 REM APPENDED GOSUBS FOLLOW 3440 REM 5000 REM BETA DISTRIBUTION 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 REM GOSUB 5340 DENSITY 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