Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod67.bas
There are 2 other files named cmod67.bas in the archive. Click here to see a list.
00020REM *********************************************************************
00030REM CMOD67 CMOD67 CMOD67 CMOD67 CMOD67 CMOD67 CMOD67
00040REM *********************************************************************
00050REM
00060REM POSTERIOR PREDICTIVE DISTRIBUTION - SIMPLE LINEAR REGRESSION
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,67,I2,I3
00150RESTORE#3
00151 INPUT# 3,G,W5,W6,W7,W8,W9,W0
00152IF G>0 THEN 160
00154PRINT "YOU HAVE NOT YET SPECIFIED A PRIOR DISTRIBUTION. YOU MUST DO THIS"
00156PRINT "BEFORE YOU CAN USE THIS MODULE. TYPE '1' TO CONTINUE."
00158GOSUB 9000
00159CHAIN "RSTRT"
00160 DIM T(10),Y(10),X(10)
00170X=0
00180MAT T=ZER
00190 DIM V(6)
00210MAT X=ZER
00220MAT Y=ZER
00230GOSUB 6060
00240 E1$="REENTER. INPUT MUST BE NUMBER OF OPTION."
00250 S1$="------------------------------------------------------------"
00260PRINT L$
00270PRINT " EVALUATION OF POSTERIOR DISTRIBUTIONS ON Y AND THE MEAN OF Y"
00280PRINT
00290PRINT "THIS MODULE ALLOWS YOU TO EXAMINE THE POSTERIOR DISTRIBUTIONS ON"
00300PRINT "Y AND THE MEAN OF Y GIVEN X."
00305PRINT
00310IF G>0 THEN 680
00320PRINT
00330PRINT "INPUT THE PARAMETERS OF THE REGRESSION LINE."
00340PRINT
00350PRINT "DEGREES OF FREEDOM=";
00360GOSUB 9000
00370IF O1>6 THEN 400
00380PRINT "REENTER. DEGREES OF FREEDOM MUST BE GREATER THAN 6."
00390GOTO 360
00400G=O1
00410PRINT
00420PRINT "MEAN OF PREDICTOR=";
00430GOSUB 9000
00440W7=O1
00450PRINT
00460PRINT "MEAN OF CRITERION=";
00470GOSUB 9000
00480W5=O1
00490PRINT
00500PRINT "SLOPE=";
00510GOSUB 9000
00520W6=O1
00530PRINT
00540PRINT "TOTAL OBSERVATIONS=";
00550GOSUB 9000
00560IF O1>G THEN 590
00570PRINT "REENTER. MUST BE LARGER THAN THE DEGREES OF FREEDOM= ";G
00580GOTO 550
00590W8=O1
00600PRINT
00610PRINT "SCALE FACTOR=";
00620GOSUB 9000
00630IF O1>0 THEN 660
00640PRINT "REENTER. MUST BE POSITIVE."
00650GOTO 620
00660W0=O1
00670PRINT
00672PRINT "WEIGHTING FACTOR=";
00673GOSUB 9000
00674IF O1>0 THEN 678
00675PRINT "REENTER. MUST BE POSITIVE."
00676GOTO 673
00678W9=1/O1
00680PRINT "DISTRIBUTION ON ( Y=1 MEAN OF Y=2 ) OR EXIT=0 ";
00690GOSUB 9000
00700IF O1=0 THEN 900
00710IF O1=1 THEN 750
00720IF O1=2 THEN 750
00730PRINT E1$
00740GOTO 690
00750B3=O1
00760PRINT
00770PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
00780PRINT " 1. SPECIFIED PERCENTILES FOR SPECIFIED X VALUES"
00790PRINT " 2. SPECIFIED HDR'S FOR SPECIFIED X VALUES"
00800PRINT " 3. PROBABILITY LESS THAN SPECIFIED VALUES FOR"
00810PRINT " SPECIFIED X VALUES."
00820PRINT " 4. CHANGE TYPE OF DISTRIBUTION (E.G. FROM Y TO MEAN OF Y)"
00822PRINT " 5. EXIT MODULE"
00830GOSUB 9000
00840IF O1=1 THEN 2210
00850IF O1=2 THEN 910
00860IF O1=4 THEN 680
00870IF O1=3 THEN 1540
00872IF O1 <> 5 THEN 880
00873CHAIN "RSTRT"
00880PRINT E1$
00890GOTO 830
00900CHAIN "RSTRT"
00910PRINT L$
00920PRINT " OPTION 2: P% HDR'S FOR SPECIFIED X VALUES"
00930PRINT
00940PRINT "YOU CAN SPECIFY 1 OR 2 DIFFERENT P% HDR'S AND FROM 1"
00950PRINT "THRU 10 DIFFERENT VALUES."
00960PRINT
00970GOSUB 1010
00980GOTO 1150
00990REM*************************************************************
01000REM INPUT THE P% HDR ' S
01010PRINT "INPUT 1 OR 2 FOR THE NUMBER OF P% HDR'S.";
01020GOSUB 9000
01030IF O1=1 THEN 1070
01040IF O1=2 THEN 1070
01050PRINT "REENTER. INPUT MUST BE 1 OR 2."
01060GOTO 1020
01070B1=O1
01080FOR I=1 TO B1
01090PRINT "INPUT P%.";
01100GOSUB 9000
01110Y(I)=O1
01120NEXT I
01130RETURN
01140REM**********************************************************************
01150PRINT
01160GOSUB 2710
01170PRINT L$
01180PRINT " P% HIGHEST DENSITY REGIONS"
01190PRINT
01200GOSUB 3720
01210 ONB1 GOTO 1260,1230
01220: X P%=##.# P%=##.#
01230PRINT USING 1220,Y(1),Y(2)
01240GOTO 1270
01250: X P%=##.#
01260PRINT USING 1250,Y(1)
01270PRINT S1$
01280FOR I=1 TO B0
01290GOSUB 3680
01300FOR J=1 TO B1
01310P0=Y(J)/200+.5
01320GOSUB 3810
01330T(J*2)=J2*SQR(S0/G)+M0
01340T(J*2-1)=M0-J2*SQR(S0/G)
01350NEXT J
01360GOSUB 3400
01370NEXT I
01380REM*****************************************************************
01390PRINT S1$
01400PRINT "CHANGE( P% VALUES=1 X VALUES=2 BOTH=3) OR EXIT=0 ";
01410GOSUB 9000
01420O1=O1+1
01430 ONO1 GOTO 1460,1460,1460,1460
01440PRINT E1$
01450GOTO 1410
01460PRINT L$
01470 ONO1 GOTO 770,1480,1160,970
01480GOSUB 1010
01490GOTO 1170
01500PRINT L$
01510GOTO 1600
01520REM*****************************************************************
01530REM
01540PRINT L$
01550PRINT " OPTION 3: PROBABILITY LESS THAN Y0"
01560PRINT
01570PRINT "YOU CAN SPECIFY FROM 1 THRU 4 DIFFERENT Y0 VALUES AND"
01580PRINT "FROM 1 THRU 10 DIFFERENT X VALUES."
01590PRINT
01600GOSUB 1640
01610GOTO 1760
01620REM******************************
01630REM INPUT Y0 VALUES
01640PRINT "INPUT THE NUMBER OF DIFFERENT Y0 VALUES.";
01650GOSUB 9000
01660IF O1<1 THEN 2180
01670IF O1>4 THEN 2180
01680B1=O1
01690FOR I=1 TO B1
01700PRINT "INPUT Y0.";
01710GOSUB 9000
01720Y(I)=O1
01730NEXT I
01740RETURN
01750REM******************************************************************
01760PRINT
01770GOSUB 2710
01780PRINT L$
01790PRINT " OPTION 3: PROBABILITY LESS THAN Y0."
01800PRINT
01810GOSUB 3720
01820 ONB1 GOTO 1840,1870,1900,1930
01830: X #######.##
01840PRINT USING 1830,Y(1)
01850GOTO 1940
01860: X #######.## ######.##
01870PRINT USING 1860,Y(1),Y(2)
01880GOTO 1940
01890: X #######.## ######.## ######.##
01900PRINT USING 1890,Y(1),Y(2),Y(3)
01910GOTO 1940
01920: X #######.## ######.## ######.## ######.##
01930PRINT USING 1920,Y(1),Y(2),Y(3),Y(4)
01940PRINT S1$
01950FOR I=1 TO B0
01960GOSUB 3680
01970FOR J=1 TO B1
01980J2=ABS((Y(J)-M0)*SQR(G/S0))
01990GOSUB 6000
02000IF Y(J) >= M0 THEN 2030
02010T(J)=1-P
02020GOTO 2040
02030T(J)=P
02040NEXT J
02050GOSUB 3490
02060NEXT I
02070PRINT S1$
02080PRINT "CHANGE ( Y0 VALUES=1 X VALUES=2 BOTH=3) OR EXIT=0 ";
02090GOSUB 9000
02100O1=O1+1
02110 ONO1 GOTO 2140,2140,2140,2140
02120PRINT E1$
02130GOTO 2090
02140PRINT L$
02150 ONO1 GOTO 770,2160,1770,1600
02160GOSUB 1640
02170GOTO 1780
02180PRINT "REENTER. MINIMUM IS 1 AND MAXIMUM IS 4."
02190GOTO 1650
02200REM*******************************************************************
02210PRINT L$
02220PRINT " OPTION 1: PERCENTILES FOR SPECIFIED X VALUES"
02230PRINT
02240PRINT "YOU CAN SPECIFY FROM 1 THRU 4 DIFFERENT PERCENTILES"
02250PRINT "AND FROM 1 THRU 10 DIFFERENT X VALUES."
02260PRINT
02270GOSUB 2310
02280GOTO 2570
02290REM*************************************************************
02300REM INPUT PERCENTILES
02310PRINT "INPUT NUMBER OF DIFFERENT PERCENTILES.";
02320GOSUB 9000
02330IF O1<1 THEN 2350
02340IF O1 <= 4 THEN 2370
02350PRINT "REENTER. MINIMUM IS 1 AND MAXIMUM IS 4."
02360GOTO 2320
02370B1=O1
02380FOR I=1 TO B1
02390PRINT "INPUT PERCENTILE.";
02400GOSUB 9000
02410IF O1<1 THEN 3190
02420IF O1>99 THEN 3190
02430P(I)=O1
02440IF O1=50 THEN 2530
02450IF O1>50 THEN 2500
02460P0=1-O1/100
02470GOSUB 3810
02480J2=-J2
02490GOTO 2540
02500P0=O1/100
02510GOSUB 3810
02520GOTO 2540
02530J2=0
02540V(I)=J2
02550NEXT I
02560RETURN
02570PRINT
02580GOSUB 2710
02590GOTO 2820
02600REM**************************************************************
02610REM INPUT THE X VALUES
02620PRINT L$
02630IF B4=0 THEN 2700
02640PRINT "IF YOU WANT TO CHANGE THE X VALUES TYPE '1', ELSE '0'.";
02650GOSUB 9000
02660IF O1=0 THEN 2810
02670IF O1=1 THEN 2710
02680PRINT E1$
02690GOTO 2650
02700B4=1
02710PRINT "INPUT THE NUMBER OF DIFFERENT X VALUES.";
02720GOSUB 9000
02730IF O1>10 THEN 3210
02740IF O1<1 THEN 3210
02750B0=O1
02760FOR I=1 TO B0
02770PRINT "INPUT X.";
02780GOSUB 9000
02790X(I)=O1
02800NEXT I
02810RETURN
02820PRINT L$
02830PRINT " OPTION 1: PERCENTILES"
02840PRINT
02850GOSUB 3720
02860 ONB1 GOTO 2880,2910,2940,2970
02870: X #######.##
02880PRINT USING 2870,P(1)
02890GOTO 3000
02900: X #######.## ######.##
02910PRINT USING 2900,P(1),P(2)
02920GOTO 3000
02930: X #######.## ######.## ######.##
02940PRINT USING 2930,P(1),P(2),P(3)
02950GOTO 3000
02960: X #######.## ######.## ######.## ######.##
02970PRINT USING 2960,P(1),P(2),P(3),P(4)
02980GOTO 3000
02990: X #######.## ######.## ######.## ######.##
03000PRINT S1$
03010FOR I=1 TO B0
03020GOSUB 3680
03030FOR J=1 TO B1
03040T(J)=V(J)*SQR(S0/G)+M0
03050NEXT J
03060GOSUB 3250
03070NEXT I
03080PRINT S1$
03090PRINT "CHANGE (PERCENTILES=1 X VALUES=2 BOTH=3) OR EXIT=0 ";
03100GOSUB 9000
03110O1=O1+1
03120 ONO1 GOTO 3150,3150,3150,3150
03130PRINT E1$
03140GOTO 3100
03150PRINT L$
03160 ONO1 GOTO 770,3170,2580,2270
03170GOSUB 2310
03180GOTO 2820
03190PRINT "REENTER. MINIMUM PERCENTILES IS 1 AND MAXIMUM IS 99."
03200GOTO 2400
03210PRINT "REENTER. MINIMUM IS 1 AND MAXIMUM IS 10."
03220GOTO 2720
03230REM****************************************************************
03240REM PERCENTILES
03250 ONB1 GOTO 3270,3300,3330,3360
03260:######.## ######.##
03270PRINT USING 3260,X(I),T(1)
03280GOTO 3370
03290:######.## ######.## ######.##
03300PRINT USING 3290,X(I),T(1),T(2)
03310GOTO 3370
03320:######.## ######.## ######.## ######.##
03330PRINT USING 3320,X(I),T(1),T(2),T(3)
03340GOTO 3370
03350:######.## ######.## ######.## ######.## ######.##
03360PRINT USING 3350,X(I),T(1),T(2),T(3),T(4)
03370RETURN
03380REM********************************************************************
03390REM P% HDR ' S
03400 ONB1 GOTO 3420,3450
03410:#######.## *#######.## #######.##
03420PRINT USING 3410,X(I),T(1),T(2)
03430RETURN
03440:#######.## *#######.## #######.## *#######.## #######.##
03450PRINT USING 3440,X(I),T(1),T(2),T(3),T(4)
03460RETURN
03470REM**************************************************************
03480REM PROB LESS THAN Y0
03490 ONB1 GOTO 3510,3540,3570,3600,3630
03500:#######.## ######.##
03510PRINT USING 3500,X(I),T(1)
03520GOTO 3640
03530:#######.## ######.## ######.##
03540PRINT USING 3530,X(I),T(1),T(2)
03550GOTO 3640
03560:#######.## ######.## ######.## ######.##
03570PRINT USING 3560,X(I),T(1),T(2),T(3)
03580GOTO 3640
03590:#######.## ######.## ######.## ######.## ######.##
03600PRINT USING 3590,X(I),T(1),T(2),T(3),T(4)
03610GOTO 3640
03620:######.## ###.## ###.## ###.## ###.## ###.##
03630PRINT USING 3620,X(I),T(1),T(2),T(3),T(4),T(5)
03640RETURN
03650REM STAN T
03660J2=(X-M0)*G/S0
03670RETURN
03680M0=W5+W6*(X(I)-W7)
03690IF B3=1 THEN 3790
03700S0=W0*(1/W8+(X(I)-W7)*(X(I)-W7)*W9)
03710RETURN
03720 ONB3 GOTO 3730,3760
03730PRINT "POSTERIOR PREDICTIVE DISTRIBUTION ON Y GIVEN X"
03740PRINT
03750RETURN
03760PRINT "POSTERIOR DISTRIBUTION ON MEAN OF Y GIVEN X."
03770PRINT
03780RETURN
03790S0=W0*(1+1/W8+(X(I)-W7)*(X(I)-W7)*W9)
03800RETURN
03810REM********************** PERCENTILE FINDER **************************
03820P5=2
03830P6=2
03840P2=P0
03850GOSUB 4110
03860IF ABS(P-P0)<.0001 THEN 4100
03870IF P>P0 THEN 3940
03880E5=J2
03890P5=P
03900IF P6 <> 2 THEN 4000
03910P2=P2+.001
03920GOSUB 4110
03930GOTO 3860
03940E6=J2
03950P6=P
03960IF P5 <> 2 THEN 4000
03970P2=P2-.001
03980GOSUB 4110
03990GOTO 3860
04000J2=.5*(E6+E5)
04010GOSUB 6000
04020IF ABS(P-P0)<.0001 THEN 4100
04030IF P>P0 THEN 4070
04040P5=P
04050E5=J2
04060GOTO 4000
04070E6=J2
04080P6=P
04090GOTO 4000
04100RETURN
04110P3=P2
04120IF P2 <= .5 THEN 4140
04130P2=1-P2
04140A1=SQR(LOG(1/P2/P2))
04150A2=2.51552+.802853*A1+.010328*A1*A1
04160A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
04170A2=A1-A2
04180J2=SQR(G*EXP(A2*(G-5/6)*A2/(G-2/3+.1/G)/(G-2/3+.1/G))-G)
04190GOSUB 6000
04200IF P3 <= .5 THEN 4230
04210P2=P3
04220GOTO 4250
04230J2=-J2
04240P=1-P
04250RETURN
04260REM ****************************************************
04270REM LOG GAMMA ROUTINE
04280REM INPUT G9
04290REM OUTPUT G0
04300G5=G9
04310IF G9 <= 1.E+30 THEN 4340
04320G0=1.E+38
04330RETURN
04340IF G9>1.E-09 THEN 4370
04350G0=0
04360RETURN
04370IF G9<1.E+10 THEN 4400
04380G0=G9*(LOG(G9)-1)
04390RETURN
04400G6=1
04410IF 18<G5 THEN 4450
04420G6=G6*G5
04430G5=G5+1
04440GOTO 4410
04450R8=1/G5/G5
04460G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
04470C1=8.33333E-02
04480C2=2.77778E-03
04490C3=7.93651E-04
04500C4=5.95238E-04
04510G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
04520RETURN
04530REM END OF LOG GAMMA ROUTINE
04540REM ****************************************************
06000REM****************************************************************
06001REM STUDENT'S T CDF ROUTINE
06002REM INPUT G J2
06003REM OUTPUT P
06004REM PRIOR GOSUB 6060
06005P=0
06006IF J2<.0001 THEN 6034
06007IF G<6 THEN 6011
06008IF J2<6 THEN 6012
06009P=1
06010GOTO 6035
06011IF J2>12 THEN 6009
06012 DIM W(16),O(16)
06013Y3=J2
06014IF G<10 THEN 6019
06015Y3=(G-2/3+.1/G)*SQR(ABS(LOG(1+J2*J2/G))/(G-5/6))
06016GOSUB 8000
06017GOTO 6035
06018REM PEIZER PRATT APPROXIMATION
06019IF G=1 THEN 6098
06020J1=0
06021N=G
06022P=0
06023GOSUB 6084
06024D0=(J2-J1)*.5
06025D1=(J1+J2)*.5
06026FOR I1=1 TO 16
06027D9=D0*O(I1)+D1
06028IF D9=0 THEN 6031
06029IF D9=1 THEN 6031
06030P=P+W(I1)*(EXP(-(N+1)/2*LOG(1+D9*D9/N)))
06031NEXT I1
06032P=P*F0
06033P=P*D0
06034P=P+.5
06035RETURN
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 4260
06088F0=G0
06090G9=N/2
06092GOSUB 4260
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*************************************************************
08000REM **********************************************************
08010REM ROUTINE CALCULATES THE CDF FOR NORMAL DISTRIBUTION
08020REM INPUT Y3
08030REM OUTPUT P
08040REM
08050Y4=ABS(Y3)
08060X1=X
08070X=Y3
08080T=1/(1+.231642*Y4)
08090D=.398942*EXP(-X*X/2)
08100C1=1.33027
08110C2=1.82126
08120C3=1.78148
08130C4=.356564
08140C5=.319382
08150P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
08160IF X >= 0 THEN 8180
08170P=1-P
08180X=X1
08190RETURN
08200REM
08210REM END OF NORMAL CDF ROUTINE
08220REM **********************************************************
09000REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.
09010INPUT O1
09020IF O1=-9999 THEN 9040
09030RETURN
09040CHAIN "RSTRT"
09050REM*************END ROUTINE
09060REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED. 2 INPUTS
09070INPUT O1,O2
09080IF O1=-9999 THEN 9110
09090IF O2=-9999 THEN 9110
09100RETURN
09110CHAIN "RSTRT"
09120REM*************END ROUTINE
09999END