Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-04 - 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