1/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ 2/* [ Created with wxMaxima version 20.03.1-DevelopmentSnapshot ] */ 3/* [wxMaxima: title start ] 4Fitting equations to real measurement data 5 [wxMaxima: title end ] */ 6 7 8/* [wxMaxima: comment start ] 9One place Maxima can show its full power at is where symbolical mathematics that describes a phenomenon is used for understanding real measurement data. 10 [wxMaxima: comment end ] */ 11 12 13/* [wxMaxima: comment start ] 14As a first step normally the real measurement data is loaded from a .csv file using read_matrix: 15 [wxMaxima: comment end ] */ 16 17 18/* [wxMaxima: input start ] */ 19? read_matrix; 20/* [wxMaxima: input end ] */ 21 22 23/* [wxMaxima: comment start ] 24For the sake of generating a self-contained example we generate this data synthetically. For the first example a parabola with a bit of measurement noise might be a good start: 25 [wxMaxima: comment end ] */ 26 27 28/* [wxMaxima: input start ] */ 29fun1(x):=3*x^2-2*x+7; 30/* [wxMaxima: input end ] */ 31 32 33/* [wxMaxima: input start ] */ 34time1:makelist(i,i,-10,10,.02)$ 35data1:transpose( 36 matrix( 37 time1, 38 makelist(fun1(i)+random(32.0)-16,i,time1) 39 ) 40)$ 41wxdraw2d(grid=true,xlabel="x",ylabel="y", 42 points(data1) 43)$ 44/* [wxMaxima: input end ] */ 45 46 47/* [wxMaxima: section start ] 48The general approach 49 [wxMaxima: section end ] */ 50 51 52/* [wxMaxima: comment start ] 53Let's assume we have guessed that the curve might be something parabola-like: 54 [wxMaxima: comment end ] */ 55 56 57/* [wxMaxima: input start ] */ 58approach1:y=a*x^2+b*x+c; 59/* [wxMaxima: input end ] */ 60 61 62/* [wxMaxima: comment start ] 63Maxima now offers a simple, but powerful curve fitter that guesses the parameters a,b and c for us: 64 [wxMaxima: comment end ] */ 65 66 67/* [wxMaxima: input start ] */ 68load("lsquares")$ 69/* [wxMaxima: input end ] */ 70 71 72/* [wxMaxima: input start ] */ 73lsquares_estimates_approximate( 74 lsquares_mse( 75 data1,[x,y],approach1 76 ), 77 [a,b,c], 78 initial=[0,0,0] 79); 80params1:%[1]; 81/* [wxMaxima: input end ] */ 82 83 84/* [wxMaxima: comment start ] 85In this example the "initial=" wasn't really necessary. But sometimes using the wrong starting point means that lsquares heads for the wrong local optimum of the problem. 86 [wxMaxima: comment end ] */ 87 88 89/* [wxMaxima: input start ] */ 90rec1:subst(params1,approach1); 91/* [wxMaxima: input end ] */ 92 93 94/* [wxMaxima: input start ] */ 95wxdraw2d( 96 grid=true,xlabel="x",ylabel="y", 97 color=red,key="Reconstructed", 98 explicit(rhs(rec1),x,-10,10), 99 color=blue,key="Original", 100 explicit(fun1,x,-10,10) 101)$ 102/* [wxMaxima: input end ] */ 103 104 105/* [wxMaxima: comment start ] 106Besides lsquares_estimates_approximate the lsquares package also provides a command named lsquares_estimates that tries to find the exact optimum by running the problem through solve() before finding the answer numerically. But as it is always the case with computers if the problem that is to be solved is complex one does never know in advance how long it will take and if it ever will finish. lsquares_estimates_approximate doesn't have this drawback. 107 [wxMaxima: comment end ] */ 108 109 110/* [wxMaxima: section start ] 111Dealing with non-evenly-spaced data 112 [wxMaxima: section end ] */ 113 114 115/* [wxMaxima: comment start ] 116Sometimes the input data isn't evenly spaced. 117 [wxMaxima: comment end ] */ 118 119 120/* [wxMaxima: comment start ] 121The places with a higher density of samples have more weight when fitting data to the curves. So let's add an error to the place with the highest data density and see if we can cope with it. 122 [wxMaxima: comment end ] */ 123 124 125/* [wxMaxima: input start ] */ 126time2:makelist(i^4*i/abs(i),i,-2,2,.02)$ 127data2:transpose( 128 matrix( 129 time2, 130 makelist(fun1(i)+random(4.0)-2+50*sin(i)/i,i,time2) 131 ) 132)$ 133wxdraw2d(grid=true,xlabel="x",ylabel="y", 134 points(data2), 135 yrange=[0,800] 136)$ 137/* [wxMaxima: input end ] */ 138 139 140/* [wxMaxima: comment start ] 141Since the biggest error occurs where we have the highest density of data the result of fitting the curve to this data directly will be suboptimal: 142 [wxMaxima: comment end ] */ 143 144 145/* [wxMaxima: input start ] */ 146lsquares_estimates_approximate( 147 lsquares_mse( 148 data2,[x,y],approach1 149 ), 150 [a,b,c] 151); 152params2_1:%[1]; 153wxdraw2d( 154 grid=true,xlabel="x",ylabel="y", 155 color=red,key="Reconstructed", 156 explicit(rhs(subst(params2_1,approach1)),x,-10,10), 157 color=blue,key="Original", 158 explicit(fun1,x,-10,10) 159)$ 160/* [wxMaxima: input end ] */ 161 162 163/* [wxMaxima: subsect start ] 164Converting the data to an continuous curve 165 [wxMaxima: subsect end ] */ 166 167 168/* [wxMaxima: comment start ] 169The interpol package allows to generate continuous and half-way smooth curves from any input data. 170 [wxMaxima: comment end ] */ 171 172 173/* [wxMaxima: input start ] */ 174load("interpol")$ 175/* [wxMaxima: input end ] */ 176 177 178/* [wxMaxima: input start ] */ 179data2_cont:cspline(data2,'varname=x)$ 180/* [wxMaxima: input end ] */ 181 182 183/* [wxMaxima: input start ] */ 184wxdraw2d( 185 explicit(data2_cont,x,-10,10), 186 grid=true,xlabel="x",ylabel="y" 187)$ 188/* [wxMaxima: input end ] */ 189 190 191/* [wxMaxima: subsect start ] 192Generating evenly-spaced samples from this curve 193 [wxMaxima: subsect end ] */ 194 195 196/* [wxMaxima: comment start ] 197First we generate the new data and time vector. As always the fact that a subst() is necessary in this step is caused by a bug in makelist. 198 [wxMaxima: comment end ] */ 199 200 201/* [wxMaxima: input start ] */ 202time_i:makelist(i,i,-10,10,.05)$ 203values2_i:makelist(subst(x=i,data2_cont),i,time_i)$ 204/* [wxMaxima: input end ] */ 205 206 207/* [wxMaxima: comment start ] 208Now we generate a matrix of values lsquares can deal with: 209 [wxMaxima: comment end ] */ 210 211 212/* [wxMaxima: input start ] */ 213data2_i:float(transpose(matrix(time_i,values2_i)))$ 214/* [wxMaxima: input end ] */ 215 216 217/* [wxMaxima: comment start ] 218The result is a still noisy and distorted, but more evenly-spaced curve: 219 [wxMaxima: comment end ] */ 220 221 222/* [wxMaxima: input start ] */ 223wxdraw2d( 224 points(data2_i), 225 grid=true,xlabel="x",ylabel="y" 226)$ 227/* [wxMaxima: input end ] */ 228 229 230/* [wxMaxima: comment start ] 231Fitting this curve will yield a better result as the first attempt: 232 [wxMaxima: comment end ] */ 233 234 235/* [wxMaxima: input start ] */ 236lsquares_estimates_approximate( 237 lsquares_mse( 238 data2_i,[x,y],approach1 239 ), 240 [a,b,c] 241); 242params2_2:%[1]; 243wxdraw2d( 244 grid=true,xlabel="x",ylabel="y", 245 color=red,key="Reconstructed", 246 explicit(rhs(subst(params2_2,approach1)),x,-10,10), 247 color=blue,key="Original", 248 explicit(fun1,x,-10,10) 249)$ 250/* [wxMaxima: input end ] */ 251 252 253/* [wxMaxima: section start ] 254Fitting data to a·exp(k·t) 255 [wxMaxima: section end ] */ 256 257 258/* [wxMaxima: comment start ] 259The natural approach to fitting data to exponential curves would be: 260 [wxMaxima: comment end ] */ 261 262 263/* [wxMaxima: input start ] */ 264approach2:y=a*exp(k*t); 265/* [wxMaxima: input end ] */ 266 267 268/* [wxMaxima: comment start ] 269Unfortunately it is hard to fit to experimental data to this approach in many ways, for example: 270 * One part of this curve results in low values. Even small measurement noise on this part will yield widely incorrect results for the curve parameters as the fitter is trying to model the noise, too, and as noise in nonlinear systems tends not to be averanged out completely. 271 * Another part of this curve contains quite in high values and is sensitive to even small changes in k. As the fitter wants to keep the overall error low it will therefore respect this part much more than the lower, more good-natured part of the curve. 272 * Starting from the wrong point and optimizing a and k in the direction that reduces the error most in each step might lead to a point that is far from being the solution 273 * and trying to use lsquares_estimates to find the ideal solution often leads to numbers that exceed the floating-point range (or exact numbers longer than the computer's memory). 274 [wxMaxima: comment end ] */ 275 276 277/* [wxMaxima: comment start ] 278A better approach is therefore to use Caruana's approach for fitting: 279 [wxMaxima: comment end ] */ 280 281 282/* [wxMaxima: input start ] */ 283approach2/a; 284caruana:log(%); 285/* [wxMaxima: input end ] */ 286 287 288/* [wxMaxima: comment start ] 289This approach is much more good-natured for finding the parameters (that can then be substituted into approach2). 290 [wxMaxima: comment end ] */ 291 292 293/* [wxMaxima: comment start ] 294The problem with noise causing finding erroneous parameters is still valid, though, in this case. It can be partially eliminated by introducing a c_noise, as proposed by Guo. 295 [wxMaxima: comment end ] */ 296 297 298/* [wxMaxima: input start ] */ 299lhs(approach2)=rhs(approach2+c_noise); 300%-c_noise; 301%/a; 302approach_guo:log(%),logexpand=super; 303/* [wxMaxima: input end ] */ 304 305 306/* [wxMaxima: comment start ] 307Additionally https://ieeexplore.ieee.org/document/5999593 offers an iterative method that allows to reduce the influence of noise in each step: 308 [wxMaxima: comment end ] */ 309 310 311/* [wxMaxima: section start ] 312Using a different fitting algorithm 313 [wxMaxima: section end ] */ 314 315 316/* [wxMaxima: comment start ] 317Maxima provides a second algorithm that often produces even better results, but is a bit more complicated to use. For example it requires us to manually compile a list of the errors we want to minimize. Let's do that for data2_i and approach1: 318 [wxMaxima: comment end ] */ 319 320 321/* [wxMaxima: input start ] */ 322approach1; 323/* [wxMaxima: input end ] */ 324 325 326/* [wxMaxima: input start ] */ 327errval:lhs(approach1)-rhs(approach1); 328my_mse:makelist( 329 subst( 330 [x=i[1],y=i[2]], 331 errval 332 ), 333 i,args(data2_i) 334)$ 335/* [wxMaxima: input end ] */ 336 337 338/* [wxMaxima: comment start ] 339Hints for understanding this construct: 340 * args() converts a matrix to a list of lists, 341 * The makelist() command steps through this list and for each data point assigns the list of the x and y value to i 342 * rhs() and lhs() extract the right hand side and the left hand side of an equation (the part right and left of the "=") 343 * and since we know which element in i means which variable we can use subst() in orer to substitute the elements in i into errval, the equation that tells us how big the error in this point is. 344 [wxMaxima: comment end ] */ 345 346 347/* [wxMaxima: comment start ] 348The result is a list of error values, each in the following format: 349 [wxMaxima: comment end ] */ 350 351 352/* [wxMaxima: input start ] */ 353my_mse[1]; 354/* [wxMaxima: input end ] */ 355 356 357/* [wxMaxima: comment start ] 358Now let's load the package that contains the other fitter: 359 [wxMaxima: comment end ] */ 360 361 362/* [wxMaxima: input start ] */ 363load("minpack")$ 364/* [wxMaxima: input end ] */ 365 366 367/* [wxMaxima: comment start ] 368Fitting the data is simple: 369 [wxMaxima: comment end ] */ 370 371 372/* [wxMaxima: input start ] */ 373param_list:[a,b,c]; 374result:minpack_lsquares(my_mse,param_list,[1,1,1]); 375/* [wxMaxima: input end ] */ 376 377 378/* [wxMaxima: comment start ] 379Converting the result into a list of equations is a bit more complicated: 380 [wxMaxima: comment end ] */ 381 382 383/* [wxMaxima: input start ] */ 384params_3:map(lambda([x],x[1]=x[2]),args(transpose(matrix(param_list,result[1])))); 385/* [wxMaxima: input end ] */ 386 387 388/* [wxMaxima: comment start ] 389Hints for understanding this line: 390 * Map runs a command on each element of a list 391 * We want to provide map with a command that converts something of the type "[a,3]" to an "a=3". But we won't re-use this command so there is no need to actually give this command a name. Therefore we use lambda() in order to create a name-less command with one parameter, x. 392 * args again converts a matrix to a list of lists and 393 * we need to convert our input data into a list in order to transpose it (which means: exchange the columns by rows and vice versa) 394 [wxMaxima: comment end ] */ 395 396 397/* [wxMaxima: comment start ] 398The result isn't too bad, this time, neither: 399 [wxMaxima: comment end ] */ 400 401 402/* [wxMaxima: input start ] */ 403wxdraw2d( 404 grid=true,xlabel="x",ylabel="y", 405 color=red,key="Reconstructed", 406 explicit(rhs(subst(params_3,approach1)),x,-10,10), 407 color=blue,key="Original", 408 explicit(fun1,x,-10,10) 409)$ 410/* [wxMaxima: input end ] */ 411 412 413 414/* Old versions of Maxima abort on loading files that end in a comment. */ 415"Created with wxMaxima 20.03.1-DevelopmentSnapshot"$ 416