Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod84.bas
There are 2 other files named cmod84.bas in the archive. Click here to see a list.
00030REM*****************************************************************
00040 DIM C(5),F(5),V(5),G(5,5),H(5,5),U(1)
00050REM CMOD84 CMOD84 CMOD84 CMOD84 CMOD84
00060REM*****************************************************************
00070 DIM J(5),Q(6),X(333,5) ,N(12)
00080 DIM T(8),S(5,5),A(5,5)
00090GOSUB 5160
00100 DIM K(5)
00120PRINT L$
00130PRINT " HIGHEST DENSITY REGIONS"
00140PRINT
00141PRINT "FOR ANY SET OF VALUES FOR THE INTERCEPT AND THE REGRESSION"
00142PRINT "COEFFICIENTS THE MODULE WILL DETERMINE THE SMALLEST HDR"
00143PRINT "THAT CONTAINS THE SET OF VALUES. ANY LARGER HDR WILL ALSO"
00144PRINT "CONTAIN THE SET OF VALUES."
00145PRINT
00150PRINT "HERE ARE THE INTERCEPT AND REGRESSION COEFFICIENTS."
00160PRINT
00170 FILES RFILE1,RFILE2,RFILE3,RF4,,,RF7,RF8
00180X=0
00220RESTORE#1
00221 INPUT# 1,I1,I2,I3
00240RESTORE#2
00241 INPUT# 2,C7,G7
00245GOSUB 520
00250K=0
00260MAT T=ZER(P+3)
00270MAT Q=ZER(P+1)
00280 RESTORE#8
00281 FOR I=1 TO P+3
00282 INPUT#8,T(I)
00283 NEXT I
00290MAT J=ZER(P)
00295 K$=""
00300FOR I=1 TO P
00310IF T(I)=0 THEN 360
00320K=K+1
00330 K$=K$+MID$(V$,I*6-5,6)
00340P(K)=I
00350J(I)=1
00360NEXT I
00370I5=K
00380MAT S=ZER(P,P)
00390MAT A=ZER(I5+1,I5+1)
00400MAT F=ZER(I5+1)
00410K=1
00420FOR I=1 TO P
00430IF J(I)=0 THEN 460
00440K=K+1
00450F(K)=T(I)
00460NEXT I
00470F(1)=T(P+3)
00480MAT A=ZER(I5+1,I5+1)
00490 RESTORE #7
00491 FOR I=1 TO P
00492 FOR J=1 TO P
00493 INPUT#7,S(I,J)
00494NEXT J
00495 NEXT I
00500GOSUB 730
00510GOSUB 970
00520RESTORE#4
00521INPUT#4,N$
00522 INPUT#4,M
00523INPUT#4,P
00524 INPUT#4,G$
00525 INPUT#4,V$
00530 FORI=1 TO 12
00531 INPUT#4,N(I)
00532 NEXT I
00540IF M=0 THEN 650
00550IF G7=1 THEN 650
00560S0=0
00570FOR I=1 TO G7-1
00580S0=S0+N(I)
00590NEXT I
00600FOR I=1 TO S0*P
00610 INPUT#4,C3
00620NEXT I
00630S0=N(G7)
00640GOTO 660
00650S0=N(1)
00660MAT X=ZER(S0,P)
00670FOR I=1 TO P
00680FOR J=1 TO S0
00690 INPUT#4,X(J,I)
00700NEXT J
00710NEXT I
00720RETURN
00730K0=0
00740K4=0
00750A(1,1)=S(1,1)
00760FOR I=1 TO P
00770K1=0
00780K3=0
00790IF I <> C7 THEN 810
00800K4=1
00810IF J(I)=1 THEN 840
00820K0=K0+1
00830GOTO 950
00840FOR J=1 TO P
00850IF C7 <> J THEN 870
00860K3=K3+1
00870IF J(J)=1 THEN 900
00880K1=K1+1
00890GOTO 940
00900REM
00910A(I+1-K0,J+1-K1)=S(I+1-K4,J+1-K3)
00920A(1,J+1-K1)=S(1,J+1-K3)
00930A(I+1-K0,1)=S(I+1-K4,1)
00940NEXT J
00950NEXT I
00960RETURN
00970S9=0
00980:INTERCEPT =#######.###
00990PRINT USING 980,T(P+3)
01000FOR I=1 TO P
01010IF J(I)=0 THEN 1040
01020:'CCCCC =#######.###
01030 PRINT USING 1020,MID$(V$,I*6-5,I*6-(I*6-5)+1),T(I)
01040NEXT I
01050PRINT " "
01060MAT V=ZER(I5+1)
01070M=I5+1
01072PRINT "INPUT SET OF VALUES IN THE SAME ORDER AS THEY APPEAR ABOVE."
01080PRINT "(INTERCEPT,COEFFICIENTS)";
01090MAT INPUT V
01100MAT C=ZER(M)
01110FOR I=1 TO M
01120C(I)=V(I)-F(I)
01130NEXT I
01140MAT G=ZER(1,M)
01150MAT H=ZER(1,M)
01160MAT G=TRN(C)
01170MAT H=G*A
01190MAT U=H*C
01200G0=S0-M
01210S2=T(P+1)*T(P+1)*(G0-.6667)
01230F=U(0)/((S2/G0)*M)
01300A=M/2
01310B=G0/2
01320J1=0
01330J2=F*A/(B+F*A)
01340P4=P
01350M9=G0
01360GOSUB 5220
01370GOSUB 1500
01380G0=M9
01390IF P>0 THEN 1410
01400P=0
01410IF P <= .99 THEN 1440
01420PRINT "SET IS NOT IN THE 99% HDR."
01422PRINT
01430GOTO 1460
01440:SMALLEST HDR CONTAINING SET =###.#
01450PRINT USING 1440,P*100
01452PRINT
01460PRINT "IF YOU WISH TO ENTER ANOTHER SET TYPE '1', ELSE '0'.";
01470GOSUB 9000
01480P=P4
01490IF O1=1 THEN 1080
01492CHAIN "CMOD83"
01500REM ***************************************************
05005REM BETA CDF ROUTINE
05010REM INPUT A B J1 J2
05015REM OUTPUT P
05020REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220
05025 DIM W(16),O(16)
05030IF J1 <> 0 THEN 5035
05031IF A<5 THEN 5035
05032IF B >= 5 THEN 5400
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A>1 THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=0 THEN 5115
05105IF D9=1 THEN 5115
05107D9=LOG(D9)*(A-1)+LOG(1-D9)*(B-1)
05108IF D9<-80 THEN 5115
05110P=P+W(I1)*EXP(D9)
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=0 THEN 5155
05135A=C6
05140B=C7
05145P=1-P
05150J2=1-J2
05155RETURN
05160FOR I1=1 TO 16
05165READ W(I1),O(I1)
05170NEXT I1
05175DATA 2.71525E-02,-.989401
05180DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
05185DATA .124629,-.755404,.149596,-.617876
05190DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
05195DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
05200DATA .149596,.617876,.124629,.755404
05205DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
05210DATA .989401
05215RETURN
05220G9=A+B
05225GOSUB 5850
05230F0=G0
05235G9=A
05240GOSUB 5850
05245F0=F0-G0
05250G9=B
05255GOSUB 5850
05260F0=F0-G0
05265IF A+B>85 THEN 5275
05270F0=EXP(F0)
05275RETURN
05280W1=(B*J2)**(1/3)
05285W2=(A*(1-J2))**(1/3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)**(1/3)
05305W2=(A*(1-J1))**(1/3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05340REM 2/16/76 CHANGED TO ALL LOG
05345D2=LOG(F0)+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350RETURN
05390REM
05400REM *******************************************************
05403REM PEIZER PRATT APROXIMATION
05406D0=.02*(J2/B-(1-J2)/A+(J2-.5)/(A+B))
05409D0=B-1/3-(A+B-2/3)*(1-J2)+D0
05412Y=(B-.5)/((A+B-1)*(1-J2))
05415GOSUB 5436
05418D9=G6
05421Y=(A-.5)/((A+B-1)*J2)
05424GOSUB 5436
05427Y3=D0*SQR((1+J2*D9+(1-J2)*G6)/((A+B-5/6)*(1-J2)*J2))
05430GOSUB 8000
05433RETURN
05436IF Y<.00001 THEN 5451
05439IF ABS(Y-1)<.00001 THEN 5457
05442G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05445IF G6>-1 THEN 5460
05446G6=10**37
05448GOTO 5460
05451G6=1
05454GOTO 5460
05457G6=0
05460RETURN
05465REM END BETA CDF ROUTINE
05470REM***********************************************************
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 ****************************************************
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)
08021 IF X*X/2<80 THEN 8025
08022 D=0
08023 GOTO 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
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
09999 END