100'  NAME--STAT18
110'
120'  DESCRIPTION--COMPUTES THE ANALYSIS OF VARIANCE TABLE AND
130'  THE F-RATIO FOR TREATMENTS FOR A YOUDEN SQUARE DESIGN.
140'  SUM-OF-SQUARES FOR TREATMENTS IS ADJUSTED BECAUSE OF INCOMPLETENESS.
150'
160'  SOURCE--UNKNOWN
170'
180'  INSTRUCTIONS--PLACE DATA IN LINE 1200 AND FOLLOWING.
190'  FIRST DATA IS N,THE NUMBER OF ROWS AND TREATMENTS, THEN K,
200'  THE NUMBER OF COLUMNS AND REPLICATIONS OF EACH TREATMENT
210'  IN THE EXPERIMENT. NEXT ENTER THE M(I,H) MATRIX BY ROWS.
220'  M(I,H)=J IF TREATMENT J APPEARS IN ROW I,COLUMN H.
230'  THIS MATRIX HAS DIMENSIONS N BY K. NEXT ENTER THE N(I,J)
240'  MATRIX. N(I,J)=1 IF TREATMENT J APPEARS IN ROW I AND IS
250'  0 OTHERWISE.  THE N(I,J) MATRIX HAS DIMENSIONS N BY N.
260'  FINALLY ENTER THE MATRIX OF OBSERVATIONS X(I,H) BY ROWS.
270'  IF N>10 DIM STATEMENTS MUST BE ADDED. 
280'  SAMPLE DATA ARE IN LINES 1200 THROUGH 1320.
290'
300'
310'  *  *  *  *  *  *   MAIN PROGRAM   *  *  *  *  *  *  *  *  *
320'
330 READ N, K
340 LET L = K*(K-1)/(N-1)
350 MAT READ M(N,K)
360 MAT READ N(N,N)
370 FOR I = 1 TO N
380    FOR H = 1 TO K
390       READ X(I,H)
400       LET R(I) = R(I) + X(I,H)
410       LET C(H) = C(H) + X(I,H)
420       LET T(M(I,H)) = T(M(I,H)) + X(I,H)
430       LET S = S + X(I,H)
440       LET S2 = S2 + X(I,H)^2
450    NEXT H
460 NEXT I
470 FOR J = 1 TO N
480    FOR I = 1 TO N
490       LET P(J) = P(J) + N(I,J)*R(I)
500    NEXT I
510 LET Q(J)=K*T(J)-P(J)
520 NEXT J
530 LET C = S*S/N/K
540 FOR I = 1 TO N
550    LET R1 = R1 + R(I)^2
560    LET T1 = T1 + Q(I)^2
570 NEXT I
580 FOR H = 1 TO K
590    LET C1 = C1 + C(H)^2
600 NEXT H
610 LET R2 = R1/K - C
620 LET C2 = C1/N - C
630 LET T2 = T1/N/K/L
640 LET D = (N-1)*(K-2)
650 LET D1 = N*(K-2) + 1
660 LET E2 = S2 - R2 - C2 - T2 - C
670 PRINT "ANOVA TABLE:"
680 PRINT
690 PRINT "ITEM", "SS", "DF", "MS"
700 PRINT
710 PRINT "GRAND TOTAL", S2, N*K
720 PRINT "GRAND MEAN", C, " 1"
730 PRINT "TREATMENTS", T2, N-1, T2/(N-1)
740 PRINT "ROWS", R2, N-1"     ...SS NOT ADJUSTED..."
750 PRINT "COLUMNS", C2, K-1, C2/(K-1)
760 PRINT "ERROR", E2, D, E2/D
770 PRINT
780 PRINT
790 LET F = T2/(N-1)/(E2/D)
800 PRINT "TREATMENT F-RATIO ="F", ON"N-1"AND"D"DEGREES OF FREEDOM."
810 LET G=F
820 LET M=N-1
830 LET N=D
840 GOSUB 900
850 PRINT
860 PRINT "IF MSC/MSE ="C2/(K-1)/(E2/D)" IS NOT SIGNIFICANT, IT MAY BE"
870 PRINT "DESIRABLE TO POOL COLUMN AND ERROR SS TO OBTAIN AS AN"
880 PRINT"ERROR MS ESTIMATE";(C2+E2)/D1;"WITH";D1;"DEGREES OF FREEDOM."
890 STOP
900 REM THE SUBROUTINE FOR COMPUTATION OF THE F PROBABILITIES WAS
910 REM PROGRAMMED BY VICTOR E. MCGEE, PSYCHOLOGY DEPARTMENT, 646-2771
920 LET P=1
930 IF G<1 THEN 980
940 LET A=M
950 LET B=N
960 LET F=G
970 GO TO 1010
980 LET A=N
990 LET B=M
1000 LET F=1/G
1010 LET A1=2/(9*A)
1020 LET B1=2/(9*B)
1030 LET Z=ABS((1-B1)*F^(.333333)-1+A1)
1040 LET Z=Z/SQR(B1*F^(.666667)+A1)
1050 IF B<4 THEN 1090
1060 LET P=(1+Z*(.196854+Z*(.115194+Z*(.000344+Z*.019527))))^4
1070 LET P=.5/P
1080 GO TO 1110
1090 LET Z=Z*(1+.08*Z^4/B^3)
1100 GO TO 1060
1110 IF G<1 THEN 1130
1120 GO TO 1150
1130 LET P=1-P
1140 GO TO 1150
1150 PRINT
1160 LET P = INT(100000*P)/100000
1170 PRINT "EXACT PROB. OF F=";G;"WITH ( "M;", "N;" ) D.F. IS ";P
1180 PRINT
1190 RETURN
1200 DATA 4,3
1210 DATA 1,2,3
1220 DATA 4,1,2
1230 DATA 2,3,4
1240 DATA 3,4,1
1250 DATA 1,1,1,0
1260 DATA 1,1,0,1
1270 DATA 0,1,1,1
1280 DATA 1,0,1,1
1290 DATA 2,1,0
1300 DATA -2,2,2
1310 DATA -1,-1,-3
1320 DATA 0,-4,2
1330END