Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-01 - decus/20-0025/polfit.bas
There are 2 other files named polfit.bas in the archive. Click here to see a list.
```10	GO TO 1750
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
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

```