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