1 3REM P O L F I T
  2 4REM
  3 7REM       (DO NOT RESEQUENCE)
  4 8 PRINT"POLFIT"
  5 10 GO TO 1890
  6 19 READ M,N
  7 20 DIM A(15),B(15),S(15),G(16),U(15)
  8 25 DIM Q(100),P(100),X(100),Y(100),C(100)
  9 30 LET Z=0
 10 33 LET O = 1
 11 34 DEF FNL(P) = .43429448 * LOG(P)
 12 35 LET O9 = 1E-19
 13 36 LET O8 = 1E19
 14 37 LET O7 = 1E38
 15 38 LET O6 = 1E-38
 16 39 DEF FNR(P) = INT(P*1000000+.5)/1000000
 17 40 LET K=12
 18 45 LET N=N+1
 19 50 IF N > 12 THEN 1870
 20 55 IF M <  N THEN 2170
 21 60 IF M >100 THEN 1840
 22 70 LET T7=Z
 23 75 LET T8=Z
 24 80 LET W7=Z
 25 100 DATA 1,38.1,1.5,14.67,2,12.76,2.5,13.15,3,11.78
 26 101 DATA 3.5,1.67,4,5.35,4.5,14.6,5,5.3,5.5,1.67,6
 27 102 DATA 8.91,6.5,15.67
 28 300 FOR I=1 TO M
 29 310 READ X(I),Y(I)
 30 320 LET W7=W7+X(I)
 31 330 LET T7=T7+Y(I)
 32 340 LET T8=T8+Y(I)^2
 33 350 NEXT I
 34 360 LET T9=(M*T8-T7^2)/(M^2-M)
 35 370 PRINT
 36 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"
 37 390 PRINT
 38 391 PRINT "VERSION OF 9/16/69"
 39 392 PRINT
 40 400 PRINT "       NUMBER OF POINTS =";M
 41 410 PRINT "        MEAN VALUE OF X =";W7/M
 42 420 PRINT "        MEAN VALUE OF Y =";T7/M
 43 430 PRINT "         STD ERROR OF Y =";SQR(T9)
 44 440 PRINT
 45 450 PRINT "   NOTE:  CODE FOR 'WHAT NEXT?' IS:"
 46 460 PRINT
 47 470 PRINT "           0 = STOP PROGRAM"
 48 480 PRINT "           1 = COEFFICIENTS ONLY"
 49 490 PRINT "           2 = ENTIRE SUMMARY"
 50 500 PRINT "           3 = FIT NEXT HIGHER DEGREE"
 51 510 PRINT
 52 520 PRINT
 53 530 FOR I=1 TO M
 54 540 LET P(I) = Z
 55 550 LET Q(I) = O
 56 560 NEXT I
 57 570 FOR I = 1 TO 11
 58 580 LET A(I) = Z
 59 590 LET B(I) = Z
 60 600 LET S(I) = Z
 61 610 NEXT I
 62 620 LET E1=Z
 63 630 LET F1=Z
 64 640 LET W1=M
 65 650 LET N4=K
 66 660 LET I=1
 67 670 LET K1=2
 68 680 IF N=0 THEN 700
 69 690 LET K1=N4
 70 700 LET W=Z
 71 710 FOR L=1 TO M
 72 720 LET W=W+Y(L)*Q(L)
 73 730 NEXT L
 74 740 LET S(I)=W/W1
 75 750 IF I-N4>=0 THEN 1090
 76 760 IF I-M>=0 THEN 1090
 77 770 LET E1=Z
 78 780 FOR L=1 TO M
 79 790 LET A9 = ABS(Q(L))
 80 800 IF A9 < O9 THEN 860
 81 801 LET X9 = ABS(X(L))
 82 802 IF X9 < O6 THEN 844
 83 803 REM: SO  X & Q  NOT TOO CLOSE TO 0.0 FOR LOG
 84 810 LET L2 = FNL(X9) + 2*FNL(A9)
 85 820 IF L2 < 38 THEN 850
 86 821 REM: SO  X & Q  TOO BIG FOR  X*Q^2
 87 830 LET E1 = O7
 88 831 LET L2 = 38
 89 840 GO TO 870
 90 844 IF A9 < 1  THEN 860
 91 846 IF A9 < O8 THEN 850
 92 847 REM: SO  Q  TOO BIG FOR  Q^2
 93 849 GO TO 830
 94 850 LET E1 = E1 + X(L) * A9^2
 95 860 NEXT L
 96 870 IF L2 - FNL(W1) > (-38) THEN 900
 97 880 LET E1 = 0
 98 890 GO TO 910
 99 900 LET E1 = E1/W1
