 Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-01 - decus/20-0020/stat20.sta
There are 2 other files named stat20.sta in the archive. Click here to see a list.
```100'  NAME--STAT20
110'
120'  DESCRIPTION--MULTIPLE LINEAR REGRESSION ACCORDING TO
130'  EFROYMSON'S ALGORITHM.
140'
150'  SOURCE--RALSTON AND WILF,"MATHEMATICAL METHODS FOR DIGITAL
160'  COMPUTERS",P.191.
170'
180'  INSTRUCTIONS--ENTER DATA STARTING IN LINE 1610.
190'  FIRST DATA IS N, THE NUMBER OF INDEPENDENT VARIABLES, THEN
200'  P, THE NUMBER OF DEPENDENT VARIABLES, THEN M, THE NUMBER
210'  OF DATA SETS, THEN F1,THE VALUE OF F FOR ENTERING A
220'  A VARIABLE, AND THEN F2 THE VALUE OF F FOR REMOVING A VARIABLE.
230'  THEN ENTER THE DATA BY SET,INDEPENDENT VARIABLES BEFORE
240'  DEPENDENT VARIABLES. IF M>50 OR P>7 OR N>8 THE DIM STATEMENTS IN
250'  LINE 300 SHOULD BE CHANGED. SAMPLE DATA IS IN LINES 1610-1760.
260'  BE SURE TO REMOVE SAMPLE DATA BEFORE RUNNING PROGRAM.
270'
280'  *  *  *  *  *  *  *  MAIN PROGRAM  *  *  *  *  *  *  *  *  *  *
290'
300  DIM S(15,15),A(50,8),M(8),D(8),B(7),E(7)
310    READ N, P, M, F1, F2
320  LET T1 = 1E-6
330  LET N1 = N + P
340  LET N2 = N1 + N
350  LET K = 1
360  LET D1 = M - 1
380  FOR I = 1 TO M
390     LET A(I,0)=1
400  NEXT I
410    PRINT "MEANS "
420  FOR I = 0 TO N1
430     FOR J = I TO N1
440        LET S = 0
450        FOR L = 1 TO M
460           LET S = S + A(L,I) * A(L,J)
470        NEXT L
480        LET S(I,J) = S
490     NEXT J
500     LET M(I) = S(0,I) / S(0,0)
510     IF I = 0 THEN 530
520       PRINT M(I),
530  NEXT I
540    PRINT
550    PRINT
560    PRINT "STANDARD DEVIATIONS"
570  LET M1 = M * D1
580  FOR I = 1 TO N1
590     FOR J = I TO N1
600        LET S(I,J) = ( M*S(I,J) - S(0,I)*S(0,J) )/M1
610     NEXT J
620     LET D(I) = SQR( S(I,I) )
630       PRINT D(I),
640  NEXT I
650    PRINT
660    PRINT
670    PRINT "CORRELATION COEFFICIENTS"
680  FOR I = 1 TO N2
690     FOR J = I TO N2
700        IF J > N1 THEN 740
710        LET S(I,J) = S(I,J) / D(I) / D(J)
720        LET S(J,I) = S(I,J)
730        GO TO 800
740        IF I <> J - N1 THEN 780
750        LET S(I,J) = 1
760        LET S(J,I) = -1
770        GO TO 800
780        LET S(I,J) = 0
790        LET S(J,I) = 0
800     NEXT J
810  NEXT I
820  FOR I = 1 TO N1
830     FOR J = 1 TO N1
840          PRINT S(I,J),
850     NEXT J
860       PRINT
870       PRINT
880  NEXT I
890    PRINT
900  LET K1 = K + N
910    PRINT "DEPENDENT VARIABLE",K
920  FOR I = 1 TO N
930     LET B(I) = 0
940  NEXT I
950  LET S8 = D(K1) * SQR( S(K1,K1) / D1 )
960 LET V1 = 1E35
970  LET V2 = 0
980  LET N3 = 0
990  LET N4 = 0
1000  FOR I = 1 TO N
1010     IF ABS( S(I,I) ) <= T1 THEN 1140
1020     LET V0 = S(I,K1) * S(K1,I) / S(I,I)
1030     ON SGN(V0)+2 GOTO 1080,1140,1040
1040     IF V0 <= V2 THEN 1140
1050     LET V2 = V0
1060     LET N4 = I
1070     GO TO 1140
1080     LET I1 = I + N1
1090     LET B(I) = S(I1,K1) * D(K1) / D(I)
1100     LET E(I) = S8 / D(I) * SQR( S(I1,I1) )
1110     IF ABS( V0 ) >= ABS( V1 ) THEN 1140
1120     LET V1 = V0
1130     LET N3 = I
1140  NEXT I
1150  LET S = 0
1160  FOR I = 1 TO N
1170     LET S = S + B(I) * M(I)
1180  NEXT I
1190  LET B(0) = M(K1) - S
1200    PRINT "INDEX ", "B  ", "STD. DEV.", "T-RATIO  "
1210    PRINT 0, B(0)
1220  FOR I = 1 TO N
1230     IF B(I) = 0 THEN 1250
1240       PRINT I, B(I), E(I), B(I) / E(I)
1250  NEXT I
1260    PRINT "STANDARD ERROR OF Y  ", S8 * SQR(M-1)
1270    PRINT
1280  LET F = ABS(V1) * D1 / S(K1,K1)
1290  IF F < F2 THEN 1360
1300  LET F = V2 * (D1 - 1) / ( S(K1,K1) - V2 )
1310  IF F <= F1 THEN 1570
1320  LET Q = N4
1330  LET D1 = D1 - 1
1340    PRINT "VARIABLE ENTERING ", Q
1350  GO TO 1390
1360  LET Q = N3
1370  LET D1 = D1 + 1
1380    PRINT "VARIABLE LEAVING  ", Q
1390    PRINT "F-LEVEL  ", F
1400  LET Y = 1 / S(Q,Q)
1410  FOR J = 1 TO N2
1420     LET S(Q,J) = S(Q,J) * Y
1430  NEXT J
1440  LET Y = -Y
1450  FOR I = 1 TO N2
1460     IF I = Q THEN 1540
1470     LET X = -S(I,Q)
1480     FOR J = 1 TO N2
1490        IF J = Q THEN 1520
1500        LET S(I,J) = S(I,J) + X * S(Q,J)
1510        GO TO 1530
1520        LET S(I,J) = S(I,J) * Y
1530     NEXT J
1540  NEXT I
1550  LET S(Q,Q) = -Y
1560  GOTO 950
1570  LET K = K + 1
1580  IF K <= P THEN 900
1590    PRINT "   *****  END OF PROBLEM  *****"
1600  STOP
1610  DATA  3, 1, 15, 2.500, 2.500
1620  DATA  32, 48, 54, 15
1630  DATA  36, 33, 19, 16
1640  DATA   3, 28, 30, 14
1650  DATA  12, 33, 64, 22
1660  DATA  36, 34, 60, 24
1670  DATA  24, 36, 53, 19
1680  DATA  19, 42, 29, 13
1690  DATA  20, 33, 55, 15
1700  DATA  27, 36, 62, 23
1710  DATA  15, 22, 33, 12
1720  DATA  45, 46, 68, 25
1730  DATA   9, 28, 42, 17
1740  DATA  11, 32, 45, 18
1750  DATA  33, 34, 39, 19
1760  DATA  21, 45, 39, 18
1770END

```