Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50144/polfit.bas
There are 2 other files named polfit.bas in the archive. Click here to see a list.
10 GO TO 1750
15 READ M,N
20 DIM A(15),B(15),S(15),G(15),U(15)
25 DIM Q(100),P(100),X(100),Y(100),C(100)
30 LET Z=0
35 LET O=1
40 LET K=12
45 LET N=N+1
50 IF N>11 THEN 1730
55 IF M<N THEN 1930
60 IF M>100 THEN 1700
70 LET T7=Z
75 LET T8=Z
80 LET W7=Z
100 DATA 1,-5,2,-17,3,5,4,145,5,535,6,1355,7,2833,8,5245,9,8915,10,14215
300 FOR I=1 TO M
310 READ X(I),Y(I)
320 LET W7=W7+X(I)
330 LET T7=T7+Y(I)
340 LET T8=T8+Y(I)^2
350 NEXT I
360 LET T9=(M*T8-T7^2)/(M^2-M)
370 PRINT
380 PRINT "L E A S T - S Q U A R E S P O L Y N O M I A L S"
390 PRINT
400 PRINT " NUMBER OF POINTS =";M
410 PRINT " MEAN VALUE OF X =";W7/M
420 PRINT " MEAN VALUE OF Y =";T7/M
430 PRINT " STD ERROR OF Y =";SQR(T9)
440 PRINT
450 PRINT " NOTE: CODE FOR 'WHAT NEXT?' IS:"
460 PRINT
470 PRINT " 0 = STOP PROGRAM"
480 PRINT " 1 = COEFFICIENTS ONLY"
490 PRINT " 2 = ENTIRE SUMMARY"
500 PRINT " 3 = FIT NEXT HIGHER DEGREE"
510 PRINT
520 PRINT
530 FOR I=1 TO M
540 LET P(I) = Z
550 LET Q(I) = O
560 NEXT I
570 FOR I = 1 TO 11
580 LET A(I) = Z
590 LET B(I) = Z
600 LET S(I) = Z
610 NEXT I
620 LET E1=Z
630 LET F1=Z
640 LET W1=M
650 LET N4=K
660 LET I=1
670 LET K1=2
680 IF N=0 THEN 700
690 LET K1=N4
700 LET W=Z
710 FOR L=1 TO M
720 LET W=W+Y(L)*Q(L)
730 NEXT L
738 IF W1<>0 THEN 740
739 W1=1
740 LET S(I)=W/W1
750 IF I-N4>=0 THEN 990
760 LET E1=Z
770 FOR L=1 TO M
780 IF(ABS(Q(L)))<1E-20 THEN 800
790 LET E1=E1+X(L)*Q(L)*Q(L)
800 NEXT L
810 IF W1=0 THEN 830
820 LET E1=E1/W1
830 LET A(I+1)=E1
840 LET W=Z
850 FOR L=1 TO M
860 LET V=(X(L)-E1)*Q(L)-F1*P(L)
870 LET P(L)=Q(L)
880 LET Q(L)=V
890 IF(ABS(V))<1E-20 THEN 910
900 LET W=W+V*V
910 NEXT L
920 IF W1<>0 THEN 940
930 W1=1
940 LET F1= W/W1
950 LET B(I+2)=F1
960 LET W1=W
970 LET I=I+1
980 GOTO 700
990 FOR L=2 TO 11
1000 LET G(L)=Z
1010 NEXT L
1020 LET G(1)=O
1030 FOR J=1 TO N
1040 LET S1 =Z
1050 FOR L=1 TO N
1060 IF L=1 THEN 1080
1070 LET G(L)=G(L)-A(L)*G(L-1)-B(L)*G(L-2)
1080 LET S1=S1+S(L)*G(L)
1090 NEXT L
1100 LET U(J)=S1
1110 LET L=N
1120 FOR I2=2 TO N
1130 LET G(L)=G(L-1)
1140 LET L=L-1
1150 NEXT I2
1160 LET G(1)=Z
1170 NEXT J
1180 PRINT
1190 LET T=Z
1200 FOR L=1 TO M
1210 LET C(L)=Z
1220 LET J=N
1230 FOR I2=1 TO N
1240 LET C(L)=C(L)*X(L)+U(J)
1250 LET J=J-1
1260 NEXT I2
1270 LET T3=Y(L)-C(L)
1280 LET T=T+T3^2
1290 NEXT L
1300 IF M<>N THEN 1330
1310 LET T5=0
1320 GOTO 1340
1330 LET T5=T/(M-N)
1340 LET Q7=1-(T5/T9)
1350 PRINT
1360 PRINT "POLYFIT OF DEGREE";N-1;
1370 PRINT "INDEX OF DETERM =";Q7;
1380 GOSUB 1960
1390 PRINT
1400 PRINT
1410 IF R=0 THEN 1990
1420 IF R=3 THEN 1670
1430 PRINT "TERM","COEFFICIENT"
1440 PRINT
1450 FOR J=1 TO N
1460 LET I2=J-1
1470 PRINT I2,U(J)
1480 NEXT J
1490 IF R=1 THEN 1640
1500 PRINT
1510 PRINT "X-ACTUAL","Y-ACTUAL","Y-CALC","DIFF","PCT-DIFF"
1520 PRINT
1530 FOR L=1 TO M
1540 LET Q8=Y(L)-C(L)
1550 PRINT X(L),Y(L),C(L),Q8,
1560 IF C(L)=0 THEN 1590
1570 PRINT 100*Q8/C(L)
1580 GOTO 1600
1590 PRINT "INFINITE"
1600 NEXT L
1610 PRINT
1620 PRINT " STD ERROR OF ESTIMATE FOR Y =";SQR(T5)
1630 IF K=N THEN 1990
1640 PRINT
1650 GOSUB 1960
1660 GOTO 1410
1670 LET N=N+1
1680 IF M<N THEN 1930
1690 GOTO 990
1700 PRINT
1710 PRINT "PROGRAM SIZE LIMIT IS 100 DATA POINTS."
1720 GOTO 1990
1730 PRINT "ELEVENTH DEGREE IS THE LIMIT."
1740 GOTO 1990
1750 PRINT
1760 PRINT "THIS PROGRAM FITS LEAST-SQUARES POLYNOMIALS TO BIVARIATE"
1770 PRINT "DATA, USING AN ORTHOGONAL POLYNOMIAL METHOD. LIMITS ARE"
1780 PRINT "11-TH DEGREE FIT AND A MAX OF 100 DATA POINTS. PROGRAM"
1790 PRINT "ALLOWS USER TO SPECIFY THE LOWEST DEGREE POLYNOMIAL TO BE"
1800 PRINT "FIT, AND THEN FITS THE POLYNOMIALS IN ORDER OF ASCENDING"
1810 PRINT "DEGREE. AT EACH STAGE, THE INDEX OF DETERMINATION IS"
1820 PRINT "PRINTED, AND THE USER HAS THE CHOICE OF GOING TO THE NEXT"
1830 PRINT "HIGHER DEGREE FIT, SEEING EITHER OF TWO SUMMARIES OF FIT"
1840 PRINT "AT THAT STAGE, OR OF STOPPING THE PROGRAM. TO USE, TYPE:"
1850 PRINT
1860 PRINT " 10 DATA N, D"
1870 PRINT " (WHERE N = NUMBER OF DATA POINTS TO BE READ"
1880 PRINT " AND D = INITIAL (LOWEST) DEGREE TO BE FIT)"
1890 PRINT " 100 DATA X(1),Y(1),X(2),Y(2),....,X(N),Y(N)"
1900 PRINT " (CONTINUATION ON LINES 101-299 AS NEEDED)"
1910 PRINT " RUN"
1920 GOTO 1990
1930 PRINT
1940 PRINT "TOO FEW POINTS FOR FITTING DEGREE";N-1
1950 GO TO 1990
1960 PRINT" WHAT NEXT";
1970 INPUT R
1980 RETURN
1990 END