Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod5.bas
There are 2 other files named cmod5.bas in the archive. Click here to see a list.
00020REM *********************************************************************
00030REM CMOD5 CMOD5 CMOD5 CMOD5 CMOD5 CMOD5 CMOD5
00040REM *********************************************************************
00050REM
00060REM PRIOR FOR MEAN OF TWO PARAMETER NORMAL
00070REM
00080REM******************************************************************
00090FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00130RESTORE#1
00131 INPUT# 1,I1,I2,I3
00140SCRATCH#1
00141 PRINT # 1,5,I2,I3
00150GOSUB 6060
00160S9=0
00170S6=0
00180 DIM F(4,2)
00190X=0
00200RESTORE#2
00201 INPUT# 2,N1,L1,O1,O1,Q
00210REM *********************************************************************
00220REM NOTATION N1=DEGREES OF FREEDOM OF INVERSE CHI PRIOR ON ST> DEV>
00230REM L1=SCALE PARAMETER
00240REM G=HYPOTHETICAL SAMPLE SIZE
00250REM Q=MEDIAN
00260REM****************************************************************
00270PRINT L$
00280PRINT " PRIOR DISTRIBUTION ON THE MEAN"
00290PRINT
00300PRINT "THIS MODULE WILL ASSIST YOU IN SPECIFYING A PRIOR DISTRIBUTION"
00310PRINT "ON THE MEAN OF A NORMAL DISTRIBUTION."
00320PRINT
00330IF N1=0 THEN 360
00340G=N1+1
00350GOTO 560
00360PRINT "INPUT THE DEGREES OF FREEDOM OF YOUR PRIOR DISTRIBUTION ON"
00370PRINT "THE STANDARD DEVIATION.";
00380GOSUB 9000
00390N1=O1
00400G=N1+1
00410PRINT
00420IF N1>2 THEN 460
00430PRINT
00440PRINT "REENTER. DEGREEES OF FREEDOM MUST BE GREATER THAN 2."
00450GOTO 360
00460PRINT "INPUT THE SCALE PARAMETER.";
00470GOSUB 9000
00480L1=O1
00490G=O1+1
00500PRINT
00510IF L1>0 THEN 540
00520PRINT "REENTER. SCALE PARAMETER MUST BE POSITIVE."
00530GOTO 460
00540Q=L1/SQR(N1-2/3)
00550PRINT L$
00560:SUPPOSE THE POPULATION STANDARD DEVIATION IS######.##.
00570 PRINT USING 560,Q
00580PRINT "SPECIFY THE 25TH, 50TH, AND 90TH PERCENTILES OF YOUR PRIOR"
00590PRINT "DISTRIBUTION ON THE POPULATION MEAN."
00600PRINT
00610PRINT "SPECIFY 25TH. YOUR BETTING ODDS ARE 3 TO 1 THAT THE MEAN IS MORE"
00620PRINT "THAN THIS VALUE.";
00630GOSUB 9000
00640Q1=O1
00650PRINT
00660PRINT "SPECIFY 50TH. YOUR ODDS ARE EVEN THAT THE MEAN IS LESS THAN THIS"
00670PRINT "VALUE.";
00680GOSUB 9000
00690IF O1 <= Q1 THEN 780
00700Q2=O1
00710PRINT
00720PRINT "SPECIFY 90TH. YOUR ODDS ARE 9 TO 1 THAT THE MEAN IS LESS THAN THIS"
00730PRINT "VALUE.";
00740GOSUB 9000
00750IF O1 <= Q2 THEN 780
00760Q3=O1
00770GOTO 830
00780PRINT "90TH MUST BE LARGER THAN 50TH AND 50TH LARGER THAN 25TH."
00790GOTO 600
00800REM************************************************************
00810REM STANDARD DEVIATION IS ESTIMATED IN DIFFERENT WAYS
00820REM **********************************************************************
00830K5=1
00840PRINT L$
00850PRINT "HERE ARE THE PERCENTILES OF FOUR NORMAL DISTRIBUTIONS FITTED"
00860PRINT "TO YOUR PERCENTILE SPECIFICATIONS."
00870PRINT
00880REM
00890PRINT" 10TH 25TH 50TH 75TH 90TH"
00900M=Q2
00910S=(Q2-Q1)/.6745
00920GOSUB 1270
00930S=(Q3-Q2)/1.282
00940K5=2
00950GOSUB 1270
00960S=(Q3-Q1)/1.957
00970M=Q3-1.282*S
00980K5=3
00990GOSUB 1270
01000M=.5*(M+Q2)
01010K5=4
01020GOSUB 1270
01030PRINT
01040PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE "
01050PRINT "WHICH DISTRIBUTION MOST CLOSELY CORRESPONDS TO YOUR PRIOR"
01060PRINT "BELIEFS. YOU CAN EITHER TENTATIVELY ACCEPT ONE OF THESE"
01070PRINT "DISTRIBUTIONS OR RESPECIFY THE PERCENTILES."
01080PRINT
01090PRINT "IF YOU WANT ONE OF THE DISTRIBUTIONS TYPE ITS NUMBER."
01100PRINT "IF YOU WANT TO RESPECIFY THE PERCENTILES TYPE '0'."
01110GOSUB 9000
01120IF O1=0 THEN 1380
01130IF O1=1 THEN 1200
01140IF O1=2 THEN 1200
01150IF O1=3 THEN 1200
01160IF O1=4 THEN 1200
01170PRINT "REENTER. INPUT MUST BE 0 OR NUMBER OF DISTRIBUTION.";
01180GOSUB 9000
01190GOTO 1120
01200M=F(O1,1)
01210S=F(O1,2)
01220GOTO 1530
01230PRINT
01240PRINT " 10TH 25TH 50TH 75TH 90TH"
01250RETURN
01260GOTO 1290
01270F1=M-1.2816*S
01280F2=M-.6745*S
01290F3=M
01300F4=M+.6745*S
01310F5=M+1.2816*S
01320 PRINT USING 1330,K5,F1,F2,F3,F4,F5
01330:######## ######.## ######.## ######.## ######.## ######.##
01340: 10TH 25TH 50TH 75TH 90TH
01350F(K5,1)=M
01360F(K5,2)=S
01370RETURN
01380PRINT
01390PRINT "RESPECIFY YOUR PERCENTILES."
01400PRINT "25TH.";
01410GOSUB 9000
01420Q1=O1
01430PRINT "50TH.";
01440GOSUB 9000
01450Q2=O1
01460PRINT "90TH.";
01470GOSUB 9000
01480Q3=O1
01490IF Q1 >= Q2 THEN 1510
01500IF Q3>Q2 THEN 830
01510PRINT "REENTER. PERCENTILES MUST BE INCREASING."
01520GOTO 1390
01530REM ****************************************************************
01540REM **** SWITCH FROM PERCENTILES TO HYPO SAMPLE SIZE ********
01550REM ******************************************************************
01560N2=(Q/S)*(Q/S)
01570PRINT L$
01580IF S9=0 THEN 1620
01590PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE NORMAL"
01600PRINT "DISTRIBUTION YOU ARE NOW CONSIDERING."
01610GOTO 1670
01620PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE NORMAL"
01630PRINT "DISTRIBUTION YOU ARE NOW CONSIDERING. THIS IS A"
01640PRINT "CONDITIONAL DISTRIBUTION SINCE IT IS ASSUMED THAT"
01650 PRINT USING 1660,Q
01660:THE POPULATION STANDARD DEVIATION IS####.##.
01670F2=M
01680PRINT
01690 PRINT USING 1700,N2
01700:HYPOTHETICAL SAMPLE SIZE (M)######.##
01710 PRINT USING 1720,F2
01720:MEAN=MODE=MEDIAN #########.##
01730S=Q/SQR(N2)
01740 PRINT USING 1770,S
01750F4=F2+.6745*S
01760F1=F2-.6745*S
01770:STANDARD DEVIATION #########.##
01780F5=F2+1.2816*S
01790F6=F2-1.2816*S
01800:10TH PERCENTILE #########.##
01810 PRINT USING 1800,F6
01820:25TH PERCENTILE #########.##
01830 PRINT USING 1820,F1
01840:75TH PERCENTILE #########.##
01850 PRINT USING 1840,F4
01860:90TH PERCENTILE #########.##
01870 PRINT USING 1860,F5
01890PRINT
01900IF S9=1 THEN 2230
01910S9=1
01920PRINT "IF YOU DO NOT FEEL THAT THIS VALUE OF M REFLECTS YOUR"
01930PRINT "PRIOR INFORMATION ABOUT THE MEAN YOU CAN SPECIFY A"
01940PRINT "DIFFERENT M. A SMALLER M WILL GIVE LONGER INTERVALS"
01950PRINT "AND A LARGER M SHORTER INTERVALS."
01960PRINT
01970PRINT "IF YOU WANT TO CHANGE M TYPE THE NEW VALUE ELSE '0'.";
01980GOSUB 9000
01990IF O1=0 THEN 2050
02000IF O1>0 THEN 2030
02010PRINT "INPUT AGAIN."
02020GOTO 1980
02030N2=O1
02040GOTO 1570
02050PRINT
02060IF S6=1 THEN 2100
02070PRINT "YOU CAN CHANGE THE CENTERING OF THE DISTRIBUTION BY"
02080PRINT "SPECIFYING A DIFFERENT MEDIAN. THIS WILL NOT AFFECT"
02090PRINT "THE HYPOTHETICAL SAMPLE SIZE."
02100PRINT
02110PRINT "IF YOU WANT TO SPECIFY A DIFFERENT MEDIAN TYPE '1'."
02120PRINT "IF YOU DO NOT TYPE '0'."
02130GOSUB 9000
02140IF O1=1 THEN 2180
02150IF O1=0 THEN 2310
02160PRINT "REENTER. INPUT MUST BE 0 OR 1."
02170GOTO 2130
02180PRINT "SPECIFY DIFFERENT MEDIAN.";
02190S6=1
02200GOSUB 9000
02210M=O1
02220GOTO 1570
02230PRINT "IF YOU WANT TO CHANGE M TYPE THE NEW VALUE ELSE '0'.";
02240GOSUB 9000
02250IF O1=0 THEN 2050
02260IF O1<0 THEN 2290
02270N2=O1
02280GOTO 1570
02290PRINT "REENTER. M MUST BE GREATER THAN 0."
02300GOTO 2240
02310PRINT L$
02320PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE PRIOR MARGINAL"
02330PRINT "DISTRIBUTION ON THE MEAN."
02340PRINT
02350PRINT " STUDENT'S T DISTRIBUTION"
02360:DEGREES OF FREEDOM #########.##
02370:STANDARD DEVIATION #########.##
02380:SCALE PARAMETER #########.##
02390PRINT
02400:50% HDR #########.## TO#########.##
02410:75% HDR #########.## TO#########.##
02420:95% HDR #########.## TO#########.##
02430 PRINT USING 2360,N1
02440 PRINT USING 2380,L1*L1/N2
02450 PRINT USING 1720,M
02460S=L1/SQR(N2*(N1-2))
02470 PRINT USING 2370,S
02480J5=.5
02490GOSUB 7500
02500GOSUB 2520
02510GOTO 2560
02520J2=J2*L1/SQR(N2*N1)
02530J1=M-J2
02540J2=M+J2
02550RETURN
02560 PRINT USING 2400,J1,J2
02570J5=.75
02580GOSUB 7500
02590GOSUB 2520
02600 PRINT USING 2410,J1,J2
02610J5=.95
02620GOSUB 7500
02630GOSUB 2520
02640 PRINT USING 2420,J1,J2
02650PRINT
02660PRINT "THIS COMPLETES THE SPECIFICATION OF A PRIOR ON THE MEAN."
02670PRINT "IF YOU ARE NOT GOING TO CONTINUE ANALYSIS, YOU NOW"
02680PRINT "SHOULD COPY DOWN THE DEGREES OF FREEDOM, MEAN, AND"
02690PRINT "SCALE PARAMETER OF THE PRIOR MARGINAL DISTRIBUTIONS."
02700PRINT "IF YOU WANT TO CONTINUE YOUR ANALYSIS TYPE '1'."
02710PRINT "IF YOU DO NOT TYPE '0'.";
02720GOSUB 9000
02730I=O1
02740IF I=1 THEN 2850
02750IF I=0 THEN 2890
02760PRINT "REENTER. INPUT MUST BE 0 OR 1."
02770GOTO 2720
02780REM *********************************************************
02790REM
02800REM N1 DEGREES OF FREEDOM
02810REM L1 SCALE FOR INV CHI
02820REM F2 MEAN OF T DIST ON POPULKATION MEAN
02830REM LAST SCALE FOR T ON MEAN
02840REM
02850SCRATCH#2
02851 PRINT # 2,N1,L1,F2,L1*L1/N2,Q,0,N2
02860CHAIN "CMOD2A"
02870REM
02880REM *******************************************************
02890CHAIN "RSTRT"
02900REM ********************************************************
02910REM APPENDED GOSUBS
02920REM
02930REM *********************************************************
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*************************************************************
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