Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod6.bas
There are 2 other files named cmod6.bas in the archive. Click here to see a list.
00020REM ****************************************************************
00030REM CMOD6 CMOD6 CMOD6 CMOD6 CMOD6 CMOD6
00040REM *****************************************************************
00050REM
00060REM POSTERIOR FOR TWO PARAMETER NORMAL
00070REM
00080REM***************************************************************
00090GOSUB 6060
00100FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00110S9=0
00150RESTORE#1
00151 INPUT# 1,I1,I2,I3
00160SCRATCH#1
00161 PRINT # 1,6,I2,I3
00170RESTORE#2
00171 INPUT# 2,V1,L1,Q3,M0,Q4,V8,Q2
00180REM ********************************************************************
00190REM Q4=MEDIAN OF PRIOR ON STANDARD DEVIATION
00200REM Q1=HYPO. SAMPLE SIZE ON " "
00210REM
00220REM Q3=MEDIAN OF CONDITIONAL PRIOR ON MEAN
00230REM Q2=HYPO. SAMPLE SIZE ON MEAN
00240REM
00250REM *********************************************************************
00260PRINT L$
00270PRINT " POSTERIOR ANALYSIS FOR THE TWO PARAMETER NORMAL MODEL"
00280PRINT
00290IF V1 <> 0 THEN 660
00300PRINT "INPUT THE DEGREES OF FREEDOM OF THE PRIOR DISTRIBUTION ON THE "
00310PRINT "STANDARD DEVIATION.";
00320GOSUB 9000
00330V1=O1
00340Q1=V1+1
00350PRINT
00360IF V1>2 THEN 390
00370PRINT "REENTER. DEGREES OF FREEDOM MUST BE GREATER THAN 2."
00380GOTO 300
00390PRINT "INPUT THE SCALE PARAMETER OF THE PRIOR DISTRIBUTION ON THE "
00400PRINT "STANDARD DEVIATION.";
00410GOSUB 9000
00420L1=O1
00430PRINT
00440IF L1>0 THEN 470
00450PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00460GOTO 390
00470PRINT "INPUT THE MEAN OF PRIOR ON POPULATION MEAN.";
00480GOSUB 9000
00490Q3=O1
00500PRINT
00510PRINT "INPUT THE SCALE PARAMETER OF PRIOR ON THE MEAN.";
00520GOSUB 9000
00530M0=O1
00540PRINT
00550IF M0>0 THEN 580
00560PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00570GOTO 510
00580REM
00590Q4=L1/SQR(V1-2/3)
00600Q2=M0/(L1*L1)
00610Q2=1/Q2
00620REM
00630REM WHILE SAMPLE DATA IS BEING INPUTTED WE ARE SETTING UP PRINT
00640REM FOR SUMMARY OF ANALYSIS
00650REM
00660PRINT "HOW MANY OBSERVATIONS ARE THERE IN YOUR SAMPLE?";
00670G=V1
00680Q5=L1*L1
00690J5=.5
00700Q1=V1+1
00710GOSUB 9000
00720Q6=O1
00730GOSUB 7300
00740PRINT
00750IF Q6 >= 2 THEN 780
00760PRINT "REENTER. YOU MUST HAVE AT LEAST 2 OBSERVATIONS."
00770GOTO 660
00780PRINT "WHAT IS THE MEAN OF YOUR SAMPLE?";
00790REM
00800REM H1 AND H2 ARE THE ENDPOINTS OF PRIOR INV CHI 50% HDR
00810REM
00820H1=J1*L1
00830H2=J2*L1
00840GOSUB 9000
00850PRINT
00860J5=.75
00870Q8=O1
00880PRINT "WHAT IS THE STANDARD DEVIATION OF YOUR SAMPLE?";
00890GOSUB 9000
00900PRINT
00910F=O1
00920IF F>0 THEN 950
00930PRINT "THE STANDARD DEVIATION MUST BE GREATER THAN 0."
00940GOTO 880
00950GOSUB 7300
00960REM
00970REM H3 AND H4 ARE ENDPOINTS OF PRIOR INV CHI 75% HDR
00980REM
00990H3=L1*J1
01000H4=L1*J2
01010REM
01020REM SEE NOVICK AND JACKSON PAGE 224
01030REM
01040REM POSTERIOR K
01050Q7=Q2+Q6
01060REM POSTERIOR H
01070Q9=(Q2*Q3+Q6*Q8)/Q7
01080REM S(SQUARED) PART OF POSTERIOR LAMBDA
01090F1=Q6*F*F
01100REM 3RD TERM IN POSTERIOR LAMBDA
01110F2=Q2*Q6*ABS(Q8-Q3)**2/Q7
01120REM ST. DEV. PART OF JOINT POS. MODE
01130F3=SQR((Q5+F1+F2)/(Q6+Q1+1))
01140REM POSTERIOR NU (DEGREES OF FREEDOM)
01150F4=Q1+Q6-1
01160REM POSTERIOR LAMBDA
01170F5=Q5+F1+F2
01180REM POSTERIOR INV CHI MODE
01190F6=SQR(F5/(F4+1))
01200REM POSTERIOR INV CHI MEAN
01210F8=SQR(F5/(F4-1.5))
01220REM POSTERIOR INVERSE CHI MEDIAN
01230F7=SQR(F5/(F4-2/3))
01240PRINT L$
01250PRINT "THE JOINT POSTERIOR MODE FOR THE POPULATION MEAN AND"
01260PRINT "STANDARD DEVIATION IS THE POINT ON THE PLANE AROUND WHICH"
01270PRINT "THE PROBABILITY IS MOST HIGHLY CONCENTRATED. HERE IS THE "
01280PRINT "JOINT MODE FOR YOUR POSTERIOR DISTRIBUTION."
01290PRINT
01300 PRINT USING 1310,Q9
01310: MEAN ########.##
01320PRINT
01330 PRINT USING 1340,F3
01340: STANDARD DEVIATION########.##
01350J5=.95
01360GOSUB 7300
01370REM H5 AND H6 PRIOR INVI CHI 75% HDR
01380H5=J1*L1
01390REM
01400REM SWITCH TO COMPUTING POSTERIOR INV CHI HDRS
01410REM
01420H6=J2*L1
01430G=F4
01440L2=SQR(F5)
01450GOSUB 7300
01460I5=J1*L2
01470I6=J2*L2
01480J5=.75
01490GOSUB 7300
01500I4=L2*J2
01510PRINT
01520PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01530I3=L2*J1
01540GOSUB 9000
01550J5=.5
01560GOSUB 7300
01570I1=J1*L2
01580I2=J2*L2
01590PRINT L$
01600PRINT "*** SUMMARY OF ANALYSIS ON THE STANDARD DEVIATION ****"
01610PRINT
01620IF S9 <> 0 THEN 1670
01630REM
01640PRINT" MARGINAL INVERSE CHI DISTRIBUTIONS"
01650GOTO 1680
01660: MARGINAL STUDENT'S DISTRIBUTIONS
01670PRINT" MARGINAL STUDENT'S DISTRIBUTIONS"
01680PRINT
01690PRINT" PRIOR POSTERIOR"
01700: PRIOR POSTERIOR
01710PRINT "--------------------------------------------------------"
01720:DEGREES OF FREEDOM##########.## ########.##
01730:SCALE PARAMETER#############.## ########.##
01740:MEAN ########.## ########.##
01750:MODE ########.## ########.##
01760:MEDIAN ########.## ########.##
01770:50% HDR ########.## TO########.## ########.## TO########.##
01780:75% HDR ########.## TO########.## ########.## TO########.##
01790:95% HDR ########.## TO########.## ########.## TO########.##
01800REM
01810REM S9=1 IF T DISTRIBUTION (POSTERIOR ON MEAN) IS BEING SET UP
01820REM S9=0 IF INV CHI (POSTERIOR ON ST. DEV)
01830REM
01840IF S9 <> 0 THEN 1870
01850J5=.5
01860GOSUB 7300
01870 PRINT USING 1720,V1,F4
01880IF S9 <> 0 THEN 1960
01890I8=J1*L2
01900I2=J2*L2
01910GOTO 1990
01920REM **************************************************************
01930REM
01940REM SET SCALE PARAMETERS FOR THE PRIOR AND POSTEIOR T'S
01950REM
01960L1=L1*L1/Q2
01970L7=L2
01980L2=F9*F9*(F4-2)
01990 PRINT USING 1730,L1,L2
02000IF S9=1 THEN 2060
02010REM
02020REM *******************************************************************
02030 PRINT USING 1740,L1/SQR(V1-3/2),L2/SQR(F4-3/2)
02040 PRINT USING 1750,L1/SQR(V1+1),L2/SQR(F4+1)
02050GOTO 2080
02060 PRINT USING 1760,Q3,Q9
02070GOTO 2090
02080 PRINT USING 1760,Q4,L2/SQR(F4-2/3)
02090 PRINT USING 1770,H1,H2,I8,I2
02100 PRINT USING 1780,H3,H4,I3,I4
02110 PRINT USING 1790,H5,H6,I5,I6
02120PRINT "----------------------------------------------------------"
02130IF S9=1 THEN 2540
02140F9=SQR(F5/(Q7*(F4-2)))
02150B=SQR(F5/Q7/F4)
02160J5=.5
02170GOSUB 7500
02180J2=J2*B
02190I8=Q9-J2
02200I2=J2+Q9
02210J5=.75
02220GOSUB 7500
02230J2=J2*B
02240I3=Q9-J2
02250I4=J2+Q9
02260J5=.95
02270GOSUB 7500
02280J2=J2*B
02290I5=Q9-J2
02300I6=Q9+J2
02310G=V1
02320B=SQR(L1*L1/Q2/V1)
02330GOSUB 7500
02340J2=J2*B
02350H5=Q3-J2
02360H6=Q3+J2
02370J5=.75
02380GOSUB 7500
02390J2=J2*B
02400H4=Q3+J2
02410H3=Q3-J2
02420J5=.5
02430GOSUB 7500
02440J2=J2*B
02450PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
02460GOSUB 9000
02470H1=Q3-J2
02480H2=J2+Q3
02490PRINT L$
02500S9=1
02510PRINT L$
02520PRINT "********* SUMMARY OF ANALYSIS ON THE MEAN ************"
02530GOTO 1610
02540PRINT "THIS COMPLETES THE POSTERIOR ANALYSIS."
02550PRINT "IF YOU WANT TO EVALUATE THE POSTERIOR ON THE MEAN TYPE '1'."
02560PRINT "IF YOU WANT TO EVALUATE THE POSTERIOR ON ST. DEV. TYPE '2'."
02570PRINT "IF NOT, TYPE '0'."
02580GOSUB 9000
02590I=O1
02600IF I=0 THEN 2770
02610REM INFORMATION PASSED TO SUBSEQUENT MODULES
02620SCRATCH#3
02621 PRINT # 3,F4,L7,Q9,L2
02630REM ***************************************************************
02640REM
02650REM F4 DEGREES OF FREEDOM FOR INVERSE CHI AND T - POSTERIORS
02660REM L7 SCALE PARAMETER INV CHI
02670REM Q9 MEAN OF POSTERIOR T
02680REM L2 SCALE PARAMETER POST. T
02690REM
02700REM ***************************************************************
02710IF I=2 THEN 2760
02720IF I=1 THEN 2750
02730PRINT "REENTER. INPUT MUST BE 0,1, OR 2."
02740GOTO 2580
02750CHAIN "CMODC"
02760CHAIN "CMODD"
02770CHAIN "RSTRT"
05500REM ********************************************************
05501REM CHI-SQUARE CDF ROUTINE-LOWER TAIL
05502REM INPUT G X
05503REM OUTPUT P
05504REM PRIOR GOSUB 5775
05505REM
05506REM
05510IF G>30 THEN 5795
05515T0=.5*G-1
05520X=X*.5
05525T=81.4983
05530P=5.575E-35*(T+X)**T0
05535T=69.9622
05540P=P+4.0883E-30*(T+X)**T0
05545T=61.0585
05550P=P+2.45182E-26*(T+X)**T0
05555T=53.6086
05560P=P+3.60577E-23*(T+X)**T0
05565T=47.1531
05570P=P+2.01052E-20*(T+X)**T0
05575T=41.4517
05580P=P+5.35019E-18*(T+X)**T0
05585T=36.3584
05590P=P+7.8198E-16*(T+X)**T0
05595T=31.776
05600P=P+6.89418E-14*(T+X)**T0
05605T=27.6359
05610P=P+3.91774E-12*(T+X)**T0
05615T=23.8873
05620P=P+1.50701E-10*(T+X)**T0
05625T=20.4915
05630P=P+4.07286E-09*(T+X)**T0
05635T=17.418
05640P=P+7.96081E-08*(T+X)**T0
05645T=14.6427
05650P=P+1.15132E-06*(T+X)**T0
05655T=12.1461
05660P=P+1.25447E-05*(T+X)**T0
05665T=9.9121
05670P=P+1.04461E-04*(T+X)**T0
05675T=7.92754
05680P=P+6.72163E-04*(T+X)**T0
05685T=6.18154
05690P=P+3.36935E-03*(T+X)**T0
05695T=4.66508
05700P=P+.013226*(T+X)**T0
05705T=3.37077
05710P=P+4.07325E-02*(T+X)**T0
05715T=2.29256
05720P=P+9.81663E-02*(T+X)**T0
05725T=1.4256
05730P=P+.183323*(T+X)**T0
05735T=.766097
05740P=P+.258807*(T+X)**T0
05745T=.311239
05750P=P+.258774*(T+X)**T0
05755T=5.90199E-02
05760P=P+.142812*(T+X)**T0
05765P=1-(T1*P*EXP(-X))
05770RETURN
05775G9=.5*G
05778IF G>30 THEN 5790
05780GOSUB 5850
05785T1=1/EXP(G0)
05790RETURN
05795X2=((X/G)**(1/3)-(1-2/9/G))/SQR(2/9/G)
05800Y3=X2
05805GOSUB 8000
05810RETURN
05812REM
05813REM END OF CHI-SQUARE CDF ROUTINE
05815REM *********************************************************
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 ****************************************************
06000REM****************************************************************
06002REM STUDENT'S T CDF ROUTINE
06004REM INPUT G J2
06006REM OUTPUT P
06008REM PRIOR GOSUB 6060
06010IF J2<6 THEN 6016
06012P=1
06014GOTO 6058
06016 DIM W(16),O(16)
06018Y3=J2
06020IF G <= 120 THEN 6026
06022GOSUB 8000
06024GOTO 6058
06026IF G=1 THEN 6098
06028J1=0
06030N=G
06032P=0
06034GOSUB 6084
06036D0=(J2-J1)*.5
06038D1=(J1+J2)*.5
06040FOR I1=1 TO 16
06042D9=D0*O(I1)+D1
06044IF D9=0 THEN 6050
06046IF D9=1 THEN 6050
06048P=P+W(I1)*(EXP(-(N+1)/2*LOG(1+D9*D9/N)))
06050NEXT I1
06052P=P*F0
06054P=P*D0
06056P=P+.5
06058RETURN
06060FOR I1=1 TO 16
06062READ W(I1),O(I1)
06064NEXT I1
06066DATA 2.71525E-02,-.989401
06068DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
06070DATA .124629,-.755404,.149596,-.617876
06072DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
06074DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
06076DATA .149596,.617876,.124629,.755404
06078DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
06080DATA .989401
06082RETURN
06084G9=(N+1)/2
06086GOSUB 5850
06088F0=G0
06090G9=N/2
06092GOSUB 5850
06094F0=EXP(F0-G0)/SQR(3.14159*N)
06096RETURN
06098REM FOLLOWING FOR NU=1
06100P=.5+1/3.14159*ATN(Y3)
06102RETURN
06104REM END OF STUDENT'S T CDF ROUTINE
06106REM*************************************************************
07300REM *******************************************************
07301REM INVERSE CHI HIGHEST DENSITY REGION ROUTINE
07302REM INPUTS G J5
07303REM OUTPUTS J1 J2
07304REM
07305IF G>30 THEN 7307
07306GOSUB 5775
07307J8=1
07308GOSUB 7337
07309X=1/(J1*J1)
07310GOSUB 5500
07311J3=P
07312X=1/(J2*J2)
07313GOSUB 5500
07314J3=J3-P
07315IF ABS(J3-J5)>.0001 THEN 7317
07316RETURN
07317IF J3>J5 THEN 7320
07318J8=J8+1
07319GOTO 7308
07320J9=J8-1
07321J0=J8
07322J8=(J0+J9)/2
07323GOSUB 7337
07324X=1/(J1*J1)
07325GOSUB 5500
07326J3=P
07327X=1/(J2*J2)
07328GOSUB 5500
07329J3=J3-P
07330IF ABS(J3-J5)>.0001 THEN 7332
07331RETURN
07332IF J3>J5 THEN 7335
07333J9=J8
07334GOTO 7322
07335J0=J8
07336GOTO 7322
07337J=J8*(1+EXP(-2*J8/(G+1)))/(1-EXP(-2*J8/(G+1)))
07338J1=1/SQR(J+J8)
07339J2=1/SQR(J-J8)
07340RETURN
07341REM
07342REM END OF ROUTINE FOR INVERSE CHI HDR'S
07343REM**************************************************
07500REM ************************************************************
07501REM STUDENT'S T DISTRIBUTION HIGHEST DENSITY REGIONS
07502REM INPUTS G J5
07503REM J2
07504REM
07505Z8=.5
07506N=G
07507X9=1
07508J1=0
07509J2=X9
07510GOSUB 6000
07511P=2*P-1
07512Z9=P
07513IF P>J5 THEN 7517
07514X9=X9+2
07515Z8=Z9
07516GOTO 7508
07517X0=X9-2
07518X2=X9
07519X9=X0+(J5-Z8)*(X2-X0)/(Z9-Z8)
07520J1=0
07521J2=X9
07522GOSUB 6000
07523P=2*P-1
07524IF ABS(P-J5)<.0001 THEN 7541
07525IF P<J5 THEN 7538
07526X2=X9
07527Z9=P
07528X9=(J5-Z8)/(Z9-Z8)
07529GOTO 7536
07530IF X9<.85 THEN 7533
07531X9=X0+.85*(X2-X0)
07532GOTO 7520
07533IF X9>.15 THEN 7536
07534X9=X0+.15*(X2-X0)
07535GOTO 7520
07536X9=X0+X9*(X2-X0)
07537GOTO 7520
07538X0=X9
07539Z8=P
07540GOTO 7528
07541J2=X9
07542RETURN
07543REM
07544REM END OF STUDENT'S T HDR ROUTINE
07545REM *********************************************************
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)
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