share/lsquares/rtest_lsquares.mac

   1 /* TODO: Do all the tolerances I have chosen here make any sense?*/
   2 (kill (all),
   3  load (lsquares),
   4  save_tolerance : float_approx_equal_tolerance,
   5  float_approx_equal_tolerance : 1e-8,
   6  0);
   7 0$
   8
   9 /* Example (1) from lsquares.mac comment header */
  10
  11 (M1 : matrix ([1, -3, 2, 1], [-2, 1, 3, -1], [4, 0, -2, 1], [1, 2, 0, -1], [0, 2, 1, -2]),
  12  mse : lsquares_mse (M1, [w, x, y, z], w = a*x + b*y + c*z + d));
  13 ('sum((-c*'M1[i,4]-b*'M1[i,3]-a*'M1[i,2]+'M1[i,1]-d)^2,i,1,5))/5;
  14
  15 ev (mse, nouns);
  16 ''(((-d+2*c-b-2*a)^2+(-d+c-3*b-a-2)^2+(-d+c-2*a+1)^2+(-d-c+2*b+4)^2 +(-d-c-2*b+3*a+1)^2) /5);
  17
  18 a1 : lsquares_estimates (M1, [w, x, y, z], w = a*x + b*y + c*z + d, [a, b, c, d]);
  19 [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]];
  20
  21 lsquares_residuals (M1, [w, x, y, z], w = a*x + b*y + c*z + d, first (a1));
  22 [3/425, -4/425, -11/850, 23/850, -1/85]$
  23
  24 lsquares_residual_mse (M1, [w, x, y, z], w = a*x + b*y + c*z + d, first (a1));
  25 1/4250;
  26
  27 /* Example (2a) from lsquares.mac comment header */
  28
  29 (M2a : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2]),
  30  mse : lsquares_mse (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C));
  31 ('sum(((D+'M2a[i,1])^2-C-'M2a[i,3]*B-'M2a[i,2]*A)^2,i,1,4))/4;
  32
  33 ev (mse, nouns);
  34 ''((((D+3)^2-C-2*B-2*A)^2+((D+9/4)^2-C-B-2*A)^2 +((D+3/2)^2-C-2*B-A)^2+((D+1)^2-C-B-A)^2) /4);
  35
  36 a2a : lsquares_estimates (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
  37 [[A = -75/8, B = -33/8, C = 2089/64, D = -43/8]];
  38
  39 lsquares_residuals (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2a));
  40 [0, 0, 0, 0];
  41
  42 lsquares_residual_mse (M2a, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2a));
  43 0;
  44
  45 /* Example (2b) from lsquares.mac comment header */
  46 (M2b : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]),
  47  mse : lsquares_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C));
  48 ('sum(((D+'M2b[i,1])^2-C-'M2b[i,3]*B-'M2b[i,2]*A)^2,i,1,5))/5;
  49
  50 ev (mse, nouns);
  51 ''((((D+3)^2-C-2*B-2*A)^2+((D+9/4)^2-C-B-2*A)^2+((D+2)^2-C-B-2*A)^2 +((D+3/2)^2-C-2*B-A)^2+((D+1)^2-C-B-A)^2) /5);
  52
  53 a2b : lsquares_estimates (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
  54 [[A = -59/16, B = -27/16, C = 10921/1024, D = -107/32]];
  55
  56 lsquares_residual_mse (M2b, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a2b));
  57 169/2560;
  58
  59 /* Example (2c) from lsquares.mac comment header */
  60 (a2c : lsquares_estimates_approximate (mse, [A, B, C, D], iprint = [-1, 0],tol=1e-11),0);
  61 0$
  62
  63 float_approx_equal(subst(a2c,A),-3.6875);
  64 true$
  65
  66 float_approx_equal(subst(a2c,B),-1.6875);
  67 true$
  68
  69 float_approx_equal(subst(a2c,C),10.66503906);
  70 true$
  71
  72 float_approx_equal(subst(a2c,D),-3.34375);
  73 true$
  74
  75 /* Example (3) from lsquares.mac comment header */
  76 (M3 : matrix ([1, 1], [2, 7/4], [3, 11/4], [4, 13/4]),
  77  mse : lsquares_mse (M3, [x, y], y = a*x^b + c));
  78 ('sum(('M3[i,2]-a*'M3[i,1]^b-c)^2,i,1,4))/4;
  79
  80 ev (mse, nouns);
  81 ''(((-c-a*4^b+13/4)^2+(-c-a*3^b+11/4)^2+(-c-a*2^b+7/4)^2 +(-c-a+1)^2) /4);
  82
  83 (a3 : lsquares_estimates (M3, [x, y], y = a*x^b + c, [a, b, c], initial = [3, 3, 3], iprint = [-1, 0],tol=1e-11),0);
  84 0;
  85
  86 float_approx_equal(subst(a3[1],a),1.38889001796731);
  87 true$
  88
  89 float_approx_equal(subst(a3[1],b),0.7106139927150336);
  90 true$
  91
  92 float_approx_equal(subst(a3[1],c),-0.4158391087045305);
  93 true$
  94
  95 (a3:lsquares_residuals (M3, [x, y], y = a*x^b + c, first (a3)),0);
  96 0$
  97
  98 float_approx_equal(a3[1],0.02694909073722052);
  99 true$
 100
 101 float_approx_equal(a3[2],-0.1070800206548264);
 102 true$
 103
 104 float_approx_equal(a3[3],0.1339249827507731);
 105 true$
 106
 107 float_approx_equal(a3[4],-0.05379405286758576);
 108 true$
 109
 110 /* Example (4) from lsquares.mac comment header */
 111 (M1_padded : transpose (apply (matrix, apply (append, map (lambda ([r], [r, 2*r]), args (transpose (M1)))))),
 112  lsquares_estimates (M1_padded, [w, w2, x, x2, y, y2, z, z2], w = a*x + b*y + c*z + d, [a, b, c, d],tol=1e-11));
 113 [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]];
 114
 115 /* Example (5) from lsquares.mac comment header */
 116 (p : [4, 2, 1, 3],
 117  M1_permutation : transpose (apply (matrix, block ([L : args (transpose (M1))], makelist (L[i], i, p)))),
 118  lsquares_estimates (M1_permutation, [z, x, w, y], w = a*x + b*y + c*z + d, [a, b, c, d]));
 119 [[a = -19/34, b = -499/425, c = -181/850, d = 798/425]];
 120
 121 /* Example (6a) from lsquares.mac comment header */
 122 (A : matrix ([1, 2, 0], [3, 5, 4], [4, 7, 9], [5, 8, 10]),
 123  soln : lsquares_estimates (A, [x, y, z], z = a*x*y + b*x + c*y + d, [a, b, c, d]));
 124 [[a = 1/6, b = -29/6, c = 23/6, d = -19/6]];
 125
 126 ratsimp (ev (z = a*x*y + b*x + c*y + d, soln [1]));
 127 z = ((x + 23)*y - 29*x - 19)/6;
 128
 129 lsquares_residual_mse (A, [x, y, z], z = a*x*y + b*x + c*y + d, soln [1]);
 130 0;
 131
 132 /* Example (6b) from lsquares.mac comment header */
 133 (kill (n, p, a0, a1, a2, a3, a4),
 134  A : matrix ([0, 0], [1, 0], [2, 0], [3, 8], [4, 44]),
 135  soln : lsquares_estimates (A, [n, p], p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, [a0, a1, a2, a3, a4]));
 136 [[a0 = 0, a1 = -1/3, a2 = 3/2, a3 = -5/3, a4 = 1/2]];
 137
 138 ratsimp (ev (p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, soln [1]));
 139 p = (3*n^4 - 10*n^3 + 9*n^2 - 2*n)/6;
 140
 141 lsquares_residual_mse (A, [n, p], p = a4*n^4 + a3*n^3 + a2*n^2 + a1*n + a0, soln [1]);
 142 0;
 143
 144 /* Example (6c) from lsquares.mac comment header */
 145 (A : matrix ([1, 7], [2, 13], [3, 25]),
 146  soln : lsquares_estimates (A, [x, y], (y + c)^2 = a*x + b, [a, b, c]));
 147 [[a = - 216, b = 657, c = - 28]];
 148
 149 ev ((y + c)^2 = a*x + b, soln [1]);
 150 (y - 28)^2 = 657 - 216*x;
 151
 152 lsquares_residual_mse (A, [x, y], (y + c)^2 = a*x + b, soln [1]);
 153 0;
 154
 155 /* Example (6d) from lsquares.mac comment header */
 156 (A : matrix ([1, 7], [2, 13], [3, 25], [4, 49]),
 157  soln : lsquares_estimates (A, [x, y], y = a*b^x + c, [a, b, c]));
 158 [[a = 3, b = 2, c = 1]];
 159
 160 ev (y = a*b^x + c, soln [1]);
 161 y = 3*2^x + 1;
 162
 163 lsquares_residual_mse (A, [x, y], y = a*b^x + c, soln [1]);
 164 0;
 165
 166 /* Example (7a) from lsquares.mac comment header */
 167 ((B : matrix ([1.1, 7.1], [2.1, 13.1], [3.1, 25.1], [4.1, 49.1]),
 168  soln : lsquares_estimates (B, [x, y], y = a*b^x + c, [a, b, c], initial = [5, 5, 5], tol = 1e-11, iprint = [-1, 0])),0);
 169 0$
 170 float_approx_equal(subst(soln[1],a),2.799098974610445);
 171 true$
 172
 173 float_approx_equal(subst(soln[1],b),2.0);
 174 true$
 175
 176 float_approx_equal(subst(soln[1],c),1.1);
 177 true$
 178
 179 /* Example (7b) from lsquares.mac comment header */
 180 ((B : matrix ([1.1, 4.1], [4.1, 7.1], [9.1, 10.1], [16.1, 13.1]),
 181  soln : lsquares_estimates (B, [x, y], y = a*x^b + c, [a, b, c], initial = [4, 1, 2], tol = 1e-11, iprint = [-1, 0])),0);
 182 0$
 183
 184 float_approx_equal(subst(soln[1],a),3.177315891104065);
 185 true$
 186
 187 float_approx_equal(subst(soln[1],b),0.4878659755913501);
 188 true$
 189
 190 float_approx_equal(subst(soln[1],c),.7723843418856658);
 191 true$
 192
 193 /* Example (1) from lsquares.mac comment header */
 194 ((kill (A, B),
 195  data : matrix ([0, 2, 4], [3, 3, 5], [8, 6, 6]),
 196  soln : lsquares_estimates (data, [m, n, y], y = (A*m + B*n)^(1/3) + C, [A, B, C], initial = [3, 3, 3], tol = 1e-11, iprint = [-1, 0])),0)
 197 ;
 198 0;
 199
 200 float_approx_equal(subst(soln[1],A),5);
 201 true$
 202
 203 float_approx_equal(subst(soln[1],B),4);
 204 true$
 205
 206 float_approx_equal(subst(soln[1],C),1.999999999931501);
 207 true$
 208
 209 /* SF bug #3282: "lsquares and lists of list of data to be fitted on" */
 210
 211 ((kill(a, b, c, x, y),
 212  data:[matrix([0,0],[1,1],[2,2]),
 213        matrix([2,2],[1,1],[2,2])],
 214  eqtn:y=a*x^2+b*x+c,
 215  mse_data_1 : lsquares_mse(data[1],[x,y],eqtn),
 216  soln:lsquares_estimates_approximate(%%,[a,b,c])),0);
 217 0;
 218
 219 float_approx_equal(subst(soln[1],a),-3.405898581949557*10^-4);
 220 true$
 221
 222 float_approx_equal(subst(soln[1],b),1.000676542954337);
 223 true$
 224
 225 float_approx_equal(subst(soln[1],c),-1.342533507359797*10^-4);
 226 true$
 227
 228 lsquares_estimates_exact (mse_data_1, [a, b, c]);
 229 [[a = 0, b = 1, c = 0]]; /* preceding approximate result matches this */
 230
 231 ((foo(k, z) := data[k]*z,
 232  mse_foo_2_1_5 : lsquares_mse(foo(2, 1.5), [x, y], eqtn),
 233  soln:lsquares_estimates_approximate(%%, [a, b, c])),0);
 234 0$
 235
 236 float_approx_equal(subst(soln[1],a),0.132534631873083);
 237 true$
 238
 239 float_approx_equal(subst(soln[1],b),0.4036035745637151);
 240 true$
 241
 242 float_approx_equal(subst(soln[1],c),0.5963736563696966);
 243 true$
 244
 245 (reset (%rnum, %rnum_list), /* ensure result contains %r1 */
 246  lsquares_estimates_exact (mse_foo_2_1_5, [a, b, c]));
 247 [[a = (2*%r1)/9,b = 1-%r1,c = %r1]]$ /* preceding approximate result matches this */
 248
 249