Trailing-Edge
-
PDP-10 Archives
-
decuslib20-03
-
decus/20-0096/tbf.bas
There is 1 other file named tbf.bas in the archive. Click here to see a list.
01010 REM "TBF" DENSITY & CUMUL FOR -T- OR BEHRENS DISTRIBUTION
01020 REM
01030 REM WRITTEN BY
01040 REM JAMES FENNESSEY AND SUSAN RADIUS
01050 REM DEPARTMENT OF SOCIAL RELATIONS
01060 REM THE JOHNS HOPKINS UNIVERSITY
01070 REM BALTIMORE, MARYLAND 21218
01080 REM
01090 REM 17 NOVEMBER 1975 0800 VERSION
01100 REM
01110 REM THIS PROGRAM HAS BEEN CHECKED WITH REASONABLE CARE,
01120 REM BUT IS NOT GUARANTEED. PLEASE NOTIFY THE AUTHORS
01130 REM OF ANY ERRORS.
01140 REM
01150 FILE#1,"TBFOUT.BAS"
01160 SCRATCH #1
01170 MARGIN 80 '80 CHARACTER PRINT LINE USED
01180 MARGIN#1,80 '80 CHARACTER PRINT LINE USED
01190 PAGE (60)
01200 PAGE#1,(60)
01210 PRINT <PA>
01220 PRINT#1<PA>
01230 PRINT "'TBF' DISTRIBUTION EVALUATOR"
01240 PRINT#1,"'TBF' DISTRIBUTION EVALUATOR"
01250 PRINT "TYPE 2 FOR PLOT ONLY; OTHERWISE TYPE 1"
01260 INPUT S9
01270 GOSUB 3450'PRINTS A DASHED LINE
01280 PRINT "TYPE LABEL FOR PROBLEM (TO STOP, TYPE 'END')"
01290 INPUT A$
01300 IF A$="END" THEN 2680
01310 PRINT A$
01320 PRINT#1,A$
01330 PRINT "TYPE 1 IF SINGLE -T-, 2 IF DUAL -T-, 3 IF BEHRENS"
01340 INPUT S8
01350 ON S8 GO TO 1360,1470,1470
01360 PRINT "TYPE N, MEAN, VARIANCE"
01370 INPUT N,M,V
01380 D=N-1 'D=DEGREES OF FREEDOM
01390 S2=V/N 'S2=SQUARED SCALE FACTOR OF -T- DIST.
01400 L$="SINGLE -T-"
01410 PRINT L$," N"," M"," V"
01420 PRINT#1,L$," N"," M"," V"
01430 PRINT ,N,M,V
01440 PRINT#1, ,N,M,V
01450 GO TO 1910
01460 REM
01470 PRINT "TYPE N, MEAN, VARIANCE OF GROUP 1"
01480 INPUT N1,M1,V1
01490 PRINT "TYPE N, MEAN, VARIANCE OF GROUP 2"
01500 INPUT N2,M2,V2
01510 IF S8=3 THEN 1640
01520 L$="DUAL -T-"
01530 PRINT L$," N"," M"," V"
01540 PRINT#1,L$," N"," M"," V"
01550 PRINT "GRP 1",N1,M1,V1
01560 PRINT#1,"GRP 1",N1,M1,V1
01570 PRINT "GRP 2",N2,M2,V2
01580 PRINT#1,"GRP 2",N2,M2,V2
01590 D=N1+N2-2 'D=DEGREES OF FREEDOM
01600 M=M2-M1 'M=MEAN OF -T- DIST FOR DIFFERENCES
01610 S2=(1/N1+1/N2)*(((N1-1)*V1+(N2-1)*V2)/D) 'S2=SQUARED SCALE FACTOR
01620 GO TO 1910
01630 REM
01640 REM SEE G.E.P. BOX & G.C. TIAO, BAYESIAN INFERENCE IN
01650 REM STATISTICAL ANALYSIS (ADDISON-WESLEY, 1973) P. 107
01660 C9=(V2/N2)/(V1/N1+V2/N2) 'C9=(COS(PHI))^2
01670 C8=1-C9 'C8=(SIN(PHI))^2
01680 C7=(((N2-1)/((N2-1)-2))*C9)+(((N1-1)/((N1-1)-2))*C8)
01690 C6=(((N2-1)^2)/(((N2-1)-2)^2*((N2-1)-4)))*C9^2
01700 C5=(((N1-1)^2)/(((N1-1)-2)^2*((N1-1)-4)))*C8^2
01710 C4=C6+C5
01720 C3=(V1/N1)+(V2/N2)
01730 C2=4+(C7^2)/C4
01740 C1=((C2-2)/C2)*C7
01750 C1=C1*C3
01760 D=C2 'D=DEGREES OF FREEDOM IN -T- APPROXIMATION
01770 M=M2-M1 'M=MEAN OF BEHRENS DIST
01780 S2=C1 'S2=SQUARED SCALE FACTOR FOR -T- APPROXIMATION
01790 S7=SQR((V1/N1)/(V2/N2))
01800 S6=ATN(S7) 'S6=ANGLE PHI OF BEHRENS DISTRIBUTION (IN RADIANS)
01810 PRINT "PHI=";180/3.14159265*S6;" TAN(PHI)=";S7
01820 PRINT#1,"PHI=";180/3.14159265*S6;" TAN(PHI)=";S7
01830 L$="BEHRENS"
01840 PRINT L$," N"," M"," V"
01850 PRINT#1,L$," N"," M"," V"
01860 PRINT "GRP 1",N1,M1,V1
01870 PRINT#1,"GRP 1",N1,M1,V1
01880 PRINT "GRP 2",N2,M2,V2
01890 PRINT#1,"GRP 2",N2,M2,V2
01900 REM END OF BEHRENS PREPARATION
01910 S=SQR(S2) 'S=-T- SCALE FACTOR
01920 F9=INT(S*10^4+.5)/10^4
01930 D9=INT(D*10^4+.5)/10^4
01940 M9=INT(M*10^4+.5)/10^4
01950 PRINT
01960 PRINT#1
01970 IF S8<>3 THEN 2000
01980 PRINT " FOR -T- DIST USED AS APPROXIMATION,"
01990 PRINT#1," FOR -T- DIST USED AS APPROXIMATION,"
02000 PRINT "D.F. =";D9;" MEAN =";M9;" SCALE FACTOR =";F9
02010 PRINT#1,"D.F. =";D9;" MEAN =";M9;" SCALE FACTOR =";F9
02020 REM CALCULATE K (= LN(CONSTANT) IN -T- DISTRIBUTION)
02030 G9=(D+1)/2
02040 GOSUB 2710
02050 S5=G0
02060 G9=D/2
02070 GOSUB 2710
02080 S4=G0
02090 K=S5-(S4+LN(S*SQR(3.141592654*D))) 'K NOW READY
02100 L=M-(3.5*F9)
02110 U=L+2*(M-L)
02120 X=L
02130 GOSUB 3550'SUBROUTINE SETS C=CUMUL PROB AT X
02140 P=FNT(X) 'FUNCTION T SETS P=DENSITY AT X
02150 PRINT "L = ";L,"U = ";U
02160 PRINT USING 2170,L,P,C
02170 :AT L = ###.#### DENSITY =##.#### CUMUL PROB =##.####
02180 PRINT "TYPE NEW LOWER BOUND OR ELSE TYPE ";M9
02190 INPUT L1
02200 IF ABS(L1-M9)<.00005 THEN 2280
02210 IF L1>M THEN 2260
02220 PRINT "TYPE NEW UPPER BOUND"
02230 INPUT U
02240 L=L1
02250 GO TO 2120
02260 PRINT "L MUST BE < M"
02270 GO TO 2180
02280 PRINT "L = ";L,"U = ";U
02290 PRINT#1,"L = ";L,"U = ";U
02300 PRINT "WIDTH OF AN INTERVAL"
02310 INPUT W
02320 W1=INT((U-L)/W)+1
02330 IF W1<=200 THEN 2360
02340 PRINT "TOO MANY POINTS"
02350 GO TO 2300
02360 PRINT "VALUES FOR ";W1;" POINTS."
02370 PRINT#1,"VALUES FOR ";W1;" POINTS."
02380 IF S9<>2 THEN 2410
02390 GOSUB 3030'SUBROUTINE PLOTS THE DENSITY
02400 GOTO 1210'GO AND BEGIN NEW PROBLEM
02410 PRINT
02420 PRINT#1
02430 PRINT " X-AXIS"," ";"DENSITY"," CUMUL. PROB."
02440 PRINT#1," X-AXIS"," ";"DENSITY"," CUMUL. PROB."
02450 PRINT
02460 PRINT#1
02470 FOR J=1 TO W1 'BEGIN MAIN LOOP OF TABULATION
02480 X=L+(J-1)*W
02490 GOSUB 3550
02500 P=FNT(X)
02510 PRINT USING 2530,J,X,P,C
02520 PRINT#1,USING 2530,J,X,P,C
02530 :### ####.#### ##.#### ##.####
02540 NEXT J
02550 PRINT
02560 PRINT#1
02570 PRINT "END OF TABULATION"
02580 PRINT#1,"END OF TABULATION"
02590 GOSUB 3450
02600 PRINT "DO YOU WANT THE PLOT"
02610 INPUT B$
02620 IF LEFT$(B$,1)="Y" THEN 2640
02630 GO TO 1210
02640 GOSUB 2960'SUBROUTINE PLOTS THE DENSITY
02650 PRINT
02660 PRINT#1
02670 GO TO 1210
02680 STOP
02690 REM END OF MAIN PROGRAM; SUBROUTINES FOLLOW
02700 REM ***************************************
02710 REM THIS SUBROUTINE SETS G0=LN(GAMMA(G9)), USING STIRLING FORMULA
02720 IF G9<=1E30 THEN 2750
02730 G0=1E38
02740 RETURN
02750 IF G9>1E-9 THEN 2780
02760 G0=0
02770 RETURN
02780 IF G9<1E10 THEN 2810
02790 G0=G9*(LOG(G9)-1)
02800 RETURN
02810 G8=G9 'HANDLE INPUT > 10^-9 AND < 10^10
02820 G7=1
02830 IF G8>18 THEN 2870
02840 G7=G7*G8 'ADJUST VIA G7 IF INPUT < 18
02850 G8=G8+1
02860 GO TO 2830
02870 G6=1/G8^2 'SET UP VALUES FOR STIRLING ASYMPTOTIC FORMULA
02880 G5=1/12
02890 G4=1/360
02900 G3=1/1260
02910 G2=1/1680
02920 G1=LOG(SQR(2*3.141592654))
02930 G0=G1+(G8-.5)*LOG(G8)-G8+1/G8*(G5-(G6*(G4+(G6*(G3-G6*G2)))))-LOG(G7)
02940 RETURN
02950 REM
02960 REM THIS SUBROUTINE PLOTS Y VERSUS X
02970 PRINT "TYPE LOWER AND UPPER BOUNDS AND STEP WIDTH"
02980 INPUT L,U,W
02990 W1=INT((U-L)/W)+1
03000 IF W1<=200 THEN 3030
03010 PRINT "TOO MANY POINTS"
03020 GO TO 2970
03030 X=M
03040 P=FNT(X)
03050 Y9=P/50 'Y9 IS DENSITY VALUE = ONE PRINTER COLUMN
03060 PRINT <PA>
03070 PRINT#1,<PA>
03080 GOSUB 3450
03090 PRINT A$
03100 PRINT#1,A$
03110 PRINT L$;" DENSITY ON ";D9;" D.F., MEAN =";M9;" SCALE FACTOR =";F9
03120 PRINT#1,L$;" DENSITY ON ";D9;" D.F., MEAN =";M9;" SCALE FACTOR =";F9
03130 Y8=INT(Y9*10^4+.5)/10^4
03140 PRINT "Y-AXIS IS * IF DENSITY <";Y8;"; OTHERWISE IS I"
03150 PRINT#1,"Y-AXIS IS * IF DENSITY <";Y8;"; OTHERWISE IS I"
03160 PRINT " TICK MARKS ON Y-AXIS EVERY";5*Y8
03170 PRINT#1," TICK MARKS ON Y-AXIS EVERY";5*Y8
03180 PRINT
03190 PRINT#1
03200 PRINT TAB(16);"0";TAB(64);50*Y8
03210 PRINT#1,TAB(16);"0";TAB(64);50*Y8
03220 PRINT TAB(16);"I----I----I----I----I----I----I----I----I----I----I----I"
03230 PRINT#1,TAB(16);"I----I----I----I----I----I----I----I----I----I----I----I"
03240 FOR J=1 TO W1 'BEGIN MAIN LOOP OF PLOT
03250 B$="*"
03260 C$="I"
03270 X=L+(J-1)*W
03280 P=FNT(X)
03290 IF P>=Y8 THEN 3320
03300 B$=" "
03310 C$="*"
03320 Y=INT(P/Y9+.5)
03330 X=INT(X*10^4+.5)/10^4
03340 PRINT J;TAB(6);X;TAB(16);C$;TAB(16+Y);B$
03350 PRINT#1,J;TAB(6);X;TAB(16);C$;TAB(16+Y);B$
03360 NEXT J
03370 PRINT
03380 PRINT#1
03390 GOSUB 3450
03400 PRINT
03410 PRINT#1
03420 RETURN
03430 REM
03440 REM THIS SUBROUTINE PRINTS A DASHED LINE
03450 FOR I=1 TO 39
03460 PRINT "-";
03470 PRINT#1,"-";
03480 NEXT I
03490 PRINT "-"
03500 PRINT#1,"-"
03510 RETURN
03520 REM
03530 DEF FNT(X)=EXP(K+LN((1+(((X-M)^2)*(D+1))/(D*S2*(D+1)))^(-(D+1)/2)))
03540 REM
03550 REM THIS SUBROUTINE GETS T CUMUL PROB USING BETA DISTRIBUTION
03560 REM THIS SUBROUTINE WAS ADAPTED WITH PERMISSION FROM THE MANECON
03570 REM PACKAGE OF PROGRAMS DEVELOPED AT THE HARVARD BUSINESS SCHOOL,
03580 REM AND DESCRIBED IN THE BOOK, "COMPUTER PROGRAMS FOR ELEMENTARY
03590 REM DECISION ANALYSIS", BY ROBERT SCHLAIFER (1971).
03600 REM THE MANECON ROUTINES USED AS THE BASIS FOR THIS ROUTINE ARE:
03610 REM "STUCUM","DBETCU","SPCASE"
03620 REM INPUTS X,M,S,D; OUTPUT C
03630 H1=1
03640 T9=.0001 'TOLERANCE VALUE FOR BETA EVALUATION
03650 T8=-88
03660 T7=+88
03670 T6=100 'MAX TERMS OF SOLUTION, REGARDLESS OF T9
03680 T5=.5*D
03690 A0=-1
03700 T=(X-M)/S
03710 Z2=.5*(1+T/SQR(D+T*T))
03720 Z=Z2
03730 IF Z2<.5 THEN 3750
03740 Z=1-Z2
03750 REM
03760 IF Z<0 THEN 4380
03770 IF Z>1 THEN 4380
03780 IF T5<0 THEN 4380
03790 IF A0=T5 THEN 3890
03800 G9=T5
03810 GOSUB 2710
03820 D3=G0
03830 G9=D
03840 GOSUB 2710
03850 D4=G0
03860 A1=2*D3-D4
03870 A2=LN(D)
03880 A0=T5
03890 A3=LN(Z)
03900 A4=LN(1-Z)
03910 A5=(T5-1)*(A3+A4)-A1
03920 IF A5<T8 THEN 4380
03930 IF A5>T7 THEN 4380
03940 H2=EXP(A5)
03950 IF H1<>0 THEN 3970
03960 RETURN
03970 REM INITIALIZE FOR CUMUL
03980 C=0
03990 A6=Z/(1-Z)
04000 A7=H2*(1-Z)/T5
04010 REM EVALUATE BINOMIAL SERIES
04020 E1=INT(T5)
04030 IF E1<1 THEN 4180
04040 E2=T6
04050 IF E1>T6 THEN 4070
04060 E2=E1
04070 A8=0
04080 FOR I=1 TO E2
04090 E3=(T5-A8)/(T5+A8)*A6
04100 IF A7>T9 THEN 4120
04110 IF A7<(T9/E3-T9) THEN 4330
04120 A7=E3*A7
04130 C=C+A7
04140 A8=A8+1
04150 NEXT I
04160 IF E2<E1 THEN 4310
04170 REM EVALUATE POWER SERIES
04180 E4=E1
04190 IF E4=T5 THEN 4330
04200 A7=A7*(T5-E4)/(T5+E4)*Z/(1-Z)^(T5-E4)
04210 C=C+A7
04220 E5=T9/A6
04230 A8=E4+1
04240 E6=E1+1
04250 FOR I=E6 TO T6
04260 A7=(T5+A8-1)/(T5+A8)*(A8-T5)/(A8-E4)*Z*A7
04270 C=C+A7
04280 IF A7<E5 THEN 4330
04290 A8=A8+1
04300 NEXT I
04310 PRINT "ITERATION DID NOT CONVERGE; CUM PROB MAY BE INACCURATE"
04320 PRINT#1,"ITERATION DID NOT CONVERGE; CUM PROB MAY BE INACCURATE"
04330 IF Z2<.5 THEN 4350
04340 C=1-C
04350 RETURN
04360 REM THIS IS NORMAL RETURN; C = CUMUL PROB
04370 REM HANDLE ILLEGAL AND SPECIAL CASES
04380 PRINT "SPECIAL CASE ROUTINE IN USE"
04390 PRINT#1,"SPECIAL CASE ROUTINE IN USE"
04400 A5=0
04410 Z2=Z
04420 IF Z<.5 THEN 4460
04430 Z2=1-Z
04440 GOSUB 4490
04450 RETURN
04460 GOSUB 4490
04470 C=1-C
04480 RETURN
04490 REM SPECIAL CASES SUBROUTINE
04500 C=0
04510 H2=0
04520 IF A5<0 THEN 4550
04530 IF A5=0 THEN 4580
04540 IF A5>0 THEN 4590
04550 IF Z2>T5 THEN 4570
04560 RETURN
04570 C=1
04580 RETURN
04590 PRINT "ARGUMENTS BEYOND CAPACITY; SPURIOUS PROBABILITY"
04600 PRINT#1,"ARGUMENTS BEYOND CAPACITY; SPURIOUS PROBABILITY"
04610 H2=2E22
04620 RETURN
04630 END