Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-04 - decus/20-0113/cmod65.bas
There are 2 other files named cmod65.bas in the archive. Click here to see a list.
00015REM *********************************************************************
00020REM  CMOD65   CMOD65   CMOD65    CMOD65   CMOD65    CMOD65   CMOD65
00025REM
00030REM      PRIOR FOR SIMPLE LINEAR REGRESSION
00035REM
00040REM******************************************************************
00045FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00065RESTORE#1
00066  INPUT#  1,I1,I2,I3
00070SCRATCH#1
00071  PRINT #  1,65,I2,I3
00075  DIM K(11),W(3),Y(3),A(3),B(3),C(3),V(3),L(3)        ,H(5)
00080  DIM D(2,5),E(5)
00085  S1$="------------------------------------------------------------"
00090H(1)=2.1318
00095H(3)=1.8595
00100H(5)=1.7247
00105E(1)=4
00110E(2)=6
00115E(3)=8
00120  DIM Z(7)
00125E(4)=12
00130E(5)=20
00135D(1,1)=.7407
00140D(1,2)=.7176
00150D(1,3)=.7064
00155D(1,4)=.6995
00160D(1,5)=.687
00165D(2,1)=1.5332
00170  E1$="REENTER.  75TH PERCENTILE MUST BE GREATER THAN ESTIMATE."
00175D(2,2)=1.4398
00180D(2,3)=1.3968
00185D(2,4)=1.3562
00190D(2,5)=1.3253
00195DATA 2.45262,-2.47084,2.01371,-2.0064
00200READ A1,B1,A2,B2
00205S7=0
00210PRINT L$
00215PRINT "    SPECIFICATION OF A PRIOR FOR SIMPLE LINEAR REGRESSION"
00220PRINT
00225PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A PRIOR DISTRIBUTION "
00230PRINT "TO YOUR BELIEFS ABOUT THE PARAMETERS OF THE SIMPLE  LINEAR "
00235PRINT "REGRESSION MODEL. "
00240PRINT
00245PRINT "VARIABLES:      X = PREDICTOR OR INDEPENDENT VARIABLE."
00250PRINT "                Y = CRITERION OR DEPENDENT VARIABLE."
00255PRINT
00260PRINT "MODEL: Y GIVEN X IS NORMALLY DISTRIBUTED WITH"
00265PRINT
00270PRINT "             MEAN = ALPHA + BETA * X"
00275PRINT "         ST. DEV. = SIGMA"
00280PRINT
00285PRINT "PARAMETERS: ALPHA = Y INTERCEPT AT X = 0"
00290PRINT "            BETA  = SLOPE OF THE REGRESSION LINE"
00295PRINT "            SIGMA = RESIDUAL STANDARD DEVIATION"
00300PRINT
00305PRINT "THE X VALUES MAY BE FIXED OR RANDOMLY SAMPLED.  IF SAMPLED,"
00310PRINT "THIS MAY BE SEPARATE FROM, OR JOINTLY WITH Y."
00320PRINT
00325PRINT "WHEN YOU ARE READY TO CONTINUE, TYPE '1'";
00330GOSUB 9000
00335PRINT L$
00340PRINT "THE FIRST STEP IN THE SPECIFICATION OF A PRIOR CONSISTS "
00345PRINT "OF FITTING A LINE TO YOUR BEST ESTIMATES OF Y GIVEN X"
00350PRINT "FOR SEVERAL DIFFERENT VALUES OF X."
00355PRINT
00360PRINT "LET X2 BE THE VALUE OF X FOR WHICH YOU FEEL MOST CONFIDENT"
00365PRINT "EXPRESSING YOUR BELIEFS ABOUT Y GIVEN X."
00370PRINT "X2 =";
00375GOSUB 9000
00380X(2)=O1
00385PRINT
00390PRINT "LET X1 BE THE SMALLEST VALUE OF X FOR WHICH YOU FEEL YOU CAN"
00395PRINT "COMFORTABLY EXPRESS YOUR BELIEFS ABOUT Y GIVEN X."
00400PRINT "X1 =";
00405GOSUB 9000
00410IF O1<X(2) THEN 425
00415PRINT "REENTER.  X1 MUST BE SMALLER THAN X2."
00420GOTO 405
00425X(1)=O1
00430PRINT
00435PRINT "LET X3 BE THE LARGEST VALUE OF X FOR WHICH YOU FEEL COMFORTABLE"
00440PRINT "EXPRESSING YOUR BELIEFS ABOUT Y GIVEN X."
00445PRINT "X3 =";
00450GOSUB 9000
00455IF X(2)<O1 THEN 470
00460PRINT "REENTER.  X3 MUST BE GREATER THAN X2."
00465GOTO 450
00470X(3)=O1
00475PRINT L$
00480PRINT
00485PRINT "THESE THREE VALUES OF X WILL SERVE AS REFERENCE POINTS IN"
00490PRINT "OBTAINING YOUR BEST ESTIMATE OF THE TRUE REGRESSION LINE."
00495PRINT
00500FOR I=1 TO 3
00505:          X## = #######.##
00510PRINT  USING 505,I,X(I)
00515NEXT I
00520PRINT
00525PRINT "     X                 BEST ESTIMATE OF Y"
00530FOR I=1 TO 3
00535GOSUB 550
00540NEXT I
00545GOTO 635
00550PRINT"    ";X(I);"                    ";
00560GOSUB 9000
00565IF I <> 2 THEN 585
00570IF Y(1) <> O1 THEN 585
00575PRINT "REENTER.  BEST ESTIMATES MUST BE DIFFERENT FOR DIFFERENT X'S."
00580GOTO 550
00585Y(I)=O1
00590IF I <> 3 THEN 630
00595IF Y(1)>Y(2) THEN 615
00600IF Y(3)>Y(2) THEN 630
00605PRINT "REENTER.  REGRESSION LINE IS INCREASING AS INCREASES."
00610GOTO 550
00615IF Y(2)>Y(3) THEN 630
00620PRINT "REENTER.  REGRESSION LINE IS DECREASING AS X INCREASES."
00625GOTO 550
00630RETURN
00635PRINT L$
00640PRINT "YOUR PRIOR DISTRIBUTION ON Y GIVEN X IS ASSUMED TO BE"
00645PRINT "SYMMETRIC ABOUT YOUR BEST ESTIMATE OF Y.  YOUR BEST "
00650PRINT "ESTIMATE IS THE 50TH PERCENTILE OF YOUR PRIOR."
00655PRINT
00660PRINT "A FURTHER ASSUMPTION IS THAT THE CLOSER X IS TO X2 THE"
00665PRINT "SMALLER IS THE ABSOLUTE DIFFERENCE BETWEEN THE 50TH AND"
00670PRINT "75TH PERCENTILES OF YOUR PRIOR DISTRIBUTION ON Y GIVEN X."
00675PRINT
00680PRINT "INPUT THE 75TH PERCENTILE OF YOUR PRIOR DISTRIBUTION ON"
00685PRINT "Y GIVEN EACH OF THE FOLLOWING X VALUES."
00690PRINT
00695PRINT"                         ABS(X-X2)  50TH         75TH"
00700K4=2
00705K5=2
00710GOSUB 735
00715K4=1
00720K5=3
00725GOSUB 735
00730GOTO 935
00735FOR I=K4 TO K5 STEP 2
00736PRINT"    X";
00737SCRATCH#9
00738PRINT#9,I
00739RESTORE#9
00740INPUT#9,Z$
00741PRINTZ$;
00742PRINT" = ";
00743SCRATCH#9
00744PRINT#9,X(I)
00745RESTORE#9
00746INPUT#9,Z$
00747Z$=RIGHT$("         "+Z$,10)
00748PRINTZ$;
00749SCRATCH#9
00750PRINT#9,ABS(X(2)-X(I))
00751RESTORE#9
00752INPUT#9,Z$
00753Z$=RIGHT$("          "+Z$,10)
00754PRINTZ$;
00755SCRATCH#9
00756PRINT#9,Y(I)
00757RESTORE#9
00758INPUT#9,Z$
00759Z$=RIGHT$("         "+Z$,10)
00760PRINTZ$;
00761PRINT"          ";
00762GOSUB9000
00763IF O1>Y(I)THEN770
00764PRINT E1$
00765 GOTO 762
00770W(I)=O1-Y(I)
00775IF I=2 THEN 845
00780IF W(I)>W(2) THEN 800
00785PRINT "REENTER.  DIFFERENCE BETWEEN 75TH AND 50TH MUST BE GREATER"
00790PRINT "FOR THIS X THAN IT IS FOR X2."
00795 GOTO 762
00800IF I=1 THEN 845
00805IF ABS(X(1)-X(2))=ABS(X(2)-X(3)) THEN 845
00810IF ABS(X(1)-X(2))<ABS(X(2)-X(3)) THEN 840
00815IF W(1) >= W(3) THEN 845
00820PRINT "REENTER. THE ASSUMPTION IS THAT THE FARTHER X IS FROM"
00825PRINT "X2 THE GREATER THE DIFFERENCE BETWEEN THE 75TH AND"
00830PRINT "50TH PERCENTILES."
00835GOTO 750
00840IF W(3) <= W(1) THEN 820
00845NEXT I
00850RETURN
00855GOTO 935
00935FOR I=1 TO 3
00940V(I)=1/W(I)/W(I)
00945NEXT I
00950N=3
00955X0=0
00960Y0=0
00965V0=0
00970FOR I=1 TO N
00975V=V(I)
00980X0=X0+V*X(I)
00985Y0=Y0+V*Y(I)
00990V0=V0+V
00995NEXT I
01000X0=X0/V0
01005Y0=Y0/V0
01010S1=0
01015S3=0
01020FOR I=1 TO N
01025V=V(I)
01030X=X(I)-X0
01035Y=Y(I)-Y0
01040S1=S1+V*X*X
01045S3=S3+V*X*Y
01050NEXT I
01055B=S3/S1
01060A=Y0-B*X0
01065PRINT L$
01070PRINT "HERE ARE PERCENTILES OF THE DISTRIBUTIONS FITTED TO YOUR"
01075PRINT "SPECIFIED PERCENTILES."
01080PRINT
01085PRINT "                                Y GIVEN X"
01090PRINT "             ------------------------------------------------"
01095PRINT"                 FITTED    SPECIFIED FITTED SPECIFIED FITTED"
01100PRINT"                  25TH       50TH    50TH    75TH     75TH"
01105PRINT S1$
01110FOR I=1 TO 11
01115K(I)=I
01120NEXT I
01125V0=0
01130FOR I=1 TO N-1
01135IF X(K(I))>X(K(I+1)) THEN 1145
01140GOTO 1165
01145K0=K(I)
01150K(I)=K(I+1)
01155K(I+1)=K0
01160V0=1
01165NEXT I
01170IF V0 <> 0 THEN 1125
01175GOSUB 1525
01180FOR K0=1 TO N
01185I=K(K0)
01190W=SQR(C+(X(I)-X2)*(X(I)-X2)*D)
01195Y1=A+B*X(I)
01200:X##=#####.## #####.## #####.## ######.## #####.## #####.##
01205PRINT  USING 1200,I,X(I),Y1-W,Y(I),Y1,Y(I)+W(I),Y1+W
01210NEXT K0
01215E0=0
01220E1=0
01225PRINT
01230PRINT "IF YOU WANT TO CHANGE SOME OF THE SPECIFIED PERCENTILES"
01235PRINT "TYPE '1', ELSE '0'.";
01240GOSUB 9000
01245IF O1 <> 1 THEN 1370
01250PRINT
01255MAT A=V
01260MAT B=W
01265MAT C=Y
01270GOTO 1380
01275IF W(1)>W(2) THEN 1285
01280GOTO 1290
01285IF W(3)>W(2) THEN 1345
01290PRINT "THE CHANGES YOU HAVE MADE HAVE RESULTED IN A SET OF PERCENTILES"
01295PRINT "THAT VIOLATE THE ASSUMPTION THAT THE FARTHER X IS FROM X2 THE"
01300PRINT "GREATER IS THE ABSOLUTE DIFFERENCE BETWEEN THE 50TH AND 75TH."
01305PRINT "THE ABOVE PERCENTILES WILL BE REDISPLAYED AND YOU WILL HAVE"
01310PRINT "THE CHANCE TO MAKE OTHER CHANGES."
01315PRINT
01320MAT V=A
01325MAT Y=C
01330MAT W=B
01335PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
01340GOSUB 9000
01345IF Y(1)>Y(2) THEN 1360
01350IF Y(3) <= Y(2) THEN 1290
01355GOTO 955
01360IF Y(3)>Y(2) THEN 1290
01365GOTO 955
01370IF O1=0 THEN 1510
01375GOTO 1500
01380E0=1
01385PRINT
01390PRINT "INPUT X, 50TH, AND 75TH PERCENTILES (NONE=0,0,0).";
01395GOSUB 9100
01400X=O1
01405IF O2 <> O3 THEN 1415
01410IF O2=0 THEN 1275
01415FOR I=1 TO N
01420IF X=X(I) THEN 1440
01425NEXT I
01430PRINT "REENTER. X MUST BE ONE OF THOSE ORIGINALLY SUPPLIED."
01435GOTO 1395
01440IF O3>O2 THEN 1455
01445PRINT E1$
01450GOTO 1395
01455W(I)=O3-O2
01460Y(I)=O2
01465V(I)=1/W(I)/W(I)
01470GOTO 1390
01475GOTO 1485
01480GOSUB 635
01485PRINT
01490PRINT "IF YOU WOULD LIKE TO ENTER ANOTHER VALUE OF X, TYPE '0'."
01495GOTO 1240
01500PRINT "REENTER.  INPUT MUST BE 0 OR 1"
01505GOTO 1240
01510IF E0=1 THEN 955
01515GOTO 1650
01520IF E1=1 THEN 1650
01525X0=0
01530Y0=0
01535V0=0
01540X2=X(2)
01545FOR I=1 TO N
01550V=V(I)*V(I)
01555X0=X0+(X(I)-X2)*(X(I)-X2)*V
01560Y0=Y0+V(I)
01565V0=V0+V
01570NEXT I
01575X0=X0/V0
01580Y0=Y0/V0
01585S1=0
01590S3=0
01595FOR I=1 TO N
01600V=V(I)*V(I)
01605X=(X(I)-X2)*(X(I)-X2)-X0
01610Y=W(I)*W(I)-Y0
01615S1=S1+V*X*X
01620S3=S3+V*X*Y
01625NEXT I
01630D=S3/S1
01635C=Y0-D*X0
01640RETURN
01645PRINT L$
01650PRINT L$
01655PRINT "THIS MODEL ASSUMES THAT YOUR UNCERTAINTY ABOUT A RANDOMLY"
01660PRINT "SAMPLED Y GIVEN X HAS TWO COMPONENTS.  FIRST, YOU ARE NOT"
01665PRINT "SURE WHAT THE MEAN OF Y GIVEN X IS, AND SECOND, THERE IS"
01670PRINT "THE UNCERTAINTY, THE RESIDUAL UNCERTAINTY, THAT WOULD BE"
01675PRINT "THERE EVEN IF YOU KNEW THE MEAN OF Y GIVEN X. "
01680PRINT
01685PRINT "WE NOW WANT TO FIND OUT HOW MUCH OF YOUR UNCERTAINTY ABOUT"
01690PRINT "Y GIVEN X CAN BE ATTRIBUTED TO YOUR UNCERTAINTY ABOUT THE"
01695PRINT "MEAN OF Y GIVEN X."
01700PRINT
01705PRINT "THE 50TH PERCENTILES OF YOUR PRIORS ON Y GIVEN X AND THE"
01710PRINT "MEAN OF Y GIVEN X ARE ASSUMED TO BE EQUAL.  THE 75TH OF"
01715PRINT "YOUR PRIOR ON Y GIVEN X IS ASSUMED TO BE GREATER THAN THE"
01720PRINT "75TH  OF YOUR PRIOR ON THE MEAN OF Y GIVEN X, REFLECTING"
01725PRINT "YOUR GREATER UNCERTAINTY ABOUT Y GIVEN X."
01730PRINT
01735:YOUR PRIOR FOR A RANDOMLY SAMPLED Y GIVEN X2=#######.##
01740PRINT  USING 1735,X(2)
01745:HAS  50TH =#######.## AND 75TH =#######.##.
01750PRINT  USING 1745,A+B*X(2),A+B*X(2)+SQR(C)
01755PRINT
01760PRINT "INPUT THE 75TH OF YOUR PRIOR ON THE MEAN OF Y GIVEN X2.";
01765GOSUB 9000
01770IF O1>A+B*X(2) THEN 1785
01775PRINT "REENTER.  75TH MUST BE GREATER THAN 50TH."
01780GOTO 1765
01785IF O1<A+B*X(2)+SQR(C) THEN 1800
01790PRINT "REENTER.  75TH FOR Y MUST BE GREATER THAN 75TH FOR MEAN."
01795GOTO 1765
01800S=O1
01805V6=(O1-A-B*X2) ** 2/.676/.676
01810V6=C/.676/.676-V6
01815PRINT L$
01820IF S7=1 THEN 1870
01825S7=1
01830PRINT "HERE ARE SOME PERCENTILES OF YOUR PRIOR DISTRIBUTION ON"
01835PRINT "THE MEAN OF Y GIVEN X."
01840PRINT
01845:FOR X =#######.##               75TH =########.##
01850PRINT  USING 1845,X(2),O1
01855:                         75TH - 50TH =########.##
01860PRINT  USING 1855,O1-A-B*X(2)
01865PRINT
01870PRINT "                             MEAN OF Y GIVEN X"
01875PRINT "                 ---------------------------------------"
01880PRINT "     X             25TH            50TH           75TH"
01885PRINT "---------------------------------------------------------"
01890FOR K0=1 TO N
01895I=K(K0)
01900Y=A+B*X(I)
01905W=.676*SQR((C+(X(I)-X2)*(X(I)-X2)*D)/.676/.676-V6)
01910:#####.##       #####.##       #####.##       #####.##
01915PRINT  USING 1910,X(I),Y-W,Y,Y+W
01920NEXT K0
01925PRINT
01930PRINT "YOU CAN CHANGE THE 75TH PERCENTILE FOR X2 IF YOU WANT."
01935PRINT "A SMALLER VALUED 75TH PERCENTILE FOR X2 WILL RESULT IN"
01940PRINT "SMALLER VALUED 75TH PERCENTILES FOR THE OTHER X'S."
01945PRINT
01950PRINT "IF YOU WANT TO CHANGE THE 75TH FOR X2 TYPE '1', ELSE '0'.";
01955GOSUB 9000
01960IF O1=0 THEN 1980
01965IF O1=1 THEN 1755
01970PRINT "REENTER. INPUT MUST BE 0 OR 1."
01975GOTO 1955
01980REM
01985Z(1)=X2
01990Z(2)=A+B*X2
01995Z(4)=S*S/D
02000Z(5)=B*S*S/D
02005PRINT L$
02010PRINT "HERE ARE THE 50TH, 75TH, 90TH, AND 95TH PERCENTILES OF "
02015PRINT "DISTRIBUTIONS FITTED TO YOUR PERCENTILE SPECIFICATIONS "
02020PRINT "FOR THE MEAN OF Y GIVEN X2."
02025PRINT
02030Y3=A+B*X2
02035:X2 = #######.##
02040PRINT  USING 2035,X(2)
02045PRINT
02050PRINT "                        MEAN OF Y GIVEN X2"
02055PRINT "          ---------------------------------------------------"
02060PRINT "DIST.     50TH          75TH          90TH         95TH"
02065K=0
02070FOR I=1 TO 5 STEP 2
02075G=E(I)
02080V7=C*E(I)/D(1,I)/D(1,I)
02085V8=G*(S-A-B*X(2)) ** 2/D(1,I)/D(1,I)
02090P3=D(2,I)/D(1,I)*(S-A-B*X(2))
02095P4=P3*H(I)/D(1,I)
02100:##    #######.##    #######.##    #######.##    #######.##
02105K=K+1
02110PRINT  USING 2100,K,Y3,S,Y3+P3,Y3+P4
02115NEXT I
02120PRINT
02125PRINT "TYPE THE NUMBER OF THE DISTRIBUTION WHICH BEST REPRESENTS YOUR"
02130PRINT "BELIEFS ABOUT THE MEAN OF Y GIVEN X2.";
02135GOSUB 9000
02136 O1=INT(O1)
02137 IF O1<1 THEN 2145
02138 IF O1>5 THEN 2145
02140  ONO1  GOTO 2155,2155,2155,2155,2155
02145PRINT "REENTER.  INPUT MUST BE NUMBER OF DISTRIBUTION."
02150GOTO 2135
02155I=O1*2-1
02160Z(6)=E(I)
02165G=E(I)
02170S=S-A-B*X(2)
02175Z(7)=C*G/D(1,I)/D(1,I)-S*S*G/D(1,I)/D(1,I)
02180 Z(3)=Z(7)/(S*S*G/D(1,I)/D(1,I) ) 
02183 A8=(C+D)*G/D(1,I)/D(1,I)
02184 A8=A8/Z(7)-1-1/Z(3)
02185PRINT L$
02186 Z(4)=1/A8
02187 Z(5)=Z(4)*B
02190PRINT "THIS COMPLETES THE SPECIFICATION OF PRIOR DISTRIBUTIONS."
02195:  INTERCEPT  X=######.## ##### #######.##   ########.##
02200 PRINT
02205PRINT "                      DEGREES OF"
02210PRINT "                        FREEDOM    MEAN      SCALE PARAMETER"
02215PRINT S1$
02220PRINT "INVERSE CHI DISTRIBUTION"
02225:  RESIDUAL ST.DEV.     ####### #######.##   ########.##
02230PRINT  USING 2225,Z(6),SQR(Z(7))/SQR(Z(6)-1.5),SQR(Z(7))
02235PRINT
02240PRINT "STUDENT'S T DISTRIBUTIONS"
02245:  SLOPE                ####### ######.##    ########.## 
02250 PRINT USING 2245 ,Z(6),Z(5)/Z(4),Z(7)/Z(4)
02255FOR I=1 TO 3
02260M1=A+B*X(I)
02265S0=Z(7)/Z(3)+(X(I)-X(2)) ** 2/Z(4)*Z(7)
02270S1=SQR(S0/(Z(6)-2))
02275PRINT  USING 2195,X(I),Z(6),M1,S0
02280NEXT I
02285PRINT
02290PRINT "IF YOU WANT TO DO THE POSTERIOR ANALYSIS TYPE '1', ELSE '0'.";
02295GOSUB 9000
02300IF O1=1 THEN 2325
02305IF O1=0 THEN 2320
02310PRINT "REENTER.  INPUT MUST BE 1 OR 0."
02315GOTO 2295
02320GOTO 4500
02325SCRATCH#2
02326FORI=1TO7
02327PRINT#2,Z(I)
02328NEXTI
02330CHAIN "CMOD66"
03000REM***********************************************************
03010REM           CHI-SQUARE    PERCENTILE
03020REM                INPUTS P2
03030REM                OUTPUTS  J2
03040J=P0
03050P9=P0
03060P2=P0
03070GOSUB 3330
03080X=J2
03090GOSUB 5460
03100IF ABS(P-P9)<.0001 THEN 3320
03110P4=P
03120E7=J2
03130P2=P2-.009
03140GOSUB 3330
03150E1=J2
03160E2=E7
03170X=E1
03180GOSUB 5460
03190P3=P
03200X3=(P9-P3)/(P4-P3)
03210X3=X3*(E2-E1)+E1
03220X=X3
03230GOSUB 5460
03240IF ABS(P9-P)<.0001 THEN 3320
03250IF P>P9 THEN 3290
03260E1=X3
03270P3=P
03280GOTO 3200
03290E2=X3
03300P4=P
03310GOTO 3200
03320RETURN
03330P3=P2
03340IF P2 <= .5 THEN 3360
03350P2=1-P2
03360A1=SQR(LOG(1/P2/P2))
03370A2=2.51552+.802853*A1+.010328*A1*A1
03380A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
03390J2=A1-A2
03400IF P3 <= .5 THEN 3430
03410P2=P3
03420GOTO 3440
03430J2=-J2
03440J2=G*(1-2/9/G+J2*SQR(2/9/G)) ** 3
03450RETURN
03460REM
03470REM*********************  END OF CSQ PERCENTILE ********************
03480REM   ********************************************************
03490REM            SUBROUTINE FOR COMPUTING
03500REM            DENSITY OF INVERSE CHI
03510REM   ********************************************************
03520G9=.5*G
03530GOSUB 5750
03540F0=LOG(G0)
03550F0=LOG(2)-F0+.5*G*LOG(.5*S0*S0)-(G+1)*LOG(J2)-.5*S0*S0/J2 ** 2
03560D2=F0
03570RETURN
04500PRINT L$
04510PRINT "HERE IS THE INFORMATION YOU WILL NEED TO INPUT IF"
04521PRINT "YOU LATER DECIDE TO DO THE POSTERIOR ANALYSIS."
04522PRINT
04540:DEGREES OF FREEDOM          =#######.##
04550PRINT  USING 4540,Z(6)
04560:MEAN OF PREDICTOR           =#######.##
04570PRINT  USING 4560,Z(1)
04580:MEAN OF CRITERION           =#######.##
04590PRINT  USING 4580,Z(2)
04600:SLOPE                       =#######.##
04610PRINT  USING 4600,Z(5)/Z(4)
04620:TOTAL OBSERVATIONS          =#######.##
04630PRINT  USING 4620,Z(3)
04640:SCALE FACTOR                =#######.##
04650PRINT  USING 4640,Z(7)
04660:WEIGHTING FACTOR            =#######.##
04670PRINT  USING 4660,Z(4)
04680PRINT
04690PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
04700GOSUB 9000
04710CHAIN "RSTRT"
05460REM ********************************************************
05470REM          CHI-SQUARE CDF ROUTINE-LOWER TAIL
05480REM                INPUT     G        X
05490REM                OUTPUT    P
05500REM*************************
05510REM   PEIZER-PRATT-APPROXIMATION
05520REM   PROB(CHISQ<X) WITH DF=G
05530REM*************************
05540X2=X/2
05550T=X2-G*.5+1/3-.04/G
05560Y=(G*.5-.5)/X2
05570IF Y<.00001 THEN 5640
05580IF ABS(Y-1)<.00001 THEN 5660
05590G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05600IF G6>-1 THEN 5670
05610Y3=T
05620GOTO 5680
05630GOTO 5670
05640G6=1
05650GOTO 5670
05660G6=0
05670Y3=T*SQR((1+G6)/X2)
05680GOSUB 8000
05690P2=P
05700RETURN
05710REM ****************************************************
05720REM        LOG GAMMA ROUTINE
05730REM           INPUT G9
05740REM           OUTPUT G0
05750G5=G9
05760IF G9 <= 1.E+30 THEN 5790
05770G0=1.E+38
05780RETURN
05790IF G9>1.E-09 THEN 5820
05800G0=0
05810RETURN
05820IF G9<1.E+10 THEN 5850
05830G0=G9*(LOG(G9)-1)
05840RETURN
05850G6=1
05860IF 18<G5 THEN 5900
05870G6=G6*G5
05880G5=G5+1
05890GOTO 5860
05900R8=1/G5/G5
05910G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
05920C1=8.33333E-02
05930C2=2.77778E-03
05940C3=7.93651E-04
05950C4=5.95238E-04
05960G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
05970RETURN
05980REM          END OF LOG GAMMA ROUTINE
05990REM ****************************************************
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
09050REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.  2 INPUTS
09055INPUT O1,O2
09065IF O1=-9999 THEN 9080
09070IF O2=-9999 THEN 9080
09075RETURN
09080CHAIN "RSTRT"
09090REM*************END ROUTINE
09100REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.  3 INPUTS
09105INPUT O1,O2,O3
09115IF O1=-9999 THEN 9135
09120IF O2=-9999 THEN 9135
09125IF O3=-9999 THEN 9135
09130RETURN
09135CHAIN "RSTRT"
09145REM*************END ROUTINE
09999END