Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmod87.bas
There are 2 other files named cmod87.bas in the archive. Click here to see a list.
00020Q9=0
00030REM**************************************************************
00040REM CMOD87 CMOD87 CMOD87 CMOD87 CMOD87
00050REM INTERMEDIATE STEP FOR ANOVA LINEAR COMBINATIONS.
00060REM**************************************************************
00070 DIM M(12),N(12),S(12),A(12),C(12) ,P(6,6),Q(6),U(1,1),G(1,12)
00080 DIM L(12,12),I(12,12),V(12),T(12,12),B(12),E(12,12),F(12)
00090 C$="11122122"
00100GOSUB 5160
00110 FILES RFILE1,RFILE2,RFILE3,,,,RF7,RF8
00150RESTORE#1
00151 INPUT# 1,I1,I2,I3
00160SCRATCH#1
00161 PRINT # 1,87,I2,I3
00170I2=I3
00180S7=0
00190X=0
00200S9=0
00210IF I1=75 THEN 260
00220RESTORE#8
00221 INPUT# 8 ,G
00230IF G=0 THEN 390
00240I2=0
00250Q9=1
00260RESTORE#7
00261INPUT#7,M
00262INPUT#7,S1
00263INPUT#7,G
00270MAT N=ZER(M)
00280MAT B=ZER(M)
00290MAT L=ZER(M,M)
00300MAT M=ZER(M)
00310MAT A=ZER(M)
00311FOR I=1TOM
00312INPUT#7,N(I)
00313NEXT I
00314FORI=1TOM
00315INPUT#7,B(I)
00316NEXTI
00317FORI=1TOM
00318FORJ=1TOM
00319NEXTJ
00320NEXTI
00321FORI=1TOM
00322INPUT#7,M(I)
00323NEXTI
00324FORI=1TOM
00325INPUT#7,A(I)
00326 NEXT I
00330RESTORE#8
00331INPUT#8,G0
00332INPUT#8,M
00333INPUT#8,S2
00340S7=1
00350MAT E=ZER(M,M)
00352 MAT T=ZER(M,M)
00360MAT F=ZER(M)
00370FORI=1TOM
00371FORJ=1TOM
00372INPUT#8,E(I,J)
00374NEXTJ
00375NEXTI
00376FORI=1TOM
00377INPUT#8,F(I)
00378NEXT I
00380GOTO 890
00390GOSUB 1070
00400PRINT " "
00410PRINT "INPUT THE GROUP NUMBER AND COEFFICIENT YOU WISH FOR THAT GROUP."
00420MAT C=ZER(M)
00430PRINT "GROUP NUMBER,COEFFICIENT(NO MORE = 0,0)";
00440GOSUB 9050
00450IF O1=0 THEN 530
00460IF O1<1 THEN 480
00470IF O1 <= M THEN 500
00480PRINT "ILLEGAL RESPECIFY."
00490GOTO 430
00500I=INT(O1)
00510C(I)=O2
00520GOTO 430
00530FOR I=1 TO M
00540IF C(I) <> 0 THEN 600
00550NEXT I
00560PRINT "ALL COEFFICIENTS = '0'. THIS IS ILLEGAL."
00570PRINT "TO CONTINUE TYPE '1'.";
00580GOSUB 9000
00590CHAIN "CMOD69"
00600M0=0
00610S3=0
00620FOR J=1 TO M
00630M0=M0+C(J)*M(J)
00640IF Q9 <> 1 THEN 690
00650FOR I=1 TO M
00660S3=S3+E(J,I)*C(I)*C(J)
00670NEXT I
00680GOTO 700
00690S3=S3+C(J)*C(J)/N(J)
00700NEXT J
00710PRINT L$
00720PRINT " LINEAR COMBINATION"
00730PRINT "GROUP COEFFICIENT"
00740FOR J=1 TO M
00750PRINT J,C(J)
00760NEXT J
00770PRINT "MEAN=";M0,"STANDARD DEVIATION=";SQR(S1*S3/(G-2))
00780PRINT "DEGREES OF FREEDOM=";G
00790PRINT "TO EXAMINE THIS DISTRIBUTION TYPE '1'."
00800PRINT "ELSE '0'";
00810GOSUB 9000
00820IF O1=1 THEN 840
00830CHAIN "CMOD69"
00840 H$=" DISTRIBUTION OF LINEAR COMBINATION"
00850S0=SQR(S1*S3/(G-2))
00860GOSUB 3500
00870CHAIN "CMOD69"
00880CHAIN "RSTRT"
00890FOR I=1 TO M
00900FOR J=1 TO M
00910E(I,J)=E(I,J)/(S2/G0)
00920NEXT J
00930NEXT I
00940MAT T=INV(E)
00942 MAT E=T
00950GOTO 970
00960MAT T=INV(E)
00962 MAT E=T
00970MAT S=ZER(M)
00980MAT M=ZER(M)
00990MAT N=ZER(M)
01000MAT M=F
01010MAT T=INV(E)
01012 MAT E=T
01020FOR I=1 TO M
01030S(I)=SQR(S2*E(I,I)/(G0-2))
01040NEXT I
01050G=G0
01060GOTO 410
01070RESTORE#2
01071 INPUT# 2,M
01080N0=0
01090S0=0
01100MAT E=ZER(M,M)
01110MAT N=ZER(M)
01120MAT S=ZER(M)
01130MAT M=ZER(M)
01140MAT F=ZER(M)
01150FORI=1TOM
01151INPUT#2,N(I)
01152NEXTI
01153FORI=1TOM
01154INPUT#2,M(I)
01155NEXTI
01156FORI=1TOM
01157INPUT#2,S(I)
01158NEXTI
01160MAT F=M
01170FOR I=1 TO M
01180S0=S0+N(I)*S(I)*S(I)
01190N0=N0+N(I)
01200NEXT I
01210G=N0-M
01220FOR I=1 TO M
01230S(I)=SQR(S0/(G-2)/N(I))
01240E(I,I)=1/N(I)
01250NEXT I
01260S1=S0
01270S2=S1
01280G0=G
01290RETURN
03500REM
03502PRINT L$
03503PRINT "YOU HAVE THE FOLLOWING OPTIONS AVAILABLE FOR EXAMINATION OF"
03506PRINT H$
03509PRINT
03512PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
03515PRINT " 1. HIGHEST DENSITY REGION"
03518PRINT " 2. PROBABILITY LESS THAN SOME VALUE."
03521PRINT " 3. PROBABILITY BETWEEN TWO VALUES."
03524PRINT " 4. EXIT"
03527GOSUB 9000
03530IF O1 <> 4 THEN 3536
03533RETURN
03536IF O1=1 THEN 3551
03539IF O1=2 THEN 3635
03542IF O1=3 THEN 3698
03545PRINT "MUST BE 1,2,3 OR 4. RESPECIFY."
03548GOTO 3527
03551PRINT L$
03554PRINT " HIGHEST DENSITY REGIONS"
03557PRINT " "
03560GOTO 3587
03563PRINT H$
03566PRINT " "
03569:MEAN ######.## DEGREES OF FREEDOM ######
03572PRINT USING 3569,M0,G
03575:STANDARD DEV. ####.### SCALE PARAMETER ########.##
03578PRINT USING 3575,S0,S0*S0*(G-2)
03581PRINT "----------------------------------------------------------------"
03584RETURN
03587GOSUB 3563
03590PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
03593PRINT "INPUT P% AS NUMBER FROM 1 THROUGH 99."
03596PRINT "P%=";
03599GOSUB 9000
03602IF O1=0 THEN 3500
03605IF O1>99 THEN 3614
03608IF O1<1 THEN 3614
03611GOTO 3620
03614PRINT "P% MUST BE AT LEAST 1% AND NOT GREATER THAN 99%. REENTER."
03617GOTO 3599
03620P0=O1/200+.5
03623GOSUB 3810
03626: ##.#% HDR = #####.## TO ######.##
03627F0=J2*SQR(S0*S0*(G-2)/G)
03629PRINT USING 3626,O1,M0-F0,M0+F0
03632GOTO 3596
03635PRINT L$
03638PRINT " PROBABILITY LESS THAN SOME VALUE"
03641PRINT " "
03644GOSUB 3563
03647PRINT "TO EXIT ROUTINE TYPE '7777' WHEN ASKED FOR INPUT."
03650PRINT "INPUT VALUE";
03653GOSUB 9000
03656IF O1=7777 THEN 3500
03659X3=O1
03662Y0=ABS(M0-X3)
03665J2=Y0/S0/SQR((G-2)/G)
03668J1=0
03671GOSUB 6000
03674P=P-.5
03677IF X3<M0 THEN 3686
03680P=.5+P
03683GOTO 3687
03686P=.5-P
03687IF S9=0 THEN 3689
03688RETURN
03689: PROBABILITY LESS THAN ######.## =##.##
03692PRINT USING 3689,O1,P
03695GOTO 3650
03698PRINT L$
03701PRINT " PROBABILITY BETWEEN TWO VALUES"
03704PRINT
03707GOSUB 3563
03710PRINT "TO EXIT ROUTINE TYPE '0,0' WHEN ASKED FOR INPUT."
03713PRINT "INPUT TWO VALUES SEPARATED BY A COMMA. SMALLEST FIRST."
03716GOSUB 9050
03719IF O1 <> 0 THEN 3728
03722IF O2 <> 0 THEN 3728
03725GOTO 3500
03728IF O1 <= O2 THEN 3737
03731PRINT "SMALLER VALUE MUST BE ENTERED FIRST. RESPECIFY."
03734GOTO 3716
03737X3=O1
03740X4=O2
03743Y0=ABS(M0-X3)
03746J2=Y0/S0/SQR((G-2)/G)
03749J1=0
03752GOSUB 6000
03755P=P-.5
03758IF X3<M0 THEN 3767
03761P3=.5+P
03764GOTO 3770
03767P3=.5-P
03770Y0=ABS(M0-X4)
03773J2=Y0/S0/SQR((G-2)/G)
03776J1=0
03779GOSUB 6000
03782P4=P-.5
03785IF X4<M0 THEN 3794
03788P4=.5+P4
03791GOTO 3800
03794P4=.5-P4
03797: PROB(######.## < T < ######.##)=##.##
03800PRINT USING 3797,X3,X4,P4-P3
03803GOTO 3713
03810REM %ILE 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
05000REM BETA CDF
05025 DIM W(16),O(16)
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
05106D9=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
05850REM LOG GAMMA
05853REM INPUT G9
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
06000REM T CDF
06003L0=LOG(1+J2*J2/G)
06004IF L0>.000001 THEN 6006
06005L0=0
06006P=.5
06007IF J2=0 THEN 6035
06008IF J2<6 THEN 6013
06009P=1
06010GOTO 6035
06011IF J2>12 THEN 6009
06013Y3=J2
06015Y3=(G-2/3+.1/G)*SQR(L0/(G-5/6))
06016GOSUB 8000
06017GOTO 6035
06035RETURN
08000REM NORMAL CDF
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
09000REM
09005INPUT O1
09010IF O1=-9999 THEN 9020
09015RETURN
09020CHAIN "RSTRT"
09050INPUT O1,O2
09055IF O2=-9999 THEN 9020
09060GOTO 9010
09500L0=0
09505L1=1
09510FOR L3=1 TO R0
09520L4=L3
09525L5=L3+1
09530FOR L6=L5 TO R0
09535IF ABS(L(L4,L3)) >= ABS(L(L6,L3)) THEN 9545
09540L4=L6
09545NEXT L6
09550IF ABS(L(L4,L3))>10 ** (-7) THEN 9560
09552IF R0<C0 THEN 9655
09555RETURN
09560FOR L6=L3 TO C0
09565L9=L(L3,L6)
09570L(L3,L6)=L(L4,L6)
09575L(L4,L6)=L9
09580NEXT L6
09585L9=L(L3,L3)
09590L(L3,L3)=1
09595FOR L6=L5 TO C0
09600L(L3,L6)=L(L3,L6)/L9
09605NEXT L6
09610FOR L6=1 TO R0
09615IF L6=L3 THEN 9640
09620L9=L(L6,L3)
09625FOR L8=L3 TO C0
09630L(L6,L8)=L(L6,L8)-L9*L(L3,L8)
09635NEXT L8
09640NEXT L6
09645NEXT L3
09650IF R0 >= C0 THEN 9675
09655FOR L1=1 TO C0
09660IF L(R0,L1) <> 0 THEN 9675
09665NEXT L1
09670RETURN
09675L0=1
09680RETURN
09999 END