share/lsquares/lsquares.mac

   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 */
 230 fboundp (a) :=
 231     ?fboundp (a) # false
 232  or ?mget (a, ?mexpr) # false
 233  or ?get (a, ?mfexpr\*) # false
 234  or ?mget (a, ?mmacro) # false;
 235
 236 with_parameters ([L]) ::= buildq ([a : subst (":", "=",  ev (L [1])), e : rest (L)], block (a, splice (e)));
 237 /* END STUFF COPIED FROM AUGMENTED_LAGRANGIAN.MAC */
 238
 239 if not fboundp (lbfgs) then load (lbfgs);
 240
 241 lsquares_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
 252 lsquares_mse ('data%, variables, equation) := block
 253
 254    (if not atom (data%)
 255         then block ([temp : gensym ("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  *
 267 lsquares_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
 276 array_nrows (a) := 1 + first (last (apply (arrayinfo, [a])));
 277
 278 lsquares_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
 304 lsquares_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
 314 lsquares_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
 319 lsquares_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