Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod63.bas
There are 2 other files named cmod63.bas in the archive. Click here to see a list.
00020GOSUB 6060
00030X=0
00040REM****************************************************************
00050REM CMOD63 CMOD63 CMOD63 CMOD63 CMOD63
00060REM****************************************************************
00070 DIM A(10),B(4),Z(4,250),V(4,4),Q(1,1),F(4,4)
00080 DIM Y(250),N(12),X(250,4),C(4,1),D(10),E(4,4)
00090 DIM G(10,1),M(10,1),H(10,10),U(4),K(1,10),L(10,1),I(4,4),J(5)
00100D0=0
00110PRINT L$
00120PRINT " PREDICTED VALUES FOR SIMULTANEOUS ESTIMATION"
00130FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00170RESTORE#1
00171 INPUT# 1,I1,I2,I3
00180SCRATCH#1
00181 PRINT # 1,63,I2,I3
00190RESTORE#2
00191 INPUT# 2,M
00200IF M <> 0 THEN 260
00210PRINT "YOU DO NOT HAVE NECESSARY DATA IN RESTART FILE. YOU MUST"
00220PRINT "RUN MODULE 1 BEFORE YOU MAY FIND PREDICTED VALUES."
00230PRINT "TYPE '1' TO CONTINUE.";
00240GOSUB 9000
00250CHAIN "RSTRT"
00260 INPUT#2,P
00262INPUT#2,P0
00264INPUT#2,P1
00266INPUT#2,C0
00270MAT A=ZER(M)
00280MAT B=ZER(P)
00290 FOR I=1 TO M
00291INPUT#2,A(I)
00292NEXTI
00293 FOR I=1 TO P
00294 INPUT#2,B(I)
00295NEXTI
00296FORI=1TO5
00297INPUT#2,J(I)
00298NEXTI
00300RESTORE#4
00361 INPUT #4,N$
00370 INPUT#4,I2
00372 INPUT#4,I9
00374 INPUT#4,G$
00376 INPUT#4,V$
00380FORI=1TO12
00381 INPUT#4,N(I)
00382NEXTI
00390P5=P
00400S0=0
00410FOR I=1 TO M
00420IF N(I) <= 25 THEN 450
00430S0=S0+25
00440GOTO 460
00450S0=S0+N(I)
00460NEXT I
00470REM N NOW HAS GROUP COUNTS
00480MAT Y=ZER(S0)
00490I5=0
00500S3=0
00510MAT X=ZER(S0,P)
00520FOR I2=1 TO M
00530I6=0
00540FOR I1=1 TO I9
00550IF J(I1) <> 1 THEN 570
00560I6=I6+1
00570FOR I3=1 TO N(I2)
00580 INPUT#4,I4
00590IF I3>25 THEN 650
00600IF I1 <> C0 THEN 630
00610Y(S3+I3)=I4
00620GOTO 650
00630IF J(I1) <> 1 THEN 650
00640X(S3+I3,I6)=I4
00650NEXT I3
00660NEXT I1
00670IF N(I2) <= 25 THEN 690
00680N(I2)=25
00690S3=S3+N(I2)
00700NEXT I2
00710REM X NOW HAS DATA IN IT
00720S0=S3
00850MAT Z=ZER(P,S0)
00860S9=0
00870MAT M=ZER(M,1)
00880FOR I=1 TO M
00890FOR J=1 TO N(I)
00900M(I,1)=M(I,1)+Y(S9+J)
00910FOR K=1 TO P
00920Z(K,S9+J)=Z(K,S9+J)-X(S9+J,K)/S0
00930NEXT K
00940NEXT J
00950FOR J=1 TO N(I)
00960FOR K=1 TO P
00970Z(K,S9+J)=X(S9+J,K)+Z(K,S9+J)
00980NEXT K
00990NEXT J
01000S9=S9+N(I)
01010NEXT I
01020REM Z HAS ADJUSTED DATA IN IT
01030MAT U=ZER(P)
01040S3=0
01060FOR I=1 TO M
01070FOR J=1 TO N(I)
01080FOR K=1 TO P
01090U(K)=U(K)+Z(K,S3+J)
01100NEXT K
01110NEXT J
01120S3=S3+N(I)
01130NEXT I
01140FOR K=1 TO P
01150U(K)=U(K)/S0
01160NEXT K
01170FOR I=1 TO M
01180FOR K=1 TO P
01190A(I)=A(I)-U(K)*B(K)
01200NEXT K
01210NEXT I
01220REM A NOW HAS ALPHAS AT ORIGIN
01230MAT Z=TRN(X)
01240MAT E=ZER(P,P)
01250MAT V=ZER(P,P)
01260MAT E=Z*X
01270MAT V=INV(E)
01280MAT M=ZER(P,1)
01285GOTO 1300
01290PRINT L$
01300PRINT "TO OBTAIN PREDICTED VALUES FOR A GIVEN SET OF X VALUES,"
01310PRINT "TYPE '1'; OTHERWISE TYPE '0'"
01320GOSUB 9000
01330IF O1=0 THEN 1860
01335PRINT L$
01340PRINT "INPUT VALUE OF PREDICTORS"
01345K=0
01350FOR I=1 TO 5
01360IF J(I) <> 1 THEN 1400
01370 PRINT MID$(V$,(I-1)*6+1,(I-1)*6+6-((I-1)*6+1)+1);"=";
01375K=K+1
01380GOSUB 9000
01390M(K,1)=O1
01400NEXT I
01410PRINT " "
01420PRINT "GROUP PREDICTED Y"
01430PRINT " "
01440FOR I=1 TO M
01450Y=A(I)
01460K=0
01470FOR J=1 TO 5
01480IF J(J) <> 1 THEN 1510
01490K=K+1
01500Y=Y+B(K)*M(K,1)
01510NEXT J
01520L(I,1)=Y
01530IF ABS(Y)<.01 THEN 1590
01540IF ABS(Y)>9999.99 THEN 1590
01550:'CCCCC ####.##
01560 PRINT USING 1550,MID$(G$,(I-1)*6+1,(I-1)*6+6-((I-1)*6+1)+1),Y
01570GOTO 1610
01580:'CCCCC #.#####^^^^
01590 PRINT USING 1580,MID$(G$,(I-1)*6+1,(I-1)*6+6-((I-1)*6+1)+1),Y
01600L(I,1)=Y
01610NEXT I
01620MAT K=ZER(1,P)
01630MAT G=ZER(P,1)
01640MAT G=V*M
01650MAT Q=K*M
01660O9=P
01670G=S0-P
01680L0=(1+1/M+Q(1,1))*P0*(S0+3)
01685PRINT " "
01690PRINT "INPUT A VALUE Y0 FOR WHICH YOU WISH TO SEE THE PROBABILITY"
01692PRINT "FOR EACH GROUP THAT Y IS GREATER THAN Y0. FOR THE PREDICTORS"
01693K=0
01694FOR I=1 TO 5
01695IF J(I) <> 1 THEN 1698
01696K=K+1
01697 PRINT MID$(V$,(I-1)*6+1,(I-1)*6+6-((I-1)*6+1)+1);"=";M(K,1)
01698NEXT I
01700GOSUB 9000
01710Y=O1
01720PRINT "GROUP PROBABILITY Y>";Y
01730FOR I=1 TO M
01740J2=ABS(Y-L(I,1))/SQR(L0/(S0-O9))
01750GOSUB 6000
01760IF L(I,1)>Y THEN 1790
01770P=1-P
01780:'CCCCC ##.###
01790 PRINT USING 1780,MID$(G$,(I-1)*6+1,(I-1)*6+6-((I-1)*6+1)+1),P
01800NEXT I
01810PRINT "ENTER A NEW Y VALUE OR TYPE '-7777' TO CONTINUE.";
01820GOSUB 9000
01830P=O9
01840IF O1 <> -7777 THEN 1710
01850GOTO 1290
01860CHAIN "RSTRT"
05850 REM
05860G5=G9
05863IF G9 <= 1.E+30 THEN 5872
05866G0=1.E+38
05869RETURN
05872IF G9>1.E-09 THEN 5881
05875G0=C0
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
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