10X$=SYS(CHR$(7)):R$=LEFT(X$,7):L$=RIGHT(X$,8) 13 REM **************************************************************** 15 REM CMOD6 CMOD6 CMOD6 CMOD6 CMOD6 CMOD6 16 REM ***************************************************************** 20 REM 30 REM POSTERIOR FOR TWO PARAMETER NORMAL 40 REM 50 REM*************************************************************** 60 GOSUB 6054 70 DIM #1,A1(111),A2(111),A3(111) 80 S9=0 90 OPEN R$ AS FILE 1 130A1(0)=6 140V1=A2(0):L1=A2(1):Q3=A2(2):M0=A2(3):Q4=A2(4):V8=A2(5):Q2=A2(6) 150 REM ******************************************************************** 160 REM Q4=MEDIAN OF PRIOR ON STANDARD DEVIATION 170 REM Q1=HYPO. SAMPLE SIZE ON " " 180 REM 190 REM Q3=MEDIAN OF CONDITIONAL PRIOR ON MEAN 200 REM Q2=HYPO. SAMPLE SIZE ON MEAN 210 REM 220 REM ********************************************************************* 230 PRINT L$ 240 PRINT " POSTERIOR ANALYSIS FOR THE TWO PARAMETER NORMAL MODEL" 250 PRINT 260 IF V1 <> 0 THEN 610 270 PRINT "INPUT THE DEGREES OF FREEDOM OF THE PRIOR DISTRIBUTION ON THE " 280 PRINT "STANDARD DEVIATION."; 290 GOSUB 9000 300 V1=O1 305 Q1=V1+1 310 PRINT 320 IF V1>2 THEN 350 330 PRINT "REENTER. DEGREES OF FREEDOM MUST BE GREATER THAN 2." 340 GOTO 270 350 PRINT "INPUT THE SCALE PARAMETER OF THE PRIOR DISTRIBUTION ON THE " 355 PRINT "STANDARD DEVIATION."; 360 GOSUB 9000 370 L1=O1 380 PRINT 390 IF L1>0 THEN 420 400 PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE." 410 GOTO 350 420 PRINT "INPUT THE MEAN OF PRIOR ON POPULATION MEAN."; 450 GOSUB 9000 460 Q3=O1 470 PRINT 480 PRINT "INPUT THE SCALE PARAMETER OF PRIOR ON THE MEAN."; 490 GOSUB 9000 500 M0=O1 510 PRINT 520 IF M0>0 THEN 550 530 PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE." 540 GOTO 480 550 REM 560 Q4=L1/SQR(V1-2/3) 570 Q2=M0/(L1*L1) 571 Q2=1/Q2 572 REM 580 REM WHILE SAMPLE DATA IS BEING INPUTTED WE ARE SETTING UP PRINT 590 REM FOR SUMMARY OF ANALYSIS 592 REM 610 PRINT "HOW MANY OBSERVATIONS ARE THERE IN YOUR SAMPLE?"; 620 G=V1 630 Q5=L1*L1 640 J5=.5 650 Q1=V1+1 660 GOSUB 9000 670 Q6=O1 680 GOSUB 7300 690 PRINT 700 IF Q6>=2THEN 730 710 PRINT "REENTER. THERE MUST BE AT LEAST 2 OBSERVATIONS." 720 GOTO 610 730 PRINT "WHAT IS THE MEAN OF YOUR SAMPLE?"; 732 REM 733 REM H1 AND H2 ARE THE ENDPOINTS OF PRIOR INV CHI 50% HDR 734 REM 740 H1=J1*L1 750 H2=J2*L1 760 GOSUB 9000 770 PRINT 780 J5=.75 790 Q8=O1 800 PRINT "WHAT IS THE STANDARD DEVIATION OF YOUR SAMPLE?"; 810 GOSUB 9000 820 PRINT 830 F=O1 840 IF F>0 THEN 870 850 PRINT "THE STANDARD DEVIATION MUST BE GREATER THAN 0." 860 GOTO 800 870 GOSUB 7300 872 REM 874 REM H3 AND H4 ARE ENDPOINTS OF PRIOR INV CHI 75% HDR 876 REM 880 H3=L1*J1 890 H4=L1*J2 892 REM 894 REM SEE NOVICK AND JACKSON PAGE 224 895 REM 898 REM POSTERIOR K 900 Q7=Q2+Q6 908 REM POSTERIOR H 910 Q9=(Q2*Q3+Q6*Q8)/Q7 918 REM S(SQUARED) PART OF POSTERIOR LAMBDA 920 F1=Q6*F*F 928 REM 3RD TERM IN POSTERIOR LAMBDA 930 F2=Q2*Q6*ABS(Q8-Q3)^2/Q7 938 REM ST. DEV. PART OF JOINT POS. MODE 940 F3=SQR((Q5+F1+F2)/(Q6+Q1+1)) 948 REM POSTERIOR NU (DEGREES OF FREEDOM) 950 F4=Q1+Q6-1 958 REM POSTERIOR LAMBDA 960 F5=Q5+F1+F2 968 REM POSTERIOR INV CHI MODE 970 F6=SQR(F5/(F4+1)) 978 REM POSTERIOR INV CHI MEAN 980 F8=SQR(F5/(F4-1.5)) 988 REM POSTERIOR INVERSE CHI MEDIAN 990 F7=SQR(F5/(F4-2/3)) 1000 PRINT L$ 1010 PRINT "THE JOINT POSTERIOR MODE FOR THE POPULATION MEAN AND" 1020 PRINT "STANDARD DEVIATION IS THE POINT ON THE PLANE AROUND WHICH" 1030 PRINT "THE PROBABILITY IS MOST HIGHLY CONCENTRATED. HERE IS THE " 1040 PRINT "JOINT MODE FOR YOUR POSTERIOR DISTRIBUTION." 1050 PRINT 1060A$=" MEAN ############.##" 1061 PRINT USING A$ ,Q9 1080 PRINT 1090A$=" STD. DEV. ############.##" 1091 PRINT USING A$ ,F3 1110 J5=.95 1120 GOSUB 7300 1125 REM H5 AND H6 PRIOR INVI CHI 75% HDR 1130 H5=J1*L1 1135 REM 1137 REM SWITCH TO COMPUTING POSTERIOR INV CHI HDRS 1138 REM 1140 H6=J2*L1 1150 G=F4 1160 L2=SQR(F5) 1170 GOSUB 7300 1180 I5=J1*L2 1190 I6=J2*L2 1200 J5=.75 1210 GOSUB 7300 1220 I4=L2*J2 1230 PRINT 1240 PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'."; 1250 I3=L2*J1 1260 GOSUB 9000 1270 J5=.5 1280 GOSUB 7300 1290 I1=J1*L2 1300 I2=J2*L2 1310 PRINT L$ 1320 PRINT "*** SUMMARY OF ANALYSIS ON THE STANDARD DEVIATION ****" 1330 PRINT 1340 IF S9<>0 THEN 1390 1350A$=" MARGINAL INVERSE CHI DISTRIBUTIONS" 1351PRINT A$ 1370 GOTO 1400 1390A$=" MARGINAL STUDENT'S DISTRIBUTIONS" 1391PRINT A$ 1400 PRINT 1410A$=" PRIOR POSTERIOR" 1411PRINT A$ 1430 PRINT "--------------------------------------------------------" 1515 REM 1516 REM S9=1 IF T DISTRIBUTION (POSTERIOR ON MEAN) IS BEING SET UP 1518 REM S9=0 IF INV CHI (POSTERIOR ON ST. DEV) 1519 REM 1520 IF S9<>0 THEN 1550 1530 J5=.5 1540 GOSUB 7300 1550A$="DEGREES OF FREEDOM ##########.## ###############.##" 1551 PRINT USING A$ ,V1,F4 1560 IF S9<>0 THEN 1600 1570 I8=J1*L2 1580 I2=J2*L2 1590 GOTO 1620 1592 REM ************************************************************** 1593 REM 1594 REM SET SCALE PARAMETERS FOR THE PRIOR AND POSTEIOR T'S 1595 REM 1600 L1=L1*L1/Q2 1605 L7=L2 1610 L2=F9*F9*(F4-2) 1620A$="SCALE PARAMETER #########.## #############.##" 1621 PRINT USING A$ ,L1,L2 1630 IF S9=1 THEN 1670 1632 REM 1633 REM ******************************************************************* 1640A$="MEAN #########.## #############.##" 1641 PRINT USING A$ ,L1/SQR(V1-3/2),L2/SQR(F4-3/2) 1650A$="MODE #########.## #############.##" 1651 PRINT USING A$ ,L1/SQR(V1+1),L2/SQR(F4+1) 1660 GOTO 1690 1670A$="MEDIAN #########.## #############.##" 1671 PRINT USING A$ ,Q3,Q9 1680 GOTO 1700 1690A$="MEDIAN #########.## #############.##" 1691 PRINT USING A$ ,Q4,L2/SQR(F4-2/3) 1700A$="50% HDR ###############.## -######.## #######.## -#######.##" 1701 PRINT USING A$ ,H1,H2,I8,I2 1710A$="75% HDR ###############.## -######.## #######.## -#######.##" 1711 PRINT USING A$ ,H3,H4,I3,I4 1720A$="95% HDR ###############.## -######.## #######.## -#######.##" 1721 PRINT USING A$ ,H5,H6,I5,I6 1740 PRINT "----------------------------------------------------------" 1750 IF S9=1 THEN 2160 1760 F9=SQR(F5/(Q7*(F4-2))) 1770 B=SQR(F5/Q7/F4) 1780 J5=.5 1790 GOSUB 7500 1800 J2=J2*B 1810 I8=Q9-J2 1820 I2=J2+Q9 1830 J5=.75 1840 GOSUB 7500 1850 J2=J2*B 1860 I3=Q9-J2 1870 I4=J2+Q9 1880 J5=.95 1890 GOSUB 7500 1900 J2=J2*B 1910 I5=Q9-J2 1920 I6=Q9+J2 1930 G=V1 1940 B=SQR(L1*L1/Q2/V1) 1950 GOSUB 7500 1960 J2=J2*B 1970 H5=Q3-J2 1980 H6=Q3+J2 1990 J5=.75 2000 GOSUB 7500 2010 J2=J2*B 2020 H4=Q3+J2 2030 H3=Q3-J2 2040 J5=.5 2050 GOSUB 7500 2060 J2=J2*B 2070 PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'."; 2080 GOSUB 9000 2090 H1=Q3-J2 2100 H2=J2+Q3 2110 PRINT L$ 2120 S9=1 2130 PRINT L$ 2140 PRINT "********* SUMMARY OF ANALYSIS ON THE MEAN ************" 2150 GOTO 1330 2160 PRINT "THIS COMPLETES THE POSTERIOR ANALYSIS." 2170 PRINT "IF YOU WANT TO EVALUATE THE POSTERIOR ON THE MEAN TYPE '1'." 2190 PRINT "IF YOU WANT TO EVALUATE THE POSTERIOR ON ST. DEV. TYPE '2'." 2210 PRINT "IF NOT, TYPE '0'." 2220 GOSUB 9000 2230 I=O1 2240 IF I=0 THEN 2380 2250 REM INFORMATION PASSED TO SUBSEQUENT MODULES 2260A3(0)=F4:A3(1)=L7:A3(2)=Q9:A3(3)=L2 2262 REM *************************************************************** 2270 REM 2280 REM F4 DEGREES OF FREEDOM FOR INVERSE CHI AND T - POSTERIORS 2290 REM L7 SCALE PARAMETER INV CHI 2300 REM Q9 MEAN OF POSTERIOR T 2310 REM L2 SCALE PARAMETER POST. T 2315 REM 2317 REM *************************************************************** 2320 IF I=2 THEN 2370 2330 IF I=1 THEN 2360 2340 PRINT "REENTER. INPUT MUST BE 0,1, OR 2." 2350 GOTO 2220 2360 CLOSE 1: CHAIN "CMODC" 2370 CLOSE 1: CHAIN "CMODD" 2380 CLOSE 1: CHAIN "RSTRT" 5500 REM ******************************************************** 5501 REM CHI-SQUARE CDF ROUTINE-LOWER TAIL 5502 REM INPUT G X 5503 REM OUTPUT P 5504 REM PRIOR GOSUB 5775 5505 REM 5506 REM 5510 IF G>30 THEN 5795 5515 T0=.5*G-1 5520 X=X*.5 5525 T=81.4983 5530 P=5.575E-35*(T+X)^T0 5535 T=69.9622 5540 P=P+4.0883E-30*(T+X)^T0 5545 T=61.0585 5550 P=P+2.45182E-26*(T+X)^T0 5555 T=53.6086 5560 P=P+3.60577E-23*(T+X)^T0 5565 T=47.1531 5570 P=P+2.01052E-20*(T+X)^T0 5575 T=41.4517 5580 P=P+5.35019E-18*(T+X)^T0 5585 T=36.3584 5590 P=P+7.8198E-16*(T+X)^T0 5595 T=31.776 5600 P=P+6.89418E-14*(T+X)^T0 5605 T=27.6359 5610 P=P+3.91774E-12*(T+X)^T0 5615 T=23.8873 5620 P=P+1.50701E-10*(T+X)^T0 5625 T=20.4915 5630 P=P+4.07286E-09*(T+X)^T0 5635 T=17.418 5640 P=P+7.96081E-08*(T+X)^T0 5645 T=14.6427 5650 P=P+1.15132E-06*(T+X)^T0 5655 T=12.1461 5660 P=P+1.25447E-05*(T+X)^T0 5665 T=9.9121 5670 P=P+1.04461E-04*(T+X)^T0 5675 T=7.92754 5680 P=P+6.72163E-04*(T+X)^T0 5685 T=6.18154 5690 P=P+3.36935E-03*(T+X)^T0 5695 T=4.66508 5700 P=P+.013226*(T+X)^T0 5705 T=3.37077 5710 P=P+4.07325E-02*(T+X)^T0 5715 T=2.29256 5720 P=P+9.81663E-02*(T+X)^T0 5725 T=1.4256 5730 P=P+.183323*(T+X)^T0 5735 T=.766097 5740 P=P+.258807*(T+X)^T0 5745 T=.311239 5750 P=P+.258774*(T+X)^T0 5755 T=5.90199E-02 5760 P=P+.142812*(T+X)^T0 5761IF X<80 THEN 5765 5762P=1 5763RETURN 5765 P=1-(T1*P*EXP(-X)) 5770 RETURN 5775 G9=.5*G 5778 IF G>30 THEN 5790 5780 GOSUB 5850 5785 T1=1/EXP(G0) 5790 RETURN 5795 X2=((X/G)^(1/3)-(1-2/9/G))/SQR(2/9/G) 5800 Y3=X2 5805 GOSUB 8000 5810 RETURN 5812 REM 5813 REM END OF CHI-SQUARE CDF ROUTINE 5815 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 1880 THEN 6044 6043 P=P+W(I1)*EXP(-D9) 6044 NEXT I1 6046 P=P*F0 6048 P=P*D0 6050 P=P+.5 6052 RETURN 6054 FOR I1=1 TO 16 6056 READ W(I1),O(I1) 6058 NEXT I1 6060 DATA 2.71525E-02,-.989401 6062 DATA 6.22535E-02,-.944575,9.51585E-02,-.865631 6064 DATA .124629,-.755404,.149596,-.617876 6066 DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02 6068 DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017 6070 DATA .149596,.617876,.124629,.755404 6072 DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02 6074 DATA .989401 6076 RETURN 6078 G9=(N+1)/2 6080 GOSUB 5850 6082 F0=G0 6084 G9=N/2 6086 GOSUB 5850 6088 F0=EXP(F0-G0)/SQR(3.14159*N) 6090 RETURN 6092 REM FOLLOWING FOR NU=1 6094 P=.5+1/3.14159*ATN(Y3) 6096 RETURN 6097 REM END OF STUDENT'S T CDF ROUTINE 6098 REM************************************************************* 7300 REM ******************************************************* 7301 REM INVERSE CHI HIGHEST DENSITY REGION ROUTINE 7302 REM INPUTS G J5 7303 REM OUTPUTS J1 J2 7304 REM 7305 IF G>30 THEN 7307 7306 GOSUB 5775 7307 J8=1 7308 GOSUB 7337 7309 X=1/(J1*J1) 7310 GOSUB 5500 7311 J3=P 7312 X=1/(J2*J2) 7313 GOSUB 5500 7314 J3=J3-P 7315 IF ABS(J3-J5)>.0001 THEN 7317 7316 RETURN 7317 IF J3>J5 THEN 7320 7318 J8=J8+1 7319 GOTO 7308 7320 J9=J8-1 7321 J0=J8 7322 J8=(J0+J9)/2 7323 GOSUB 7337 7324 X=1/(J1*J1) 7325 GOSUB 5500 7326 J3=P 7327 X=1/(J2*J2) 7328 GOSUB 5500 7329 J3=J3-P 7330 IF ABS(J3-J5)>.0001 THEN 7332 7331 RETURN 7332 IF J3>J5 THEN 7335 7333 J9=J8 7334 GOTO 7322 7335 J0=J8 7336 GOTO 7322 7337 J=J8*(1+EXP(-2*J8/(G+1)))/(1-EXP(-2*J8/(G+1))) 7338 J1=1/SQR(J+J8) 7339 J2=1/SQR(J-J8) 7340 RETURN 7341 REM 7342 REM END OF ROUTINE FOR INVERSE CHI HDR'S 7343 REM************************************************** 7500 REM ************************************************************ 7501 REM STUDENT'S T DISTRIBUTION HIGHEST DENSITY REGIONS 7502 REM INPUTS G J5 7503 REM J2 7504 REM 7505 Z8=.5 7506 N=G 7507 X9=1 7508 J1=0 7509 J2=X9 7510 GOSUB 6000 7511 P=2*P-1 7512 Z9=P 7513 IF P>J5 THEN 7517 7514 X9=X9+1 7515 Z8=Z9 7516 GOTO 7508 7517 X0=X9-1 7518 X2=X9 7519 X9=X0+(J5-Z8)*(X2-X0)/(Z9-Z8) 7520 J1=0 7521 J2=X9 7522 GOSUB 6000 7523 P=2*P-1 7524 IF ABS(P-J5)<.0001 THEN 7534 7525 IF P= 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