1/* lsquares.mac -- fit parameters to data by the method of least squares
2 *
3 * Copyright 2007 by Robert Dodier.
4 * I release this file under the terms of the GNU General Public License.
5 *
6 * Examples:
7
8    load (lsquares);
9
10    /* (1) Linear regression. Exact solution is computed by solve. */
11
12    M1 : matrix ([1, -3, 2, 1], [-2, 1, 3, -1], [4, 0, -2, 1],
13                 [1, 2, 0, -1], [0, 2, 1, -2]);
14
15    lsquares_mse (M1, [w, x, y, z], w = a*x + b*y + c*z + d);
16        => ('sum ((-c*'M1[i,4] - b*'M1[i,3] - a*'M1[i,2] + 'M1[i,1] - d)^2, i, 1, 5))/5
17
18    ''%, nouns;
19        => ((-d + 2*c - b - 2*a)^2 + (-d + c - 3*b - a - 2)^2
20          + (-d + c - 2*a + 1)^2 + (- d - c + 2*b + 4)^2
21          + (- d - c - 2*b + 3*a + 1)^2)/5
22
23    a1 : lsquares_estimates (M1, [w, x, y, z], w = a*x + b*y + c*z + d, [a, b, c, d]);
24        => [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]]
25
26    lsquares_residuals (M1, [w, x, y, z], w = a*x + b*y + c*z + d, first (a1));
27        => [3/425, -4/425, -11/850, 23/850, -1/85]
28
29    lsquares_residual_mse (M1, [w, x, y, z], w = a*x + b*y + c*z + d, first (a1));
30        => 1/4250
31
32    /* (2a) Nonlinear regression. Exact solution (perfect fit) is computed by solve. */
33
34    M2a : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2]);
35
36    lsquares_mse (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C);
37        => ('sum (((D + 'M2a[i, 1])^2 - C - 'M2a[i, 3] * B - 'M2a[i, 2]*A)^2, i, 1, 4))/4
38
39    ''%, nouns;
40        => (((D + 3)^2 - C - 2*B - 2*A)^2 + ((D + 9/4)^2 - C - B - 2*A)^2
41          + ((D + 3/2)^2 - C - 2*B - A)^2 + ((D + 1)^2 - C - B - A)^2)/4
42
43    a2a : lsquares_estimates (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
44        => [[A = -75/8, B = -33/8, C = 2089/64, D = -43/8]]
45
46    lsquares_residuals (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2a));
47        => [0, 0, 0, 0]
48
49    lsquares_residual_mse (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2a));
50        => 0
51
52    /* (2b) Nonlinear regression. Exact solution (imperfect fit) is computed by solve. */
53
54    M2b : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
55
56    lsquares_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C);
57        => ('sum (((D + 'M2b[i, 1])^2 - C - 'M2b[i, 3]*B - 'M2b[i, 2]*A)^2, i, 1, 5))/5
58
59    ''%, nouns;
60        => (((D + 3)^2 - C - 2*B - 2*A)^2 + ((D + 9/4)^2 - C - B - 2*A)^2
61          + ((D + 2)^2 - C - B - 2*A)^2 + ((D + 3/2)^2 - C - 2*B - A)^2
62          + ((D + 1)^2 - C - B - A)^2)/5
63
64    a2b : lsquares_estimates (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
65        => [[A = -59/16, B = -27/16, C = 10921/1024, D = -107/32]]
66
67    float (%);
68        => [[A = - 3.6875, B = - 1.6875, C = 10.6650390625, D = - 3.34375]]
69
70    lsquares_residual_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2b));
71        => 169/2560
72
73    float (%);
74        => 0.066015625
75
76    /* (2c) Nonlinear regression. Approximate solution is computed by lbfgs.
77     *      Same data as for problem 2b.
78     */
79
80    MSE : lsquares_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C);
81
82    a2c : lsquares_estimates_approximate (MSE, [A, B, C, D]);
83        => [[A = - 3.678504947401859, B = - 1.683070351177876,
84            C = 10.63469950148675, D = - 3.340357993175253]]
85
86    lsquares_residual_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2c));
87        => 0.066016442308688
88
89    /* (3) Nonlinear regression. Approximate solution is computed by lbfgs. */
90
91    M3 : matrix ([1, 1], [2, 7/4], [3, 11/4], [4, 13/4]);
92
93    lsquares_mse (M3, [x, y], y = a*x^b + c);
94        => ('sum (('M3[i, 2] - a*'M3[i, 1]^b - c)^2, i, 1, 4))/4
95
96    ''%, nouns;
97        => ((- c - a*4^b + 13/4)^2 + (- c - a*3^b + 11/4)^2
98          + (- c - a*2^b + 7/4)^2 +(- c - a + 1)^2)/4
99
100    a3 : lsquares_estimates (M3, [x, y], y = a*x^b + c, [a, b, c], initial = [3, 3, 3]);
101        => [[a = 1.387365874920709, b = .7110956639593544, c = - .4142705622439636]]
102
103    lsquares_residuals (M3, [x, y], y = a*x^b + c, first (a3));
104        => [.02690468732325513, -.1069124575272924,
105            0.134080543273408, -.05376258426981284]
106
107    lsquares_residual_mse (M3, [x, y], y = a*x^b + c, first (a3));
108        => .008255535831587049
109
110    /* (4) Presence of unused matrix columns has no effect on result. */
111
112    M1_padded :
113        transpose (apply (matrix, apply (append,
114            map (lambda ([r], [r, 2*r]), args (transpose (M1))))));
115
116    lsquares_estimates (M1_padded, [w, w2, x, x2, y, y2, z, z2],
117        w = a*x + b*y + c*z + d, [a, b, c, d]);
118        => [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]]
119
120    /* (5) Permutation of matrix columns has no effect on result. */
121
122    p : [4, 2, 1, 3];
123    M1_permutation :
124        transpose (apply (matrix, block ([L : args (transpose (M1))],
125            makelist (L[i], i, p))));
126
127    lsquares_estimates (M1_permutation, [z, x, w, y],
128        w = a*x + b*y + c*z + d, [a, b, c, d]);
129        => [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]]
130
131    /* (6a) From the reference manual. Exact solution. */
132
133    A : matrix ([1, 2, 0], [3, 5, 4], [4, 7, 9], [5, 8, 10]);
134    soln : lsquares_estimates (A, [x, y, z], z = a*x*y + b*x + c*y + d, [a, b, c, d]);
135
136        => [[a = 1/6, b = -29/6, c = 23/6, d = -19/6]]
137
138    ratsimp (ev (z = a*x*y + b*x + c*y + d, soln [1]));
139
140        => z = ((x + 23)*y - 29*x - 19)/6
141
142    lsquares_residual_mse (A, [x, y, z], z = a*x*y + b*x + c*y + d, soln [1]);
143
144        => 0
145
146    /* (6b) From the reference manual. Exact solution. */
147
148    kill (values);
149    A : matrix ([0, 0], [1, 0], [2, 0], [3, 8], [4, 44]);
150    soln : lsquares_estimates (A, [n, p], p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, [a0, a1, a2, a3, a4]);
151
152        => [[a0 = 0, a1 = -1/3, a2 = 3/2, a3 = -5/3, a4 = 1/2]]
153
154    ratsimp (ev (p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, soln [1]));
155
156        => p = (3*n^4 - 10*n^3 + 9*n^2 - 2*n)/6
157
158    lsquares_residual_mse (A, [n, p], p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, soln [1]);
159
160        => 0
161
162    /* (6c) From the reference manual. Exact solution. */
163
164    A : matrix ([1, 7], [2, 13], [3, 25]);
165    soln : lsquares_estimates (A, [x, y], (y + c)^2 = a*x + b, [a, b, c]);
166
167        => [[a = - 216, b = 657, c = - 28]]
168
169    ev ((y + c)^2 = a*x + b, soln [1]);
170
171        => (y - 28)^2 = 657 - 216*x
172
173    lsquares_residual_mse (A, [x, y], (y + c)^2 = a*x + b, soln [1]);
174
175        => 0
176
177    /* (6d) From the reference manual. Exact solution. */
178
179    A : matrix ([1, 7], [2, 13], [3, 25], [4, 49]);
180    soln : lsquares_estimates (A, [x, y], y = a*b^x + c, [a, b, c]);
181
182        => [[a = 3, b = 2, c = 1]]
183
184    ev (y = a*b^x + c, soln [1]);
185
186        => y = 3*2^x + 1
187
188    lsquares_residual_mse (A, [x, y], y = a*b^x + c, soln [1]);
189
190        => 0
191
192    /* (7a) From the reference manual. Approximate solution. */
193
194    B : matrix ([1.1, 7.1], [2.1, 13.1], [3.1, 25.1], [4.1, 49.1]);
195    soln : lsquares_estimates (B, [x, y], y = a*b^x + c, [a, b, c], initial = [5, 5, 5], tol = 1e-8);
196
197        => [[a = 2.799098974486688, b = 2.000000000018375, c = 1.100000000358122]]
198
199    lsquares_residual_mse (B, [x, y], y = a*b^x + c, soln [1]);
200
201        => 7.353419577383337E-21
202
203    /* (7b) From the reference manual. Approximate solution. */
204
205    B : matrix ([1.1, 4.1], [4.1, 7.1], [9.1, 10.1], [16.1, 13.1]);
206    soln : lsquares_estimates (B, [x, y], y = a*x^b + c, [a, b, c], initial = [4, 1, 2], tol = 1e-8);
207
208        => [[a = 3.17731589660838, b = .4878659751689245, c = .7723843418856658]]
209
210    lsquares_residual_mse (B, [x, y], y = a*x^b + c, soln [1]);
211
212        => 6.822196230278462E-6
213
214    /* (7c) From the reference manual. Approximate solution. */
215
216    kill (A, B);
217    data : matrix ([0, 2, 4], [3, 3, 5], [8, 6, 6]);
218    soln : lsquares_estimates (data, [m, n, y], y = (A*m + B*n)^(1/3) + C, [A, B, C], initial = [3, 3, 3], tol = 1e-8);
219
220        => [[A = 4.999999920039424, B = 3.999999308815936, C = 2.000000105365426]]
221
222    lsquares_residual_mse (data, [m, n, y], y = (A*m + B*n)^(1/3) + C, soln [1]);
223
224        => 1.71963942036564E-16
225
226 */
227
228/* BEGIN STUFF COPIED FROM AUGMENTED_LAGRANGIAN.MAC */
229/* fboundp detects various kinds of Maxima functions */
230fboundp (a) :=
231    ?fboundp (a) # false
232 or ?mget (a, ?mexpr) # false
233 or ?get (a, ?mfexpr\*) # false
234 or ?mget (a, ?mmacro) # false;
235
236with_parameters ([L]) ::= buildq ([a : subst (":", "=",  ev (L [1])), e : rest (L)], block (a, splice (e)));
237/* END STUFF COPIED FROM AUGMENTED_LAGRANGIAN.MAC */
238
239if not fboundp (lbfgs) then load (lbfgs);
240
241lsquares_estimates (data, variables, equation, parameters, [optional_args]) := block
242
243  ([MSE : lsquares_mse (data, variables, equation)],
244
245    lsquares_estimates_exact (MSE, parameters),
246
247    if %% # []
248        then %%
249        else apply (lsquares_estimates_approximate,
250                    append ([MSE, parameters], optional_args)));
251
252lsquares_mse ('data%, variables, equation) := block
253
254   (if not atom (data%)
255        then block ([temp : ?gentemp (sconcat ("$M"))], temp :: ev(data%), data% : temp),
256
257    makelist ('data% ['i, j], j, 1, length (variables)),
258    map ("=", variables, %%),
259    psubst (%%, (lhs (equation) - rhs (equation))^2),
260    buildq
261        ([n : length (ev (data%)), summand : %%],
262        (1/n) * 'sum (summand, i, 1, n)),
263    subst (nounify (data%), nounify ('data%), %%));
264
265/* Some evaluation gyrations can be avoided by working with an array ...
266 *
267lsquares_mse_with_array (data%,variables,equation):=block
268   (makelist(data%['i,j],j,0,length(variables)-1),
269    map("=",variables,%%),
270    psubst(%%,(lhs(equation)-rhs(equation))^2),
271    buildq([n:array_nrows(data%),summand:%%],
272    ('sum(summand,i,0,n-1))/n));
273 *
274 */
275
276array_nrows (a) := 1 + first (last (apply (arrayinfo, [a])));
277
278lsquares_estimates_exact (MSE, parameters) := block
279
280  ([keepfloat : true,
281    ratprint : false,
282    realonly : true,
283    solveradcan : true,
284    %e_to_numlog : true,
285    logexpand : all,
286    solveexplicit : true,
287    equations],
288
289    MSE : ev (MSE, nouns),
290    equations : map (lambda ([x], diff (%%, x)), parameters),
291    solutions : errcatch (solve (equations, parameters)),
292
293    if solutions = []
294    then []
295    else
296        /* solve finds stationary points of the MSE.
297         * Of the points returned, find the one or ones with least MSE.
298         */
299       (solutions : solutions [1], /* errcatch puts on extra layer of []'s ... */
300        map (lambda ([x], ev (MSE, x)), solutions),
301        sublist_indices (%%, lambda ([x], x = lmin (%%))),
302        makelist (solutions [i], i, %%)));
303
304lsquares_estimates_approximate (MSE, parameters, [optional_args]) := block
305
306  ([initial : makelist (1, i, 1, length (parameters)),
307    iprint : [1, 0],
308    tol : 1e-3],
309
310    with_parameters (optional_args,
311        lbfgs (MSE, parameters, initial, tol, iprint),
312        [%%]));
313
314lsquares_residuals (data, variables, equation, parameters_estimates) :=
315
316   (equation : psubst (parameters_estimates, equation),
317    maplist (lambda ([row], (psubst (map ("=", variables, row), equation), lhs (%%) - rhs (%%))), data));
318
319lsquares_residual_mse (data, variables, equation, parameters_estimates) :=
320
321   (lsquares_residuals (data, variables, equation, parameters_estimates),
322    (1 / length (%%)) * apply ("+", map (lambda ([e], e^2), %%)));
323
324