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