Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod66.bas
There are 2 other files named cmod66.bas in the archive. Click here to see a list.
00020REM *********************************************************************
00030REM CMOD66 CMOD66 CMOD66 CMOD66 CMOD66 CMOD66
00040REM *********************************************************************
00050REM
00060REM POSTERIOR ON SLOPE AND RESIDUAL ST DEV - LINEAR REGRESSION
00070REM
00080REM******************************************************************
00090 FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00130RESTORE#1
00131 INPUT# 1,I1,I2,I3
00140SCRATCH#1
00141 PRINT # 1,66,I2,I3
00150RESTORE#2
00151 INPUT# 2,X1,Y1,M8,X3,C3,M9,R
00160M9=M9+2
00170 DIM P(5),H(3),T(5)
00180P(2)=.25
00190P(3)=.5
00200P(1)=.1
00210P(4)=.75
00220P(5)=.9
00230H(1)=.5
00240H(2)=.75
00250H(3)=.95
00270GOSUB 6060
00280T(1)=.5
00290T(2)=.75
00300T(3)=.9
00310T(4)=.95
00320T(5)=.99
00330 E1$="REENTER. INPUT MUST BE POSITIVE."
00340PRINT L$
00350PRINT " POSTERIOR ANALYIS FOR SIMPLE LINEAR REGRESSION MODEL"
00360PRINT
00362IF M9>4 THEN 370
00363GOSUB 4700
00370PRINT "INPUT THE SAMPLE DATA."
00380PRINT
00390PRINT "SAMPLE SIZE =";
00400GOSUB 9000
00410IF O1>2 THEN 440
00420PRINT "REENTER. SAMPLE SIZE MUST BE GREATER THAN 2."
00430GOTO 400
00440N0=O1
00450PRINT
00460PRINT "MEAN OF PREDICTOR (X.) =";
00470GOSUB 9000
00480X0=O1
00490PRINT
00500PRINT "SUM OF SQUARED DEVIATIONS OF PREDICTOR (X-X.)(X-X.) =";
00510GOSUB 9000
00520IF O1>0 THEN 530
00530X2=O1
00540S4=(M8*X1+N0*X0)/(M8+N0)
00550PRINT
00560PRINT "MEAN OF CRITERION (Y.) =";
00570GOSUB 9000
00580Y0=O1
00590PRINT
00600PRINT "SUM OF SQUARED DEVIATIONS OF CRITERION (Y-Y.)(Y-Y.) =";
00610GOSUB 9000
00620IF O1>0 THEN 650
00630PRINT E1$
00640GOTO 610
00650Y2=O1
00660PRINT
00670PRINT "SUM OF CROSS PRODUCTS (X-X.)(Y-Y.) =";
00680GOSUB 9000
00730C2=O1
00740 S=Y2-C2/X2*C2
00750S0=R+S+X2*X3/(X2+X3)*(C3/X3-C2/X2)*(C3/X3-C2/X2)
00760IF X1=X0 THEN 800
00770 S1=M8*N0*(X1-X0)*(X1-X0)*(X3+X2)
00775 S1=S1/(M8+N0)/(X3+X2+M8*N0/(M8+N0)*(X1-X0)*(X1-X0))
00780S1=S1*((C3+C2)/(X3+X2)-(Y1-Y0)/(X1-X0)) ** 2
00790GOTO 820
00800 S1=0
00820S0=S0+S1
00830S3=X3+X2+(M8*N0)/(M8+N0)*(X1-X0)*(X1-X0)
00840S3=1/S3
00850S4=(M8*X1+N0*X0)/(M8+N0)
00860S5=1/(M8+N0)
00870S6=C3+C2+M8*N0/(M8+N0)*(X1-X0)*(Y1-Y0)
00880S6=S6/(X3+X2+M8*N0/(M8+N0)*(X1-X0)*(X1-X0))
00890S7=(M8*Y1+N0*Y0)/(M8+N0)
00900S8=M9+N0
00910S9=M9+N0-2
00920G=S9
00930PRINT L$
00940PRINT "HERE ARE SOME CHARACTERISTICS OF THE POSTERIOR "
00950PRINT "DISTRIBUTION ON THE RESIDUAL STANDARD DEVIATION."
00960PRINT
00970PRINT " INVERSE CHI"
00980PRINT
00990: DEGREES OF FREEDOM =#######.##
01000: SCALE PARAMETER =#######.###
01010: ##TH PERCENTILE =#######.###
01020: ######.### ##% HDR #######.###
01030PRINT USING 990,G
01040PRINT USING 1000,SQR(S0)
01050PRINT
01060FOR I=1 TO 5
01070P0=1-P(I)
01080GOSUB 1590
01090PRINT USING 1010,100*P(I),SQR(S0/X)
01100NEXT I
01110PRINT
01120FOR I=1 TO 3
01130J5=H(I)
01140GOSUB 7300
01150PRINT USING 1020,J1*SQR(S0),J5*100,J2*SQR(S0)
01160NEXT I
01170PRINT
01180PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01190GOSUB 9000
01200PRINT L$
01210PRINT "HERE ARE SOME CHARACTERISTICS OF THE POSTERIOR"
01220PRINT "DISTRIBUTION ON THE SLOPE OF THE REGRESSION LINE."
01230PRINT
01240PRINT " STUDENT'S T DISTRIBUTION"
01250PRINT
01260PRINT USING 990,G
01270PRINT USING 1000,S0*S3
01280: MEAN =#######.##
01290PRINT USING 1280,S6
01300PRINT
01310FOR I=1 TO 5
01320J5=T(I)
01330P0=J5/2+.5
01340GOSUB 2170
01350F=SQR(S0*S3/G)*J2
01360PRINT USING 1020,S6-F,J5*100,S6+F
01370NEXT I
01380PRINT
01390PRINT "IF YOU WANT TO LOOK AT EITHER POSTERIOR DISTRIBUTIONS ON Y"
01400PRINT "GIVEN X OR ON THE MEAN OF Y GIVEN X, TYPE '1', ELSE '0'.";
01410GOSUB 9000
01420IF O1=1 THEN 1570
01425CHAIN "RSTRT"
01430PRINT L$
01440PRINT "HERE ARE THE PARAMETERS OF THE REGRESSION LINE. IF"
01450PRINT "YOU LATER DECIDE THAT YOU WANT TO LOOK AT POSTERIOR"
01460PRINT "DISTRIBUTIONS THESE ARE THE VALUES YOU WILL BE ASKED"
01470PRINT "INPUT."
01480:DEGREES OF FREEDOM =#######.##
01490PRINT USING 1480,S9
01500:MEAN OF PREDICTOR =#######.##
01502PRINT USING 1500,S4
01510:MEAN OF CRITERION =#######.##
01512PRINT USING 1510,S7
01520:SLOPE =#######.##
01522PRINT USING 1520,S6
01530:TOTAL OBSERVATIONS =#######.##
01532PRINT USING 1530,1/S5
01540:SCALE FACTOR =#######.##
01542PRINT USING 1540,S0
01550:WEIGHTING FACTOR =#######.##
01552PRINT USING 1550,1/S3
01554PRINT
01556PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01557GOSUB 9000
01560CHAIN "RSTRT"
01570SCRATCH#3
01571 PRINT # 3,S9,S7,S6,S4,1/S5,S3,S0
01580CHAIN "CMOD67"
01590REM***********************************************************
01600REM CHI-SQUARE PERCENTILE
01610REM INPUTS P2
01620REM OUTPUTS J2
01630J=P0
01640P9=P0
01650P2=P0
01660GOSUB 1920
01670X=J2
01680GOSUB 5500
01690IF ABS(P-P9)<.0001 THEN 1910
01700P4=P
01710E7=J2
01720P2=P2-.009
01730GOSUB 1920
01740E1=J2
01750E2=E7
01760X=E1
01770GOSUB 5500
01780P3=P
01790X3=(P9-P3)/(P4-P3)
01800X3=X3*(E2-E1)+E1
01810X=X3
01820GOSUB 5500
01830IF ABS(P9-P)<.0001 THEN 1910
01840IF P>P9 THEN 1880
01850E1=X3
01860P3=P
01870GOTO 1790
01880E2=X3
01890P4=P
01900GOTO 1790
01910RETURN
01920P3=P2
01930IF P2 <= .5 THEN 1950
01940P2=1-P2
01950A1=SQR(LOG(1/P2/P2))
01960A2=2.51552+.802853*A1+.010328*A1*A1
01970A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
01980J2=A1-A2
01990IF P3 <= .5 THEN 2020
02000P2=P3
02010GOTO 2030
02020J2=-J2
02030J2=G*(1-2/9/G+J2*SQR(2/9/G)) ** 3
02040RETURN
02050REM
02060REM********************* END OF CSQ PERCENTILE ********************
02070REM ********************************************************
02080REM SUBROUTINE FOR COMPUTING
02090REM DENSITY OF INVERSE CHI
02100REM ********************************************************
02110G9=.5*G
02120GOSUB 5860
02130F0=LOG(G0)
02140F0=LOG(2)-F0+.5*G*LOG(.5*S0*S0)-(G+1)*LOG(J2)-.5*S0*S0/J2 ** 2
02150D2=F0
02160RETURN
02170REM********************** PERCENTILE FINDER **************************
02180P5=2
02190P6=2
02200P2=P0
02210GOSUB 2470
02220IF ABS(P-P0)<.0001 THEN 2460
02230IF P>P0 THEN 2300
02240E5=J2
02250P5=P
02260IF P6 <> 2 THEN 2360
02270P2=P2+.001
02280GOSUB 2470
02290GOTO 2220
02300E6=J2
02310P6=P
02320IF P5 <> 2 THEN 2360
02330P2=P2-.001
02340GOSUB 2470
02350GOTO 2220
02360J2=(P0-P5)/(P6-P5)*(E6-E5)+E5
02370GOSUB 6000
02380IF ABS(P-P0)<.0001 THEN 2460
02390IF P>P0 THEN 2430
02400P5=P
02410E5=J2
02420GOTO 2360
02430E6=J2
02440P6=P
02450GOTO 2360
02460RETURN
02470P3=P2
02480IF P2 <= .5 THEN 2500
02490P2=1-P2
02500A1=SQR(LOG(1/P2/P2))
02510A2=2.51552+.802853*A1+.010328*A1*A1
02520A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
02530A2=A1-A2
02540J2=SQR(G*EXP(A2*(G-5/6)*A2/(G-2/3+.1/G)/(G-2/3+.1/G))-G)
02550GOSUB 6000
02560IF P3 <= .5 THEN 2590
02570P2=P3
02580GOTO 2610
02590J2=-J2
02600P=1-P
02610RETURN
02620REM ****************************************************
02630REM LOG GAMMA ROUTINE
02640REM INPUT G9
02650REM OUTPUT G0
02660G5=G9
02670IF G9 <= 1.E+30 THEN 2700
02680G0=1.E+38
02690RETURN
02700IF G9>1.E-09 THEN 2730
02710G0=0
02720RETURN
02730IF G9<1.E+10 THEN 2760
02740G0=G9*(LOG(G9)-1)
02750RETURN
02760G6=1
02770IF 18<G5 THEN 2810
02780G6=G6*G5
02790G5=G5+1
02800GOTO 2770
02810R8=1/G5/G5
02820G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
02830C1=8.33333E-02
02840C2=2.77778E-03
02850C3=7.93651E-04
02860C4=5.95238E-04
02870G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
02880RETURN
02890REM END OF LOG GAMMA ROUTINE
02900REM ****************************************************
04700PRINT
04710PRINT "INPUT THE PRIOR INFORMATION ABOUT THE REGRESSION LINE."
04720PRINT
04730PRINT "DEGREES OF FREEDOM=";
04740GOSUB 9000
04750IF O1 >= 4 THEN 4780
04760PRINT "REENTER. DEGREES OF FREEDOM MUST BE AT LEAST 4."
04770GOTO 4740
04780M9=O1+2
04790PRINT
04800PRINT "MEAN OF PREDICTOR=";
04810GOSUB 9000
04820X1=O1
04830PRINT
04840PRINT "MEAN OF CRITERION=";
04850GOSUB 9000
04860Y1=O1
04870PRINT
04880PRINT "SLOPE=";
04890GOSUB 9000
04900W6=O1
04910PRINT
04920PRINT "OBSERVATIONS ON THE MEAN=";
04930GOSUB 9000
04940IF O1>0 THEN 4970
04950PRINT "REENTER. MUST BE LARGER THAN 0."
04960GOTO 4930
04970M8=O1
04980PRINT
04990PRINT "SCALE FACTOR=";
05000GOSUB 9000
05010IF O1>0 THEN 5040
05020PRINT "REENTER. MUST BE POSITIVE."
05030GOTO 5000
05040R=O1
05050PRINT
05060PRINT "WEIGHTING FACTOR=";
05070GOSUB 9000
05080IF O1>0 THEN 5110
05090PRINT "REENTER. MUST BE POSITIVE."
05100GOTO 5070
05110X3=O1
05120C3=X3*W6
05122PRINT L$
05130RETURN
05500REM ********************************************************
05510REM CHI-SQUARE CDF ROUTINE-LOWER TAIL
05520REM INPUT G X
05530REM OUTPUT P
05540REM*************************
05550REM PEIZER-PRATT-APPROXIMATION
05560REM PROB(CHISQ<X) WITH DF=G
05570REM*************************
05580X2=X/2
05590T=X2-G*.5+1/3-.04/G
05600Y=(G*.5-.5)/X2
05610IF Y<.00001 THEN 5680
05620IF ABS(Y-1)<.00001 THEN 5700
05630G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05640IF G6>-1 THEN 5710
05650Y3=T
05660GOTO 5720
05670GOTO 5710
05680G6=1
05690GOTO 5710
05700G6=0
05710Y3=T*SQR((1+G6)/X2)
05720GOSUB 8000
05730P2=P
05740RETURN
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/G5
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
06007REM PRIOR GOSUB 6060
06008P=0
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=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
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**************************************************
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