100 910 LET A(I+1)=E1
101 920 LET W=Z
102 930 FOR L=1 TO M
103 940 LET V=(X(L)-E1)*Q(L)-F1*P(L)
104 950 LET P(L)=Q(L)
105 960 LET Q(L)=V
106 970 LET V9 = ABS(V)
107 980 IF V9 < O9 THEN 1030
108 990 IF V9 < O8 THEN 1020
109 1000 LET W = O7
110 1010 GO TO 1040
111 1020 LET W=W+V*V
112 1030 NEXT L
113 1040 LET F1= W/W1
114 1050 LET B(I+2)=F1
115 1060 LET W1=W
116 1070 LET I=I+1
117 1080 GOTO 700
118 1090 FOR L = 1 TO 13
119 1100 LET G(L)=Z
120 1110 NEXT L
121 1120 REM:
122 1130 LET G(2) = O
123 1140 FOR J=1 TO N
124 1150 LET S1 = Z
125 1160 FOR L = 2 TO N+1
126 1170 IF L=2 THEN 1190
127 1180 LET G(L) = G(L) - A(L-1)*G(L-1) - B(L-1)*G(L-2)
128 1190 LET S1 = S1 + S(L-1)*G(L)
129 1200 NEXT L
130 1210 LET U(J)=S1
131 1220 REM:
132 1230 LET L = N+1
133 1240 FOR I2=2 TO N
134 1250 LET G(L)=G(L-1)
135 1260 LET L=L-1
136 1270 NEXT I2
137 1280 LET G(2) = Z
138 1290 NEXT J
139 1300 REM:
140 1310 PRINT
141 1320 LET T=Z
142 1330 FOR L=1 TO M
143 1340 LET C(L)=Z
144 1350 LET J=N
145 1360 FOR I2=1 TO N
146 1370 LET C(L)=C(L)*X(L)+U(J)
147 1380 LET J=J-1
148 1390 NEXT I2
149 1400 LET T3=Y(L)-C(L)
150 1410 LET T=T+T3^2
151 1420 NEXT L
152 1430 IF M<>N THEN 1460
153 1440 LET T5=0
154 1450 GOTO 1470
155 1460 LET T5=T/(M-N)
156 1470 LET Q7 = 1 - T/(T9*(M-1))
157 1490 PRINT "POLYFIT OF DEGREE";N-1;
158 1500 PRINT "INDEX OF DETERM =";FNR(Q7);
159 1510 GOSUB 2200
160 1540 IF R=0 THEN 9999
161 1550 IF R=3 THEN 1810
162 1555 PRINT""
163 1560 PRINT TAB(8),  "TERM "," COEFFICIENT"
164 1570 PRINT
165 1580 FOR J=1 TO N
166 1590 LET I2=J-1
167 1600 PRINT I2,U(J)
168 1610 NEXT J
169 1620 IF R=1 THEN 1780
170 1630 PRINT
171 1640 PRINT "    X-ACTUAL","    Y-ACTUAL","      Y-CALC","        DIFF";
172 1650 PRINT "       PCT-DIFF"
173 1660 PRINT
174 1670 FOR L=1 TO M
175 1680 LET Q8=Y(L)-C(L)
176 1690 PRINT X(L),Y(L),C(L),Q8,
177 1700 IF C(L)=0 THEN 1730
178 1710 PRINT 100*Q8/C(L)
179 1720 GOTO 1740
180 1730 PRINT "INFINITE"
181 1740 NEXT L
182 1750 PRINT
183 1760 PRINT "          STD ERROR OF ESTIMATE FOR Y =";SQR(T5)
184 1770 IF K=N THEN 9999
185 1780 PRINT
186 1790 GOSUB 2200
187 1800 GOTO 1540
188 1810 LET N=N+1
189 1820 IF M<N THEN 2170
190 1830 GOTO 1090
191 1840 PRINT
192 1850 PRINT "PROGRAM SIZE LIMIT IS 100 DATA POINTS."
193 1860 GOTO 9999
194 1870 PRINT "ELEVENTH DEGREE IS THE LIMIT."
195 1880 GOTO 9999
196 1890 PRINT "MINIMUM INSTRUCTION = 0,  ALL = OTHER #, WHICH...";
197 1900 INPUT O0
198 1910 IF O0=0 THEN 2090
199 1920 PRINT
200 1930 PRINT "REVISED  9/16/69"
201 1940 PRINT
202 1950 PRINT "FITS LEAST-SQUARES POLYNOMIALS TO BIVARIATE DATA BY";
203 1960 PRINT " ORTHOGONAL POLYNMLS."
204 1970 PRINT
205 1980 PRINT "   LIMITS :   11-TH DEGREE FIT   &   MAX OF 100 DATA POINTS."
206 1985 PRINT"  BUT FITS HIGHER THAN DEGREE 5 MAY GIVE POOR RESULTS."
207 1990 PRINT "   POLFIT ALLOWS USER TO SPECIFY THE LOWEST DEGREE POLYNOM";
208 2000 PRINT "IAL TO BE FIT"
209 2010 PRINT"   AND THEN FITS POLYNOMIALS IN ORDER OF ASCENDING DEGREE."
210 2020 PRINT"   AT EACH STAGE, THE INDEX OF DETERMINATION IS PRINTED, AND"
211 2030 PRINT"   THEN THE USER HAS THE CHOICE OF:"
212 2040 PRINT"       GOING TO NEXT HIGHER FIT, OR"
213 2050 PRINT"       GETTING ONE OF TWO SUMMARIES, OR"
214 2060 PRINT"       STOPPING THE PROGRAM."
215 2070 PRINT
216 2080 PRINT "TO USE, TYPE:"
217 2090 PRINT
218 2100 PRINT "   10   DATA N,D"
219 2110 PRINT "          (WHERE N = NUMBER OF DATA POINTS TO BE READ"
220 2120 PRINT "             AND D = INITIAL [LOWEST] DEGREE TO BE FIT)"
221 2130 PRINT "   100 DATA X(1),Y(1),X(2),Y(2),....,X(N),Y(N)"
222 2140 PRINT "         (CONTINUATION ON LINES 101-299 AS NEEDED)"
223 2150 PRINT "   RUN"
224 2160 GO TO 9999
225 2170 PRINT
226 2180 PRINT "TOO FEW POINTS FOR FITTING DEGREE";N-1
227 2190 GOTO 9999
228 2200 PRINT "   WHAT NEXT";
229 2210 INPUT R
230 2220 RETURN
231 9999 END