Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod86.bas
There are 2 other files named cmod86.bas in the archive. Click here to see a list.
00020Q9=0
00030REM**************************************************************
00040REM CMOD86 CMOD86 CMOD86 CMOD86 CMOD86
00050REM INTERMEDIATE MODULE ANOVA MARGINAL DISTRIBUTIONS.
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,86,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)
00311FORI=1TOM
00312INPUT#7,N(I)
00313NEXTI
00314FORI=1TOM
00315INPUT#7,B(I)
00316NEXTI
00317FORI=1TOM
00318 FORJ=1TOM
00319INPUT#7,L(I,J)
00320NEXTJ
00321NEXTI
00322FORI=1TOM
00323INPUT#7,M(I)
00324NEXTI
00325FORI=1TOM
00326INPUT#7,A(I)
00327NEXTI
00330RESTORE#8
00331INPUT#8,G0
00332INPUT#8,M
00333INPUT#8,S2
00340S7=1
00350MAT E=ZER(M,M)
00352MAT T=ZER(M,M)
00360MAT F=ZER(M)
00370FORI=1TOM
00371FORJ=1TOM
00372INPUT#8,E(I,J)
00373NEXTJ
00374NEXTI
00375FORI=1TOM
00376INPUT#8,F(I)
00377NEXTI
00380GOTO 1060
00390GOSUB 1240
00400PRINT
00410PRINT "INPUT THE NUMBER OF THE MARGINAL YOU WANT TO EXAMINE OR '0'.";
00420GOSUB 9000
00430O1=INT(O1)
00440IF O1 <= M THEN 480
00450PRINT "NUMBER MUST BE 0 OR NUMBER FROM 1 TO ";M
00460PRINT "RESPECIFY."
00470GOTO 420
00480IF O1 <> 0 THEN 500
00490CHAIN "CMOD69"
00500IF O1<1 THEN 450
00510I=O1
00520M0=M(I)
00530S0=S(I)
00540 H$="THE MARGINAL DISTRIBUTION FOR PARAMETER "
00550 H$=H$+STR$(I)
00560 GOSUB 3500
00570CHAIN "CMOD69"
00580PRINT "INPUT THE GROUP NUMBER AND COEFFICIENT YOU WISH FOR THAT GROUP."
00590MAT C=ZER(M)
00600PRINT "GROUP NUMBER,COEFFICIENT(NO MORE = 0,0)";
00610GOSUB 9050
00620IF O1=0 THEN 700
00630IF O1<1 THEN 650
00640IF O1 <= M THEN 670
00650PRINT "ILLEGAL RESPECIFY."
00660GOTO 600
00670I=INT(O1)
00680C(I)=O2
00690GOTO 600
00700FOR I=1 TO M
00710IF C(I) <> 0 THEN 770
00720NEXT I
00730PRINT "ALL COEFFICIENTS = '0'. THIS IS ILLEGAL."
00740PRINT "TO CONTINUE TYPE '1'.";
00750GOSUB 9000
00760CHAIN "CMOD69"
00770M0=0
00780S3=0
00790FOR J=1 TO M
00800M0=M0+C(J)*M(J)
00810IF Q9 <> 1 THEN 860
00820FOR I=1 TO M
00830S3=S3+E(J,I)*C(I)*C(J)
00840NEXT I
00850GOTO 870
00860S3=S3+C(J)*C(J)/N(J)
00870NEXT J
00880PRINT L$
00890PRINT " LINEAR COMBINATION"
00900PRINT "GROUP COEFFICIENT"
00910FOR J=1 TO M
00920PRINT J,C(J)
00930NEXT J
00940PRINT "MEAN=";M0,"STANDARD DEVIATION=";SQR(S1*S3/(G-2))
00950PRINT "DEGREES OF FREEDOM=";G
00960PRINT "TO EXAMINE THIS DISTRIBUTION TYPE '1'."
00970PRINT "ELSE '0'";
00980GOSUB 9000
00990IF O1=1 THEN 1010
01000CHAIN "CMOD69"
01010 H$=" DISTRIBUTION OF LINEAR COMBINATION"
01020S0=SQR(S1*S3/(G-2))
01030GOSUB 3500
01040CHAIN "CMOD69"
01050CHAIN "RSTRT"
01060FOR I=1 TO M
01070FOR J=1 TO M
01080E(I,J)=E(I,J)/(S2/G0)
01090NEXT J
01100NEXT I
01110MAT T=INV(E)
01112MAT E=T
01120GOTO 1140
01130MAT T=INV(E)
01132 MATE=T
01140MAT S=ZER(M)
01150MAT M=ZER(M)
01160MAT N=ZER(M)
01170MAT M=F
01180MATE=T
01190FOR I=1 TO M
01200S(I)=SQR(S2*E(I,I)/(G0-2))
01210NEXT I
01220G=G0
01230GOTO 400
01240RESTORE#2
01241 INPUT# 2,M
01250N0=0
01260S0=0
01270MAT E=ZER(M,M)
01280MAT N=ZER(M)
01290MAT S=ZER(M)
01300MAT M=ZER(M)
01310MAT F=ZER(M)
01320FORI=1TOM
01321INPUT#2,N(I)
01322NEXTI
01323FORI=1TOM
01324INPUT#2,M(I)
01325NEXTI
01326FORI=1TOM
01327INPUT#2,S(I)
01328NEXTI
01330MAT F=M
01340FOR I=1 TO M
01350S0=S0+N(I)*S(I)*S(I)
01360N0=N0+N(I)
01370NEXT I
01380G=N0-M
01390FOR I=1 TO M
01400S(I)=SQR(S0/(G-2)/N(I))
01410E(I,I)=1/N(I)
01420NEXT I
01430S1=S0
01440S2=S1
01450G0=G
01460RETURN
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