064a58035ab6afc487ff4133a3d2204fefc5e4ff
[scilab.git] / scilab / modules / optimization / help / en_US / nonlinearleastsquares / datafit.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - INRIA
5  * Copyright (C) 2011 - DIGITEO - Michael Baudin
6  * Copyright (C) 2018 - Samuel GOUGEON
7  *
8  * Copyright (C) 2012 - 2016 - Scilab Enterprises
9  *
10  * This file is hereby licensed under the terms of the GNU GPL v2.0,
11  * pursuant to article 5.3.4 of the CeCILL v.2.1.
12  * This file was originally licensed under the terms of the CeCILL v2.1,
13  * and continues to be available under such terms.
14  * For more information, see the COPYING file which you should have received
15  * along with this program.
16  *
17  -->
18 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
19         xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml"
20         xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
21         xmlns:scilab="http://www.scilab.org" xml:id="datafit" xml:lang="en">
22     <refnamediv>
23         <refname>datafit</refname>
24         <refpurpose>Non linear (constrained) parametric fit of measured (weighted) data</refpurpose>
25     </refnamediv>
26     <refsynopsisdiv>
27         <title>Syntax</title>
28         <synopsis>[p, fmin, status] = datafit([iprint,] G [,DG], Data [,Wd] [,Wg][,'b',p_inf,p_sup], p0 [,algo][,stop])
29         </synopsis>
30     </refsynopsisdiv>
31     <refsection>
32         <title>Arguments</title>
33         <variablelist>
34             <varlistentry>
35                 <term>iprint</term>
36                 <listitem>
37                     <para>scalar used to set the trace mode:
38                         <table>
39                             <tr>
40                                 <th>Value</th><th>Action</th>
41                             </tr>
42                             <tr>
43                                 <td>0</td><td>nothing (except errors) is reported</td>
44                             </tr>
45                             <tr>
46                                 <td>1</td><td>initial and final reports</td>
47                             </tr>
48                             <tr>
49                                 <td>2</td><td>adds a report per iteration</td>
50                             </tr>
51                             <tr>
52                                 <td>>2</td><td>adds reports on linear search</td>
53                             </tr>
54                         </table>
55                         <warning>
56                             Most of the reports are displayed in the Scilab console.
57                         </warning>
58                     </para>
59                 </listitem>
60             </varlistentry>
61             <varlistentry>
62                 <term>p0</term>
63                 <listitem>
64                     Matrix of initial guess of the model parameters.
65                   <para/>
66                 </listitem>
67             </varlistentry>
68             <varlistentry>
69                 <term>"b"</term>
70                 <listitem>
71                     keyword heading the <literal>p_inf, p_sup</literal> sequence, where
72                     <literal>p_inf</literal> and <literal>p_sup</literal> are real matrices
73                     with the <literal>p0</literal> shape and sizes.
74                     <literal>p_inf</literal> is the vector of lower bounds of respective
75                     <literal>p</literal> components. <literal>p_sup</literal> are upper bounds.
76                     <para/>
77                 </listitem>
78             </varlistentry>
79             <varlistentry>
80                 <term>p</term>
81                 <listitem>
82                     Matrix of the best fitting model parameters, with the <varname>p0</varname> sizes.
83                     <para/>
84                 </listitem>
85             </varlistentry>
86             <varlistentry>
87                 <term>Data</term>
88                 <listitem>
89                     (ndc x nd) matrix of experimental data points, where <literal>Data(:,i)</literal>
90                     are the #nzc coordinates of the ith measurement.
91                     <para/>
92                 </listitem>
93             </varlistentry>
94             <varlistentry>
95                 <term>Wd</term>
96                 <listitem>
97                     Data weights (optional). Row vector with size(Data,2) elements.
98                     <para/>
99                 </listitem>
100             </varlistentry>
101             <varlistentry>
102                 <term>G</term>
103                 <listitem>
104                     <para>
105                         Gap function: handle of a function in Scilab language or of a compiled
106                         external function. Computed gap(s) between the fitting model and data
107                         may be of the first order, not quadratic.
108                     </para>
109                     <para>
110                         <emphasis role="bold">
111                             If <literal>G()</literal>
112                         </emphasis>
113                         does not need
114                         any input parameters (like constants), only arguments, its syntax must be
115                         <literal>g = G(p, Data)</literal>.
116                     </para>
117                     <para>
118                         <emphasis role="bold">
119                             If <literal>G()</literal>
120                         </emphasis>
121                         needs
122                         some input parameters in addition to input arguments, its syntax must be
123                         <literal>g = G(p, Data, param1, param2, .., paramN)</literal>.
124                         Then, <literal>list(G, param1, param2, .., paramN)</literal> must be provided
125                         to <literal>datafit()</literal> instead of only <literal>G</literal>.
126                     </para>
127                     <para>
128                         The following definitions will be used:
129                         <itemizedlist>
130                             <listitem>
131                                 <literal>ng = size(g,1)</literal>: number of fitting criteria (= of
132                                 gap definitions).
133                             </listitem>
134                             <listitem>
135                                 <literal>np = size(p,"*")</literal>: total number of fitting parameters.
136                             </listitem>
137                             <listitem>
138                                 <literal>nd = size(Data,2)</literal>: number of data points.
139                             </listitem>
140                             <listitem>
141                                 <literal>ndc = size(Data,1)</literal>: number of coordinates of each
142                                 data point.
143                             </listitem>
144                         </itemizedlist>
145                     </para>
146                     <para>
147                         Then, the sizes of G arguments are the following ones:
148                         <table>
149                             <tr valign="top">
150                                 <th>g</th>
151                                 <td>: (ng x nd) matrix
152                                     <itemizedlist>
153                                         <listitem>
154                                             Each column is related to a given data point and retuns its
155                                             gaps to the fitting model.
156                                         </listitem>
157                                         <listitem>
158                                             Each row is about a gap definition. <literal>G</literal> may
159                                             compute several types of gaps. For instance, the vertical
160                                             gaps may be returned in <literal>g(1,:)</literal> ; the
161                                             horizontal gaps in <literal>g(2,:)</literal> ; the nearest
162                                             distances in <literal>g(3,:)</literal> ; etc.
163                                         </listitem>
164                                     </itemizedlist>
165                                 </td>
166                             </tr>
167                             <tr valign="top">
168                                 <th>p</th>
169                                 <td>: fitting parameters to be optimized.
170                                     <para/>
171                                 </td>
172                             </tr>
173                             <tr valign="top">
174                                 <th>Data</th>
175                                 <td>: (ndc x nd) matrix: Each column provides the
176                                     <literal>#ndc</literal> coordinates of a given data point.
177                                     If <literal>G</literal> computes gaps only for a single
178                                     data point, <literal>datafit()</literal> will automatically
179                                     loop over the set of <literal>Data</literal> columns.
180                                     This makes <literal>datafit()</literal> slower.
181                                 </td>
182                             </tr>
183                         </table>
184                     </para>
185                 </listitem>
186             </varlistentry>
187             <varlistentry>
188                 <term>DG</term>
189                 <listitem>
190                     (optional) Handle of a function in Scilab language or of a compiled external
191                     function, computing the partial of <varname>G</varname> wrt <varname>p</varname>.
192                     Its syntax must be like <literal>S = DG(p, Data)</literal>. The expected sizes
193                     of the result are (ng x np x 1) (Data-unvectorized) or (ng x np x nd)
194                     (Data-vectorized).
195                     <para/>
196                 </listitem>
197             </varlistentry>
198             <varlistentry>
199                 <term>Wg</term>
200                 <listitem>
201                     Optional (ng x ng) matrix of weights of gaps definitions. Default = eye().
202                     <varname>G()</varname> defines as many types of gaps as required.
203                     Most often, a unique gap will be computed -- most often the vertical one --,
204                     but other ones can be computed (horizontal, nearest distance, etc).
205                     This matrix allows to weight the different types.
206                     <para/>
207                 </listitem>
208             </varlistentry>
209             <varlistentry>
210                 <term>algo</term>
211                 <listitem>
212                     Word <literal>'qn'</literal> (quasi-Newton, default), or
213                     <literal>'gc'</literal> (conjugate gradient), or
214                     <literal>'nd'</literal> (non-differentiable. Then <varname>p_inf</varname>
215                     and <varname>p_sup</varname> parameters bounds are not accepted).
216                     Selects the algorithm used to perform the fitting.
217                     See <function>optim()</function> for more details.
218                     <para/>
219                 </listitem>
220             </varlistentry>
221             <varlistentry>
222                 <term>stop</term>
223                 <listitem>
224                     <para>sequence of optional arguments controlling the convergence of
225                         the fitting algorithm:
226                     </para>
227                     <para>
228                         <literal>'ar',nap, [iter [,epsg [,epsf [,epsx]]]]</literal>, with
229                     </para>
230                     <table>
231                         <tr>
232                             <th>"ar"</th>
233                             <td>
234                                 reserved keyword heading the list of stopping criteria as defined below.
235                             </td>
236                         </tr>
237                         <tr>
238                             <th>nap</th>
239                             <td>
240                                 maximum number of calls to the cost f function allowed (default: 100).
241                             </td>
242                         </tr>
243                         <tr>
244                             <th>iter</th>
245                             <td>maximum number of iterations allowed (default: 100).</td>
246                         </tr>
247                         <tr>
248                             <th>epsg</th>
249                             <td>threshold on gradient norm (default: %eps).</td>
250                         </tr>
251                         <tr>
252                             <th>epsf</th>
253                             <td>
254                                 threshold controlling the decreasing of <literal>f</literal>
255                                 (defined in the description section). Default: 0.
256                             </td>
257                         </tr>
258                         <tr>
259                             <th>epsx</th>
260                             <td>
261                                 threshold controlling the variation of <varname>p</varname>.
262                                 This vector (possibly matrix) of same size as <varname>p0</varname>
263                                 can be used to scale <varname>p</varname>. Default: 0.
264                             </td>
265                         </tr>
266                     </table>
267                 </listitem>
268             </varlistentry>
269             <!--            <varlistentry>
270                 <term>fmin</term>
271                 <listitem>
272                     <para>
273                         <literal>f</literal> quadratic residu (minimal <literal>f</literal> value reached).
274                         <literal>sqrt(fmin/ng/nz)</literal> or
275                         <literal>sqrt(fmin/ng/sum(Wd))</literal> estimates the Model-to-Data average
276                         distance.
277                     </para>
278                 </listitem>
279             </varlistentry>
280 -->
281             <varlistentry>
282                 <term>dmin</term>
283                 <listitem>
284                     Average Model-to-data residual distance. It is equal to
285                     <literal>sqrt(fmin/ng/nd)</literal> or
286                     <literal>sqrt(fmin/ng/sum(Wd))</literal>, where fmin is the quadratic residue
287                     of the cost function <literal>f</literal> (see the description).
288                     <para/>
289                 </listitem>
290             </varlistentry>
291             <varlistentry>
292                 <term>status</term>
293                 <listitem>
294                   return flag = termination status (only for the 'qn' default algorithm.
295                   Set to %nan otherwise): <literal>9</literal> is successful. Possible values are:
296                   <table>
297                     <tr valign="top">
298                       <th>1</th><td>: Norm of projected gradient lower than...</td>
299                       <td>....</td>
300                       <th>6</th><td>: Optim stops: too small variations in gradient direction.</td>
301                     </tr>
302                     <tr valign="top">
303                       <th>2</th><td>: At last iteration, f decreases by less than...</td>
304                       <td></td>
305                       <th>7</th><td>: Stop during calculation of descent direction.</td>
306                     </tr>
307                     <tr valign="top">
308                       <th>3</th><td>: Optim stops: Too small variations for p.</td>
309                       <td></td>
310                       <th>8</th><td>: Stop during calculation of estimated hessian.</td>
311                     </tr>
312                     <tr valign="top">
313                       <th>4</th><td>: Optim stops: maximum number of calls to f is reached.</td>
314                       <td></td>
315                       <th>9</th><td>: Successful completion.</td>
316                     </tr>
317                     <tr valign="top">
318                       <th>5</th><td>: Optim stops: maximum number of iterations is reached.</td>
319                       <td></td>
320                       <th>10</th><td>: Linear search fails.</td>
321                     </tr>
322                   </table>
323                   <para/>
324                 </listitem>
325             </varlistentry>
326         </variablelist>
327         <para>
328             Other <function>optim()</function> input arguments such that <varname>df0</varname>,
329             <varname>mem</varname>, or <varname>'in'</varname> may be used with
330             <literal>datafit()</literal>. They will be passed as is to <function>optim()</function>.
331         </para>
332     </refsection>
333     <refsection>
334         <title>Description</title>
335         <para>
336             <literal>datafit</literal> is used to fit a parametrized model to given data.
337         </para>
338         <para>
339             A function <literal>G(p, Data)</literal> must be defined to compute the gaps between the
340             <varname>Data</varname> points and a model whose parameters to be tuned are provided
341             through the matrix <varname>p</varname>.
342
343         </para>
344         <para>
345             <literal>datafit()</literal> checks whether <literal>G()</literal> is vectorized or not
346             over <varname>Data</varname>. If <literal>G()</literal> works only with a single
347             <literal>data</literal> point as a single column of coordinates, then <literal>datafit()</literal>
348             loops over <literal>data=Data(:,i)</literal> points. Otherwise, <literal>G(p, Data)</literal>
349             is called for the whole <varname>Data</varname> array of data coordinates, which is
350             way faster if <varname>G()</varname>'s code is truly vectorized.
351         </para>
352         <para>
353             Starting from initial values <varname>p0</varname> of the model parameters,
354             <literal>datafit()</literal> varies them in order to minimize the whole Model-to-Data
355             distance defined as
356             <screen>
357 f = (G(p, Data(:,1))' * Wg * G(p, Data(:,1)))* Wd(1) + ..
358     (G(p, Data(:,2))' * Wg * G(p, Data(:,2)))* Wd(2) + ..
359     ...
360     (G(p, Data(:,nz))' * Wg * G(p,Data(:,nz)))* Wd(nz)
361 </screen>
362         </para>
363         <para>
364             When <literal>Wg = eye()</literal> (default), this definition is equivalent to
365 <screen>
366 f = sum(G(p, Data(:,1)).^2) * Wd(1) + ..
367     sum(G(p, Data(:,2)).^2) * Wd(2) + ..
368     ...
369     sum(G(p, Data(:,nz)).^2) * Wd(nz)
370 </screen>
371         </para>
372         <para>
373             If in addition <varname>G</varname> returns only one gap type, this even simplifies to
374             <screen>
375 f = G(p, Data(:,1))^2 * Wd(1) + ..
376     G(p, Data(:,2))^2 * Wd(2) + ..
377     ...
378     G(p, Data(:,nz))^2 * Wd(nz)
379 </screen>
380         </para>
381         <para>
382             <literal>datafit()</literal> calls <function>optim()</function> to minimize f.
383         </para>
384         <para>
385             Choosing appropriately <varname>p0</varname> may improve and faster
386             <literal>datafit()</literal>'s convergence to the best fit.
387         </para>
388     </refsection>
389     <refsection>
390         <title>Examples</title>
391         <table>
392             <tr>
393                 <td>
394                     <para>
395                         <emphasis role="bold">Simple example: Polynomial fitting (parabola = 3 parameters)</emphasis>.
396                     </para>
397                     <programlisting role="example"><![CDATA[
398 // Model (parabola)
399 function y = FF(x,p)
400     y = (p(1)*x + p(2)).*x + p(3)
401 endfunction
402
403 // Data creation
404 pa = [2;-4;14] // parameters of the actual law, to generate data
405 X = 0:0.1:3.5;
406 Y0 = FF(X, pa);
407 noise = 10*(rand(Y0)-.5);   // We add some noise
408 Y = Y0 + noise
409 Data = [X ; Y];
410
411 // Gap function definition
412 function e = G(p, Data)
413     e = Data(2,:) - FF(Data(1,:), p);  // vertical distance
414 endfunction
415
416 // Performing the fitting
417 // ----------------------
418 // Without weighting data:
419 p0 = [1;-1;10]
420 [p, dmin] = datafit(G, Data, p0);
421 // The uncertainty is identified to the noise value.
422 // The weight is set inversely proportional to the uncertainty
423 [pw, dmin] = datafit(G, Data, 1./abs(noise), p0);
424
425 // Display
426 // -------
427 scf(0);
428 clf()
429 xsetech([0 0 1 0.83])
430 plot2d(X, Y0, 24)       // True underlaying law
431 plot2d(X, Y, -1)        // Noisy data
432 plot2d(X, FF(X, p), 15)  // best unweighted fit
433 plot2d(X, FF(X, pw), 18) // best weighted fit
434 legend(["True law" "Noisy data" "Unweighted fit" "Weighted fit"],2)
435
436 xsetech([0 .75 1 0.25])
437 plot2d(X, FF(X,p)-Y0, 15)  // remaining gap <= best unweighted fit
438 plot2d(X, FF(X,pw)-Y0, 18) // remaining gap <= best weighted fit
439 ylabel("Data fit - True law")
440 ax = gca(); ax.x_location = "top";
441 tmp = ax.x_ticks.labels;
442 gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
443 xgrid(color("grey70"))
444          ]]></programlisting>
445                 </td>
446                 <td>
447                     <scilab:image>
448                         // Model (parabola)
449                         function y = FF(x,p)
450                         y = p(1)*x.*x + p(2)*x + p(3)
451                         endfunction
452
453                         // Data creation
454                         pa = [2;-4;14] // parameters of the actual law, used to generate data
455                         X = 0:0.1:3.5;
456                         Y0 = FF(X, pa);
457                         noise = 10*(rand(Y0)-.5);   // We add some noise
458                         Y = Y0 + noise
459                         Data = [X ; Y];
460
461                         // Gap function definition
462                         function e = G(p, Data)
463                         e = Data(2,:) - FF(Data(1,:), p);  // vertical distance
464                         endfunction
465
466                         // Performing the fitting
467                         // ----------------------
468                         // Without weighting data:
469                         p0 = [1;-1;10]
470                         [p, dmin] = datafit(G, Data, p0);
471
472                         // The uncertainty is identified to the noise value.
473                         // The weight is set inversely proportional to the uncertainty
474                         [pw, dmin] = datafit(G, Data, 1 ./abs(noise), p0);
475
476                         // Display
477                         // -------
478                         scf(0);
479                         clf()
480                         xsetech([0 0 1 0.83])
481                         plot2d(X, Y0, 24)       // True underlaying law
482                         plot2d(X, Y, -1)        // Noisy data
483                         plot2d(X, FF(X,p), 15)  // best unweighted fit
484                         plot2d(X, FF(X,pw), 18) // best weighted fit
485                         legend(["True law" "Noisy data" "Unweighted fit" "Weighted fit"],2)
486
487                         xsetech([0 .75 1 0.25])
488                         plot2d(X, FF(X,p)-Y0, 15)  // remaining gap with best unweighted fit
489                         plot2d(X, FF(X,pw)-Y0, 18) // remaining gap best weighted fit
490                         ylabel("Data fit - True law")
491                         ax = gca(); ax.x_location = "top";
492                         tmp = ax.x_ticks.labels;
493                         gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
494                         xgrid(color("grey70"))
495                         // ===============
496                         gcf().axes_size = [500 500];
497                         gcf().children(1).margins(1) = 0.07;
498                         gcf().children(2).margins(1) = 0.07;
499                     </scilab:image>
500                 </td>
501             </tr>
502             <tr>
503                 <td>
504                     <para/>
505                     <para>
506                         <emphasis role="bold">Fitting a gaussian curve + sloping base (5 parameters)</emphasis>:
507                     </para>
508                     <programlisting role="example"><![CDATA[
509 function Y = model(P, X)
510     Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
511 endfunction
512 // -------------------------------
513 // Gap function
514 function r = gap(P,XY)
515     r = XY(2,:) - model(P,XY(1,:))
516 endfunction
517 // -------------------------------
518 // True parameters
519 Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
520
521 // Generating data
522 // ---------------
523 X = 0:0.2:10;
524 Np = length(X);
525 Y0 = model(Pg, X);              // True law
526 noise = (rand(Y0)-.5)*Pg(3)/4;
527 Y = Y0 + noise                  // Noised data
528 Data = [X ; Y];
529
530 // Trying to recover original parameters from the (noisy) data:
531 // -----------------------------------------------------------
532 // Performing the fitting = parameters optimization
533 [v, k] = max(Y0);
534 P0 = [X(k) 1 v 1 1];
535 [P, dmin, s]  = datafit(gap, Data, P0);
536 Wd =  1./abs(noise);
537 [Pw, dminw, s] = datafit(gap, Data, Wd, P0);
538
539 // Computing fitting curves from found parameters
540 fY = model(P, X);
541 fYw = model(Pw, X);
542
543 // Display
544 // -------
545 scf(0);
546 clf()
547 xsetech([0 0 1 0.8])
548 plot2d(X, Y0, 24)  // True underlaying law
549 plot2d(X, Y,  -1)  // Noisy data
550 plot2d(X, fY, 15)  // best unweighted fit
551 plot2d(X, fYw,18)  // best weighted fit
552
553 // Legends: Average vertical Model-Data distance:
554 // True one v/s according to the residual cost
555 L_1 = "Unweighted fit (%5.2f /%5.2f)";
556 L_1 = msprintf(L1, mean(abs(fY-Y0)), sqrt(dmin/Np));
557 L_2 = "Weighted fit (%5.2f /%5.2f)";
558 L_2 = msprintf(L1, mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)));
559 legend(["True law", "Noisy Data", L_1, L_2],1)
560
561 // Params: row#1: true params row#2:
562 //  from unweighted fit row#3: from weighted fit
563 [xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
564 xnumb(xs, ys, [Pg ; P ; Pw]);
565 LawTxt = "$\Large p_3\,.\,"+..
566          "exp\left({-(x-p_1)^2}\over{2\,p_2} \right)+p_4.x+p_5$";
567 xstring(2, 1.3, LawTxt)
568
569 // Plotting residus:
570 xsetech([0 .75 1 0.25])
571 plot2d(X, model(P ,X)-Y0, 15) // remaining gap with best unweighted fit
572 plot2d(X, model(Pw,X)-Y0, 18) // remaining gap best weighted fit
573 ylabel("Data fit - True law")
574 ax = gca();
575 ax.x_location = "top";
576 tmp = ax.x_ticks.labels;
577 gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
578 xgrid(color("grey70"))
579          ]]></programlisting>
580                 </td>
581                 <td>
582                     <scilab:image>
583                         function Y = model(P,X)
584                         Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
585                         endfunction
586                         // -------------------------------
587                         // Gap function
588                         function r = gap(P,XY)
589                         r = XY(2,:) - model(P,XY(1,:))
590                         endfunction
591                         // -------------------------------
592                         // True parameters
593                         Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
594
595                         // Generating data
596                         // ---------------
597                         X = 0:0.2:10;
598                         Np = length(X);
599                         Y0 = model(Pg, X);              // True law
600                         noise = (rand(Y0)-.5)*Pg(3)/4;
601                         Y = Y0 + noise                  // Noised data
602                         Data = [X ; Y];
603
604                         // Let's try to find the original parameters from the (noisy) data:
605                         // ---------------------------------------------------------------
606                         // Performing the fitting = parameters optimization
607                         [v, k] = max(Y0);
608                         P0 = [X(k) 1 v 1 1];
609                         [P, dmin, s]  = datafit(gap, Data, P0);
610                         Wd =  1./abs(noise);
611                         [Pw, dminw, s] = datafit(gap, Data, Wd, P0);
612
613                         // Computing fitting curves from found parameters
614                         fY = model(P, X);
615                         fYw = model(Pw,X);
616
617                         // Display
618                         // -------
619                         scf(0);
620                         clf()
621                         xsetech([0 0 1 0.82])
622                         plot2d(X, Y0, 24)  // True underlaying law
623                         plot2d(X, Y,  -1)  // Noisy data
624                         plot2d(X, fY, 15)  // best unweighted fit
625                         plot2d(X, fYw,18)  // best weighted fit
626                         legend(["True law" "Noisy Data" ..
627                         msprintf("Unweighted fit (%5.2f /%5.2f)",mean(abs(fY-Y0)), sqrt(dmin/Np)) ..
628                         msprintf("Weighted fit (%5.2f /%5.2f)", mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)))],1)
629                         // Average vertical Model-Data distance: True one / according to the residual cost
630                         // Params: row#1: true params row#2: from unweighted fit row#3: from weighted fit
631                         [xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
632                         xnumb(xs, ys, [Pg ; P ; Pw])
633                         xstring(2,1.3,"$\Large p_3\,.\,exp\left({-(x-p_1)^2}\over{2\,p_2} \right)+p_4.x+p_5$")
634
635                         // Plotting residus:
636                         xsetech([0 .75 1 0.25])
637                         plot2d(X, model(P ,X)-Y0, 15)  // remaining gap with best unweighted fit
638                         plot2d(X, model(Pw,X)-Y0, 18)  // remaining gap best weighted fit
639                         ylabel("Data fit - True law")
640                         ax = gca();
641                         ax.x_location = "top";
642                         tmp = ax.x_ticks.labels;
643                         gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
644                         xgrid(color("grey70"))
645                         // ===============
646                         gcf().axes_size = [500 600];
647                     </scilab:image>
648                 </td>
649             </tr>
650             <tr>
651                 <td>
652                     <para/>
653                     <para>
654                         <emphasis role="bold">Extraction of 3 overlapping gaussian curves (9 parameters)</emphasis>:
655                         example with constrained bounds and a matrix of parameters.
656                     </para>
657                     <programlisting role="example"><![CDATA[
658 // Basic gaussian component:
659 function Y = gauss(X, average, stDev, ymax)
660     Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
661 endfunction
662
663 // Model: cumulated gaussian laws
664 function r = model(P, X)
665     // P(1,:): averages  P(2,:): standard deviation  P(3,:): ymax
666     r = 0;
667     for i = 1:size(P,2)
668         r = r + gauss(X, P(1,i), P(2,i), P(3,i));
669     end
670 endfunction
671
672 // Gap function:
673 function r = gap(P,XY)
674     r = XY(2,:) - model(P,XY(1,:))
675 endfunction
676
677 // Actual parameters :
678 Pa = [ 1.1  4    7    // averages
679         1  0.8  1.5   // standard deviations
680         3   2   5.5   // ymax
681       ];
682 Nc = size(Pa,2);      // Number of summed curves
683 Np = 200;             // Number of data points
684
685 // Data generation (without noise)
686 xmin = 0;
687 xmax = 10;
688 X = linspace(xmin-2, xmax+2, Np);
689 Y = model(Pa, X);
690 Data = [X ; Y];          // Curve to analyze / fit
691
692 // FITTING
693 // -------
694 // Fitting parameters:
695 // Initial values & bounds (amplitudes must be > 0)
696 Po = [ linspace(xmin,xmax,Nc)
697        ones(1,Nc)*0.5
698        ones(1,Nc)*max(Y)/2 ];
699 Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
700 Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3   ; ones(1,Nc)*max(Y)*1.2];
701 // Fitting
702 [P, dmin] = datafit(gap,Data,'b',Pinf,Psup,Po,'qn','ar',200,200);
703
704 // Display
705 // -------
706 clf()
707 xsetech([0 0 1 0.8]);
708 plot(X, Y, "-b", X, model(P,X), "-r");
709 gca().children.children(1:2).thickness = 2;
710 for i = 1:Nc
711     // True gaussian components
712     plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c");
713     // Extracted gaussian components
714     plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g");
715 end
716 legend("Data", "Parametric Model", ..
717        "True components", "Extracted components",  2);
718
719 xsetech([0 0.75 1 0.25])
720 plot(X, model(P,X)-Y);
721 gca().x_location = "top";
722 gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
723 ylabel("Model - Data")
724 xgrid(color("grey70"))
725          ]]></programlisting>
726                 </td>
727                 <td>
728                     <scilab:image>
729                         // Basic gaussian component:
730                         function Y = gauss(X, average, stDev, ymax)
731                         Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
732                         endfunction
733
734                         // Model: cumulated gaussian laws
735                         function r = model(P, X)
736                         // P(1,:): averages  P(2,:): standard deviation  P(3,:): ymax
737                         r = 0;
738                         for i = 1:size(P,2)
739                         r = r + gauss(X, P(1,i), P(2,i), P(3,i));
740                         end
741                         endfunction
742
743                         // Gap function:
744                         function r = gap(P,XY)
745                         r = XY(2,:) - model(P,XY(1,:))
746                         endfunction
747
748                         // True parameters :
749                         Pa = [ 1.1  4    7    // averages
750                         1  0.8  1.5   // standard deviations
751                         3   2   5.5   // ymax
752                         ];
753                         Nc = size(Pa,2);       // Number of summed curves
754                         Np = 200;             // Number of data points
755
756                         // Data generation (without noise)
757                         xmin = 0;
758                         xmax = 10;
759                         X = linspace(xmin-2, xmax+2, Np);
760                         Y = model(Pa, X);
761                         Data = [X ; Y];          // Curve to analyze / fit
762
763                         // FITTING
764                         // -------
765                         // Fitting parameters: initial values and bounds (amplitudes must be > 0)
766                         Po = [ linspace(xmin,xmax,Nc) ; ones(1,Nc)*0.5 ; ones(1,Nc)*max(Y)/2];
767                         Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
768                         Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3   ; ones(1,Nc)*max(Y)*1.2];
769                         // Fitting
770                         [P, dmin] = datafit(gap, Data, 'b',Pinf,Psup, Po, 'qn','ar',200,200);
771
772                         // Display
773                         // -------
774                         clf()
775                         xsetech([0 0 1 0.8]);
776                         plot(X, Y, "-b", X, model(P,X), "-r");
777                         gca().children.children(1:2).thickness = 2;
778                         for i = 1:Nc
779                         plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c"); // True gaussian components
780                         plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g");    // Extracted gaussian components
781                         end
782                         legend("Data", "Parametric Model","True components", "Extracted components",  2);
783
784                         xsetech([0 0.75 1 0.25])
785                         plot(X, model(P,X)-Y);
786                         gca().x_location = "top";
787                         gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
788                         ylabel("Model - Data")
789                         xgrid(color("grey70"))
790                         // ===============
791                         gcf().axes_size = [500 500];
792                     </scilab:image>
793                 </td>
794             </tr>
795             <tr>
796                 <td>
797                     <para>
798                         <emphasis role="bold">With a parametric curve: Fitting a off-centered tilted 2D elliptical arc (5 parameters)</emphasis>.
799                     </para>
800                     <programlisting role="example"><![CDATA[
801 function g = Gap(p, Data)
802     // p: a, b, center xc, yc, alpha tilt °
803     C = cosd(p(5)); S = -sind(p(5));
804     x = Data(1,:) - p(3);
805     y = Data(2,:) - p(4);
806     g = ((x*C - y*S )/p(1)).^2 + ((x*S + y*C)/p(2)).^2 - 1;
807 endfunction
808
809 // Generating the data
810 // -------------------
811 // Actual parameters :
812 Pa = [3, 1.3, -2, 1.5, 30];
813 Np = 30;            // Number of data points
814 Theta = grand(1,Np, "unf",-100, 170);
815 // untilted centered noisy arc
816 x = Pa(1)*(cosd(Theta) + grand(Theta, "unf", -0.07, 0.07));
817 y = Pa(2)*(sind(Theta) + grand(Theta, "unf", -0.07, 0.07));
818 // Tilting and off-centering the arc
819 A = Pa(5);
820 C = cosd(A); S = sind(A);
821 xe =  C*x + y*S + Pa(3);
822 ye = -S*x + y*C + Pa(4);
823 Data = [xe ; ye];
824
825 // Retrieving parameters from noisy data
826 // -------------------------------------
827 // Initial estimation
828 ab0 = (strange(xe)+strange(ye))/2;
829 xc0 = mean(xe);
830 yc0 = mean(ye);
831 A0 = -atand(reglin(xe,ye));
832 P0 = [ab0 ab0/2 xc0 yc0 A0];
833 // Parameters bounds
834 Pinf = [ 0     0     xc0-2*ab0, yc0-2*ab0 -90];
835 Psup = [3*ab0 3*ab0  xc0+2*ab0, yc0+2*ab0  90];// Fitting
836
837 // FITTING
838 [P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
839
840 // P(1) is the longest axis:
841 if P(1)<P(2) then
842     P([1 2]) = P([2 1]);
843     P(5) = 90 - P(5);
844 end
845
846 // DISPLAY
847 // -------
848 disp([Pa;P], dmin);
849 // A successful fit should lead to dmin ~ noise amplitude
850
851 clf
852 plot(xe, ye, "+")   // Data
853 // Model
854 Theta = 0:2:360;
855 x = P(1) * cosd(Theta);
856 y = P(2) * sind(Theta);
857 A = P(5);
858 [x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
859 plot(x, y, "r")
860 isoview
861 // Parameters values:
862 Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
863 Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
864 xstring(-3.7, 1.2, ["" "    a           b        xc       yc     A°"
865                     "Actual:" Patxt
866                     "Fit:" Ptxt ])
867 xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
868          ]]></programlisting>
869                 </td>
870                 <td>
871                     <scilab:image><![CDATA[
872 function g = Gap(p, Data)
873     // p: a, b, center xc, yc, alpha tilt °
874     C = cosd(p(5)); S = -sind(p(5));
875     x = Data(1,:) - p(3);
876     y = Data(2,:) - p(4);
877     g = ((x*C - y*S )/p(1)).^2 + ((x*S + y*C)/p(2)).^2 - 1;
878 endfunction
879
880 // Generating the data
881 // -------------------
882 // Actual parameters :
883 Pa = [3, 1.3, -2, 1.5, 30];
884 Np = 30;            // Number of data points
885 Theta = grand(1,Np, "unf",-100, 170);
886 // untilted centered noisy arc
887 x = Pa(1)*(cosd(Theta) + grand(Theta, "unf",-0.07, 0.07));
888 y = Pa(2)*(sind(Theta) + grand(Theta, "unf",-0.07, 0.07));
889 // Tilting and off-centering the arc
890 A = Pa(5);
891 C = cosd(A); S = sind(A);
892 xe =  C*x + y*S + Pa(3);
893 ye = -S*x + y*C + Pa(4);
894 Data = [xe ; ye];
895
896 // Retrieving parameters from noisy data
897 // -------------------------------------
898 // Initial estimation
899 ab0 = (strange(xe)+strange(ye))/2;
900 xc0 = mean(xe);
901 yc0 = mean(ye);
902 A0 = -atand(reglin(xe,ye));
903 P0 = [ab0 ab0/2 xc0 yc0 A0];
904 // Parameters bounds
905 Pinf = [ 0     0     xc0-2*ab0, yc0-2*ab0 -90];
906 Psup = [3*ab0 3*ab0  xc0+2*ab0, yc0+2*ab0  90];// Fitting
907
908 // FITTING
909 [P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
910
911 // P(1) is the longest axis:
912 if P(1)<P(2) then
913     P([1 2]) = P([2 1]);
914     P(5) = 90 - P(5);
915 end
916
917 // DISPLAY
918 // -------
919 clf
920 plot(xe, ye, "+")   // Data
921 // Model
922 Theta = 0:2:360;
923 x = P(1) * cosd(Theta);
924 y = P(2) * sind(Theta);
925 A = P(5);
926 [x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
927 plot(x, y, "r")
928 isoview
929 // Parameters values:
930 Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
931 Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
932 xstring(-3.7, 1.2, ["" "    a           b        xc       yc     A°"
933                     "Actual:" Patxt
934                     "Fit:" Ptxt ]);
935 xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
936 title("datafit() example #4", "fontsize", 3)
937 // ===============
938 gcf().axes_size = [500 500];
939 gcf().children(1).margins(1) = 0.05;
940 ]]>             </scilab:image>
941                 </td>
942             </tr>
943         </table>
944     </refsection>
945     <refsection role="see also">
946         <title>See also</title>
947         <simplelist type="inline">
948             <member>
949                 <link linkend="lsqrsolve">lsqrsolve</link>
950             </member>
951             <member>
952                 <link linkend="optim">optim</link>
953             </member>
954             <member>
955                 <link linkend="leastsq">leastsq</link>
956             </member>
957         </simplelist>
958     </refsection>
959     <refsection role="history">
960         <title>History</title>
961         <revhistory>
962             <revision>
963                 <revnumber>6.1</revnumber>
964                 <revdescription>
965                     <para>
966                         New option <varname>Wd</varname> to weight measured data <varname>Data</varname>.
967                     </para>
968                     <para>
969                         The gap function G() may now be vectorized over <literal>Data</literal> points.
970                     </para>
971                     <para>
972                         Initial parameters <varname>p0</varname> and their bounds may now be
973                         provided as a matrix instead of a mandatory column.
974                     </para>
975                     <para>
976                         The <literal>err</literal> output argument is replaced with
977                         <literal>dmin</literal> = average Model-to-data distance.
978                     </para>
979                     <para>
980                         New output argument <literal>status</literal> added (for the "qn" algorithm).
981                     </para>
982                 </revdescription>
983             </revision>
984         </revhistory>
985     </refsection>
986 </refentry>