Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod64.bas
There are 2 other files named cmod64.bas in the archive. Click here to see a list.
00020REM***************************************************************
00030REM CMOD64 CMOD64 CMOD64 CMOD64 CMOD64 CMOD64
00040REM******************************************************************
00050FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00090RESTORE#1
00091 INPUT# 1,I1,I2,I3
00100SCRATCH#1
00101 PRINT # 1,64,I2,I3
00105X=0
00110GOSUB 6060
00120Z9=0
00122Z8=0
00123Z7=3
00124V9=0
00130 DIM H(5)
00140H(1)=.5
00150H(2)=.75
00160H(3)=.9
00170N=5
00190H(4)=.95
00200 E1$="REENTER. INPUT MUST BE NUMBER OF OPTION."
00210 C1$="WHEN YOU ARE READY TO CONTINUE TYPE '1'."
00220H(5)=.99
00230PRINT L$
00240PRINT " ADVERSARY PREPOSTERIOR ANALYSIS FOR TW0-PARAMETER NORMAL"
00250PRINT
00260PRINT "THE PURPOSE OF AN ADVERSARY PREPOSTERIOR ANALYSIS IS TO GIVE"
00270PRINT "YOU THE OPPORTUNITY TO SEE WHAT YOUR PRIOR BELIEFS IMPLY YOU"
00280PRINT "EXPECT AN ADVERSARY (SOMEONE WITH DIFFERENT PRIOR BELIEFS) "
00290PRINT "TO BELIEVE AFTER SOME SAMPLE OBSERVATIONS ARE MADE. YOU, OF"
00300PRINT "COURSE, EXPECT THE SAMPLE DATA TO BE CONSISTENT WITH YOUR PRIOR"
00310PRINT "BELIEFS."
00320PRINT
00330PRINT "FIRST, THE MODULE ALLOWS YOU T0 SEE THE EFFECT OF DIFFERENT SAMPLE"
00340PRINT "SIZES ON THE MEAN AND STANDARD DEVIATION OF THE PREPOSTERIOR"
00350PRINT "DISTRIBUTION. THIS SHOULD PROVIDE YOU WITH A ROUGH IDEA OF THE"
00360PRINT "EXPECTED EFFECT OF DIFFERENT SAMPLE SIZES. YOU CAN THEN LOOK"
00370PRINT "MORE CLOSELY AT THE PREPOSTERIOR DISTRIBUTION FOR DIFFERENT"
00380PRINT "SAMPLE SIZES."
00390PRINT
00400 PRINT C1$;
00410GOSUB 9000
00420PRINT L$
00430PRINT "THE MODULE ALLOWS YOU TO CARRY OUT AN ANALYSIS ON THE MEAN"
00440PRINT "OR ON THE N+1ST OBSERVATION. IN OTHER WORDS, YOU CAN SEE"
00450PRINT "WHAT YOU EXPECT ADVERSARY TO BELIEVE ABOUT THE MEAN OR THE "
00460PRINT "NEXT OBSERVATION AFTER HE HAS ALREADY MADE N OBSERVATIONS."
00470PRINT
00480Z9=0
00490PRINT "PREPOSTERIOR ON: MEAN=1 N+1ST OBSERVATION=2 OR EXIT=3";
00500GOSUB 9000
00510IF O1=1 THEN 590
00520IF O1=2 THEN 580
00530IF O1=3 THEN 560
00540PRINT E1$
00550GOTO 500
00560PRINT L$
00570CHAIN "CMOD2A"
00580Z9=1
00590PRINT L$
00600IF Z8=0 THEN 650
00610PRINT
00620PRINT "CHANGE PRIORS: YOURS=1 ADVERSARY-2 BOTH=3 NEITHER=4";
00630GOSUB 9000
00640O1=INT(O1)
00642 ONO1 GOTO 646,646,646,646
00643PRINT E1$
00644GOTO 630
00646Z7=O1
00650Z8=1
00652 ONZ7 GOTO 660,810,660,1180
00660PRINT "INPUT THE PARAMETERS OF YOUR PRIOR MARGINAL DISTRIBUTION ON THE"
00670PRINT "MEAN."
00680GOSUB 860
00690V3=V0
00700L3=S0
00710M3=M0
00720PRINT
00730PRINT "INPUT THE SCALE PARAMETER OF YOUR PRIOR MARGINAL DISTRIBUTION"
00740PRINT "ON THE STANDARD DEVIATION.";
00750GOSUB 9000
00760IF O1>0 THEN 790
00770PRINT "REENTER. SCALE PARAMETER MUST BE GREATER THAN 0."
00780GOTO 750
00790H3=O1*O1
00800Q3=H3/L3
00805 ONZ7 GOTO 1180,810,810,1180
00810PRINT
00820PRINT "INPUT THE PARAMETERS OF THE ADVERSARY MARGINAL PRIOR"
00830PRINT "DISTRIBUTION ON THE MEAN."
00840GOSUB 860
00850GOTO 1040
00860PRINT
00870PRINT "INPUT THE DEGREES OF FREEDOM.";
00880GOSUB 9000
00890IF O1>2 THEN 920
00900PRINT "REENTER. DEGREES OF FREEDOM MUST BE GREATER THAN 2."
00910GOTO 880
00920V0=O1
00930PRINT "INPUT THE MEAN.";
00940GOSUB 9000
00950M0=O1
00960PRINT "INPUT THE SCALE PARAMETER.";
00970GOSUB 9000
00980S0=O1
00990IF S0>0 THEN 1020
01000PRINT "REENTER. SCALE PARAMETER MUST BE GREATER THAN 0."
01010GOTO 970
01020S0=O1
01030RETURN
01040V1=V0
01050M1=M0
01060L1=S0
01070S1=L1/(V1-2)
01080PRINT
01090PRINT "INPUT THE SCALE PARAMETER OF THE PRIOR MARGINAL DISTRIOBUTION"
01100PRINT "ON THE STANDARD DEVIATION.";
01110GOSUB 9000
01120IF O1>0 THEN 1150
01130PRINT "REENTER. SCALE PARAMETER MUST BE GREATER THAN 0."
01140GOTO 1110
01150H1=O1*O1
01160Q1=H1/L1
01180PRINT L$
01190IF V9=1 THEN 1250
01200V9=1
01210PRINT "THE MODULE WILL DISPLAY THE MEAN AND STANDARD DEVIATION OF"
01220PRINT "THE PREPOSTERIOR DISTRIBUTION FOR A GIVEN SAMPLE SIZE."
01230PRINT "YOU CAN SPECIFY AS MANY AS 5 DIFFERENT SAMPLE SIZES."
01240PRINT
01250PRINT "INPUT THE NUMBER OF SAMPLE SIZES YOU WANT TO CONSIDER.";
01260GOSUB 9000
01270IF O1 <= 5 THEN 1300
01280PRINT "REENTER. YOU CAN ONLY CONSIDER 5 AT ONE TIME."
01290GOTO 1260
01300IF O1<1 THEN 1280
01310N0=O1
01320FOR K=1 TO N0
01330PRINT "INPUT SAMPLE SIZE. (MIN=3)";
01340GOSUB 9000
01350IF O1 >= 3 THEN 1380
01360PRINT "REENTER. SAMPLE SIZE MUST BE AT LEAST 3."
01370GOTO 1340
01380N(K)=O1
01390NEXT K
01400PRINT L$
01410PRINT "HERE ARE THE MEANS AND STANDARD DEVIATIONS OF THE PREPOSTERIOR"
01420PRINT "DISTRIBUTIONS FOR DIFFERENT SAMPLE SIZES."
01430PRINT
01440IF Z9=1 THEN 1470
01450PRINT "PREPOSTERIOR FOR MEAN"
01460GOTO 1480
01470PRINT "PREPOSTERIOR FOR N+1TH OBSERVATION"
01480PRINT
01490PRINT " MEAN STANDARD DEVIATION"
01500:YOUR PRIOR ######.## #######.##
01510PRINT USING 1500,M3,SQR(H3*(1/Q3+Z9)/(V3-2))
01520:ADVERSARY PRIOR ######.## #######.##
01530PRINT USING 1520,M1,SQR(H1*(1/Q1+Z9)/(V1-2))
01540PRINT "----------------------------------------------------------------"
01550PRINT "PREPOSTERIOR"
01560FOR K=1 TO N0
01570N=N(K)
01580GOSUB 2270
01590: N =##### ######.## #######.##
01600PRINT USING 1590,N,M0,SQR(S0/(G-2))
01610NEXT K
01620PRINT
01630PRINT "IF YOU WANT TO CONSIDER OTHER N VALUES TYPE '1', ELSE '0'.";
01640GOSUB 9000
01650IF O1=1 THEN 1670
01660GOTO 1690
01670PRINT L$
01680GOTO 1250
01690PRINT L$
01700PRINT "YOU CAN NOW LOOK IN MORE DETAIL AT THE PREPOSTERIOR"
01710PRINT "DISTRIBUTION FOR ANY N YOU WANT."
01720PRINT
01730PRINT "IF YOU WANT TO DO THIS TYPE THE VALUE OF N, ELSE '0'.";
01740GOSUB 9000
01750IF O1 <> 0 THEN 1800
01760PRINT L$
01770PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
01780PRINT
01790GOTO 480
01800IF O1 >= 3 THEN 1830
01810PRINT "REENTER. MINIMUM N IS 3."
01820GOTO 1740
01830PRINT L$
01840PRINT "THE PREPOSTERIOR DISTRIBUTION CAN BE APPROXIMATED"
01850PRINT "BY A T-DISTRIBUTION WITH THESE PARAMETERS."
01860N=O1
01870GOSUB 2270
01880PRINT
01890PRINT
01900IF Z9=1 THEN 1930
01910PRINT "PREPOSTERIOR ON MEAN"
01920GOTO 1940
01930PRINT "PREPOSTERIOR ON THE N+1TH OBSERVATION"
01940PRINT
01950:DEGREES OF FREEDOM =####.## MEAN =######.##
01960PRINT USING 1950,G,M0
01970:SCALE PARAMETER=#######.##
01980PRINT USING 1970,S0
01990PRINT
02000PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THIS DISTRIBUTION."
02010 PRINT
02020: ########.## ## % HDR ########.##
02030 F=SQR(S0/G)
02040 FOR I=1 TO 5
02050 J5=H(I)/2+.5
02060 P0=J5
02070 GOSUB 5000
02080 PRINT USING 2020,M0-J2*F,H(I)*100,M0+J2*F
02090 NEXT I
02100PRINT
02110PRINT"IF YOU WANT THE PROBABILITY LESS THAN SOME VALUE TYPE THE"
02120PRINT"VALUE (EXIT=-7777).";
02130 GOSUB 9000
02140 IF O1<>-7777 THEN 2180
02150 PRINT L$
02160PRINT"IF YOU WANT TO TRY A DIFFERENT N TYPE THE VALUE, ELSE '0'.";
02170 GOTO 1740
02180 J1=0
02190 J2=ABS(O1-M0)*SQR(G/S0)
02200 GOSUB 6000
02210 IF O1>M0 THEN 2240
02220 P=1-P
02230: PROB LESS THAN ########.## = ##.##
02240 PRINT USING 2230 ,O1,P
02250 PRINT"NEXT VALUE OR '-7777'";
02260 GOTO 2130
02270 X5=H1+(N-1)*H3/(V3-2)
02280X5=X5+(Q1*N)/(Q1+N)*(H3*(N+Q3)/(V3-2)/(Q3*N)+M3*M3-2*M3*M1+M1*M1)
02290 X5=X5/(Q1+N)
02300 X5=X5/(V1+N-2)
02310 IF Z9=0 THEN 2330
02320 X5=X5*(Q1+N-1)
02330 X5=X5+N*N*H3/(Q1+N)/Q3/(V3-2)/(Q1+N)/N*(N+Q3)
02340 M0=(Q1*M1+N*M3)/(N+Q1)
02350 G=(N*V3+V1)/(N+1)
02360 S0=(G-2)*X5
02370RETURN
05000REM********************** PERCENTILE FINDER **************************
05010P5=2
05020P6=2
05030P2=P0
05040GOSUB 5300
05050IF ABS(P-P0)<.0001 THEN 5290
05060IF P>P0 THEN 5130
05070E5=J2
05080P5=P
05090IF P6 <> 2 THEN 5190
05100P2=P2+.001
05110GOSUB 5300
05120GOTO 5050
05130E6=J2
05140P6=P
05150IF P5 <> 2 THEN 5190
05160P2=P2-.001
05170GOSUB 5300
05180GOTO 5050
05190 J2=(P0-P5)/(P6-P5)*(E6-E5)+E5
05200GOSUB 6000
05210IF ABS(P-P0)<.0001 THEN 5290
05220IF P>P0 THEN 5260
05230P5=P
05240E5=J2
05250GOTO 5190
05260E6=J2
05270P6=P
05280GOTO 5190
05290RETURN
05300P3=P2
05310IF P2 <= .5 THEN 5330
05320P2=1-P2
05330A1=SQR(LOG(1/P2/P2))
05340A2=2.51552+.802853*A1+.010328*A1*A1
05350A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
05360A2=A1-A2
05370J2=SQR(G*EXP(A2*(G-5/6)*A2/(G-2/3+.1/G)/(G-2/3+.1/G))-G)
05380GOSUB 6000
05390IF P3 <= .5 THEN 5420
05400P2=P3
05410GOTO 5440
05420J2=-J2
05430P=1-P
05440RETURN
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
06008P=C0
06009IF J2<.0001 THEN 6056
06010IF G<6 THEN 6015
06011IF J2<6 THEN 6016
06012P=1
06014GOTO 6058
06015IF J2>12 THEN 6012
06016 DIM W(16),O(16)
06018Y3=J2
06020IF G<10 THEN 6026
06021Y3=(G-2/3+.1/G)*SQR(ABS(LOG(1+J2*J2/G))/(G-5/6))
06022GOSUB 8000
06024GOTO 6058
06025REM PEIZER PRATT APPROXIMATION
06026IF G=1 THEN 6098
06028J1=C0
06030N=G
06032P=C0
06034GOSUB 6084
06036D0=(J2-J1)*.5
06038D1=(J1+J2)*.5
06040FOR I1=1 TO 16
06042D9=D0*O(I1)+D1
06044IF D9=C0 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=C0
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=C0
07521J2=X9
07522GOSUB 6000
07523P=2*P-1
07524IF ABS(X2-X9)<.0001 THEN 7541
07525IF P<J5 THEN 7538
07526X2=X9
07527Z9=P
07528X9=(J5-Z8)/(Z9-Z8)
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
08000REM
08005Y4=ABS(Y3)
08010X1=X
08015X=Y3
08020T=1/(1+.231642*Y4)
08021IF X*X/2<80 THEN 8025
08022D=C0
08023GOTO 8030
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
09000REM
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09999 END