Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod88.bas
There are 2 other files named cmod88.bas in the archive. Click here to see a list.
00030 A$="-"
00040Q9=0
00050REM******************************************************************
00060REM CMOD88 CMOD88 CMOD88 CMOD88 CMOD88
00070REM INTERMEDIATE MODULE FOR ANOVA PAIRWISE COMPARISONS
00080REM******************************************************************
00090 DIM M(12),N(12),S(12),A(12),C(12) ,P(6,6),Q(6),U(1,1),G(1,12)
00100 DIM L(12,12),I(12,12),V(12),T(12,12),B(12),E(12,12),F(12)
00110 C$="11122122"
00120GOSUB 5160
00130 FILES RFILE1,RFILE2,RFILE3,,,,RF7,RF8,RF9
00170RESTORE#1
00171 INPUT# 1,I1,I2,I3
00180SCRATCH#1
00181 PRINT # 1,88,I2,I3
00190I2=I3
00200S7=0
00210X=0
00220S9=0
00230IF I1=75 THEN 280
00240RESTORE#8
00241 INPUT# 8 ,G
00250IF G=0 THEN 420
00260I2=0
00270Q9=1
00280RESTORE#7
00281INPUT#7,M
00282INPUT#7,S1
00283INPUT#7,G
00290MAT N=ZER(M)
00300MAT B=ZER(M)
00310MAT L=ZER(M,M)
00320MAT M=ZER(M)
00330MAT A=ZER(M)
00331FOR I=1TOM
00332INPUT#7,N(I)
00333NEXTI
00334FORI=1TOM
00335INPUT#7,B(I)
00336NEXTI
00337FORI=1TOM
00338FORJ=1TOM
00339INPUT#7,L(I,J)
00340NEXTJ
00341NEXTI
00342FOR I=1TOM
00343INPUT#7,M(I)
00344NEXTI
00345FORI=1TOM
00346INPUT#7,A(I)
00347NEXTI
00350RESTORE#8
00351INPUT#8,G0
00352INPUT#8,M
00353INPUT#8,S2
00360S7=1
00370MAT E=ZER(M,M)
00372MAT T=ZER(M,M)
00380MAT F=ZER(M)
00390 FORI=1TOM
00391FORJ=1TOM
00392INPUT#8,E(I,J)
00393NEXTJ
00394NEXTI
00395FORI=1TOM
00396INPUT#8,F(I)
00397NEXT I
00400GOTO 1400
00410CHAIN "RSTRT"
00420GOSUB 1570
00430PRINT L$
00440MAT P=ZER
00450PRINT " OPTION 3. PAIR-WISE COMPARISONS"
00460PRINT " "
00470PRINT "THIS OPTION PRESENTS PROBABILITIES THAT THE DIFFERENCES"
00480PRINT "BETWEEN SPECIFIED PARAMETERS ARE GREATER THAN A VALUE X0."
00490PRINT "INPUT VALUE X0.";
00500GOSUB 9000
00510X3=O1
00520IF M>6 THEN 580
00530M1=M
00540FOR I=1 TO M
00550Q(I)=I
00560NEXT I
00570GOTO 780
00580PRINT "INPUT THE NUMBER OF PARAMETERS FOR WHICH YOU WISH TO SEE"
00590PRINT "PROBABILITIES ON DIFFERENCES(MAX=6)"
00600GOSUB 9000
00610O1=INT(O1)
00620IF O1 >= 2 THEN 650
00630PRINT "MUST BE BETWEEN 2 AND 6. RESPECIFY."
00640GOTO 600
00650IF O1>6 THEN 630
00660M1=O1
00670MAT Q=ZER
00680FOR I=1 TO M1
00690PRINT "INPUT PARAMETER NUMBER FOR WHICH YOU WISH TO SEE DIFFERENCES.";
00700GOSUB 9000
00710O1=INT(O1)
00720IF O1 <= M THEN 750
00730PRINT "INVALID. RESPECIFY.";
00740GOTO 700
00750IF O1<1 THEN 730
00760Q(I)=O1
00770NEXT I
00780FOR I=1 TO M1
00790FOR J=1 TO M1
00800IF I <> J THEN 830
00810P=0
00820GOTO 980
00830M0=M(Q(I))-M(Q(J))
00840IF Q9 <> 1 THEN 870
00850S3=E(Q(J),Q(J))+E(Q(I),Q(I))-2*E(Q(I),Q(J))
00860GOTO 880
00870S3=1/N(Q(I))+1/N(Q(J))
00880S0=SQR(S1*S3/(G-2))
00890Y0=ABS(M0-X3)
00900J2=Y0/S0/SQR((G-2)/G)
00910J1=0
00920GOSUB 6000
00930P=P-.5
00940IF X3<M0 THEN 970
00950P=.5+P
00960GOTO 980
00970P=.5-P
00980P(I,J)=1-P
00982IF P(I,J) <= .99 THEN 990
00984P(I,J)=.99
00990NEXT J
01000NEXT I
01010PRINT L$
01020PRINT "PROBABILITIES THAT DIFFERENCES ARE GREATER THAN";X3
01030PRINT " "
01040PRINT "--------------------------------------------------------------"
01050:PARAMETER ## ## ## ## ## ##
01060:PARAMETER ## ## ## ## ##
01070:PARAMETER ### ### ### ###
01080:PARAMETER ## ## ##
01090:PARAMETER ## ##
01100 ONM1-1 GOTO 1110,1130,1150,1170,1190
01110PRINT USING 1090,Q(1),Q(2)
01120GOTO 1200
01130PRINT USING 1080,Q(1),Q(2),Q(3)
01140GOTO 1200
01150PRINT USING 1070,Q(1),Q(2),Q(3),Q(4)
01160GOTO 1200
01170PRINT USING 1060,Q(1),Q(2),Q(3),Q(4),Q(5)
01180GOTO 1200
01190PRINT USING 1050,Q(1),Q(2),Q(3),Q(4),Q(5),Q(6)
01200FOR I=1 TO M1
01210:" ## "
01211SCRATCH#9
01220PRINT#9USING 1210,Q(I)
01222RESTORE#9
01223INPUT#9,D1$
01224 PRINT D1$;
01230FOR J=1 TO M1
01240IF I <> J THEN 1280
01250:" 'L "
01251SCRATCH#9
01260PRINT#9USING 1250,A$
01262RESTORE#9
01264 INPUT#9,D1$
01265PRINTD1$;
01270GOTO 1300
01280:"##.## "
01281SCRATCH#9
01290PRINT#9USING 1280,P(I,J)
01292RESTORE#9
01294INPUT#9,D1$
01296 PRINT D1$;
01300NEXT J
01310PRINT " "
01320NEXT I
01330PRINT "--------------------------------------------------------------"
01340PRINT "TO EXAMINE ANOTHER VALUE TYPE '1'."
01350PRINT "ELSE '0'";
01360GOSUB 9000
01370IF O1 <> 0 THEN 1390
01380CHAIN "CMOD69"
01390GOTO 430
01400FOR I=1 TO M
01410FOR J=1 TO M
01420E(I,J)=E(I,J)/(S2/G0)
01430NEXT J
01440NEXT I
01450MAT T=INV(E)
01452 MAT E=T
01460GOTO 1470
01470MAT S=ZER(M)
01480MAT M=ZER(M)
01490MAT N=ZER(M)
01500MAT M=F
01510MAT T=INV(E)
01512 MAT E=T
01520FOR I=1 TO M
01530S(I)=SQR(S2*E(I,I)/(G0-2))
01540NEXT I
01550G=G0
01560GOTO 430
01570RESTORE#2
01571 INPUT# 2,M
01580N0=0
01590S0=0
01600MAT E=ZER(M,M)
01610MAT N=ZER(M)
01620MAT S=ZER(M)
01630MAT M=ZER(M)
01640MAT F=ZER(M)
01650FORI=1TOM
01651INPUT#2,N(I)
01652NEXT I
01653FORI=1TOM
01654NEXTI
01655FORI=1TOM
01656INPUT#2,S(I)
01657NEXT I
01660MAT F=M
01670FOR I=1 TO M
01680S0=S0+N(I)*S(I)*S(I)
01690N0=N0+N(I)
01700NEXT I
01710G=N0-M
01720FOR I=1 TO M
01730S(I)=SQR(S0/(G-2)/N(I))
01740E(I,I)=1/N(I)
01750NEXT I
01760S1=S0
01770S2=S1
01780G0=G
01790RETURN
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
09999END