</refnamediv>
<refsynopsisdiv>
<title>Syntax</title>
- <synopsis>[p, fmin, status] = datafit([iprint,] G [,DG], Data [,Wd] [,Wg][,'b',p_inf,p_sup], p0 [,algo][,stop])
+ <synopsis>
+ [p, fmin, status] = datafit(G, Data, p0)
+ [p, fmin, status] = datafit(G, Data, Wd, p0)
+ [p, fmin, status] = datafit(G, Data, Wd, 'b', p_inf, p_sup, p0)
+ [p, fmin, status] = datafit(G, Data, 'b', p_inf, p_sup, p0)
+ [p, fmin, status] = datafit(G, Data [,Wd] [,Wg][,'b',p_inf,p_sup], p0 [,algo][,stop])
+ [p, fmin, status] = datafit(G, DG, Data, ..)
+ [p, fmin, status] = datafit(iprint, ..)
</synopsis>
</refsynopsisdiv>
<refsection>
</refsection>
<refsection>
<title>Examples</title>
- <table>
- <tr>
- <td>
- <para>
- <emphasis role="bold">Simple example: Polynomial fitting (parabola = 3 parameters)</emphasis>.
- </para>
- <programlisting role="example"><![CDATA[
+ <para>
+ <emphasis role="bold">Simple example: Polynomial fitting (parabola = 3 parameters)</emphasis>.
+ </para>
+ <programlisting role="example"><![CDATA[
// Model (parabola)
function y = FF(x,p)
y = (p(1)*x + p(2)).*x + p(3)
pa = [2;-4;14] // parameters of the actual law, to generate data
X = 0:0.1:3.5;
Y0 = FF(X, pa);
-noise = 10*(rand(Y0)-.5); // We add some noise
+noise = 5*(rand(Y0)-.5); // We add some noise
Y = Y0 + noise
Data = [X ; Y];
tmp = ax.x_ticks.labels;
gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
xgrid(color("grey70"))
- ]]></programlisting>
- </td>
- <td>
- <scilab:image>
- // Model (parabola)
- function y = FF(x,p)
- y = p(1)*x.*x + p(2)*x + p(3)
- endfunction
-
- // Data creation
- pa = [2;-4;14] // parameters of the actual law, used to generate data
- X = 0:0.1:3.5;
- Y0 = FF(X, pa);
- noise = 10*(rand(Y0)-.5); // We add some noise
- Y = Y0 + noise
- Data = [X ; Y];
-
- // Gap function definition
- function e = G(p, Data)
- e = Data(2,:) - FF(Data(1,:), p); // vertical distance
- endfunction
-
- // Performing the fitting
- // ----------------------
- // Without weighting data:
- p0 = [1;-1;10]
- [p, dmin] = datafit(G, Data, p0);
-
- // The uncertainty is identified to the noise value.
- // The weight is set inversely proportional to the uncertainty
- [pw, dmin] = datafit(G, Data, 1 ./abs(noise), p0);
-
- // Display
- // -------
- scf(0);
- clf()
- xsetech([0 0 1 0.83])
- plot2d(X, Y0, 24) // True underlaying law
- plot2d(X, Y, -1) // Noisy data
- plot2d(X, FF(X,p), 15) // best unweighted fit
- plot2d(X, FF(X,pw), 18) // best weighted fit
- legend(["True law" "Noisy data" "Unweighted fit" "Weighted fit"],2)
-
- xsetech([0 .75 1 0.25])
- plot2d(X, FF(X,p)-Y0, 15) // remaining gap with best unweighted fit
- plot2d(X, FF(X,pw)-Y0, 18) // remaining gap best weighted fit
- ylabel("Data fit - True law")
- ax = gca(); ax.x_location = "top";
- tmp = ax.x_ticks.labels;
- gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
- xgrid(color("grey70"))
- // ===============
- gcf().axes_size = [500 500];
- gcf().children(1).margins(1) = 0.07;
- gcf().children(2).margins(1) = 0.07;
- </scilab:image>
- </td>
- </tr>
- <tr>
- <td>
- <para/>
- <para>
- <emphasis role="bold">Fitting a gaussian curve + sloping base (5 parameters)</emphasis>:
- </para>
- <programlisting role="example"><![CDATA[
+ ]]></programlisting>
+ <!-- Due to random input data, it is preferable to have a fixed image instead of automated.
+ <scilab:image>
+ // ===============
+ gcf().axes_size = [520 550];
+ </scilab:image>
+ -->
+ <inlinemediaobject>
+ <imageobject>
+ <imagedata fileref="../../images/datafit_A.png"/>
+ </imageobject>
+ </inlinemediaobject>
+ <para/>
+ <para>
+ <emphasis role="bold">Fitting a gaussian curve + sloping base (5 parameters)</emphasis>:
+ </para>
+ <programlisting role="example"><![CDATA[
function Y = model(P, X)
Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
endfunction
// -------------------------------
// Gap function
-function r = gap(P,XY)
+function r = gap(P, XY)
r = XY(2,:) - model(P,XY(1,:))
endfunction
// -------------------------------
// Legends: Average vertical Model-Data distance:
// True one v/s according to the residual cost
L_1 = "Unweighted fit (%5.2f /%5.2f)";
-L_1 = msprintf(L1, mean(abs(fY-Y0)), sqrt(dmin/Np));
+L_1 = msprintf(L_1, mean(abs(fY-Y0)), sqrt(dmin/Np));
L_2 = "Weighted fit (%5.2f /%5.2f)";
-L_2 = msprintf(L1, mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)));
+L_2 = msprintf(L_2, mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)));
legend(["True law", "Noisy Data", L_1, L_2],1)
// Params: row#1: true params row#2:
ylabel("Data fit - True law")
ax = gca();
ax.x_location = "top";
-tmp = ax.x_ticks.labels;
-gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+gca().x_ticks.labels = emptystr(size(ax.x_ticks.labels, "*"),1);
xgrid(color("grey70"))
- ]]></programlisting>
- </td>
- <td>
- <scilab:image>
- function Y = model(P,X)
- Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
- endfunction
- // -------------------------------
- // Gap function
- function r = gap(P,XY)
- r = XY(2,:) - model(P,XY(1,:))
- endfunction
- // -------------------------------
- // True parameters
- Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
-
- // Generating data
- // ---------------
- X = 0:0.2:10;
- Np = length(X);
- Y0 = model(Pg, X); // True law
- noise = (rand(Y0)-.5)*Pg(3)/4;
- Y = Y0 + noise // Noised data
- Data = [X ; Y];
-
- // Let's try to find the original parameters from the (noisy) data:
- // ---------------------------------------------------------------
- // Performing the fitting = parameters optimization
- [v, k] = max(Y0);
- P0 = [X(k) 1 v 1 1];
- [P, dmin, s] = datafit(gap, Data, P0);
- Wd = 1./abs(noise);
- [Pw, dminw, s] = datafit(gap, Data, Wd, P0);
-
- // Computing fitting curves from found parameters
- fY = model(P, X);
- fYw = model(Pw,X);
-
- // Display
- // -------
- scf(0);
- clf()
- xsetech([0 0 1 0.82])
- plot2d(X, Y0, 24) // True underlaying law
- plot2d(X, Y, -1) // Noisy data
- plot2d(X, fY, 15) // best unweighted fit
- plot2d(X, fYw,18) // best weighted fit
- legend(["True law" "Noisy Data" ..
- msprintf("Unweighted fit (%5.2f /%5.2f)",mean(abs(fY-Y0)), sqrt(dmin/Np)) ..
- msprintf("Weighted fit (%5.2f /%5.2f)", mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)))],1)
- // Average vertical Model-Data distance: True one / according to the residual cost
- // Params: row#1: true params row#2: from unweighted fit row#3: from weighted fit
- [xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
- xnumb(xs, ys, [Pg ; P ; Pw])
- xstring(2,1.3,"$\Large p_3\,.\,exp\left({-(x-p_1)^2}\over{2\,p_2} \right)+p_4.x+p_5$")
-
- // Plotting residus:
- xsetech([0 .75 1 0.25])
- plot2d(X, model(P ,X)-Y0, 15) // remaining gap with best unweighted fit
- plot2d(X, model(Pw,X)-Y0, 18) // remaining gap best weighted fit
- ylabel("Data fit - True law")
- ax = gca();
- ax.x_location = "top";
- tmp = ax.x_ticks.labels;
- gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
- xgrid(color("grey70"))
- // ===============
- gcf().axes_size = [500 600];
- </scilab:image>
- </td>
- </tr>
- <tr>
- <td>
- <para/>
- <para>
- <emphasis role="bold">Extraction of 3 overlapping gaussian curves (9 parameters)</emphasis>:
- example with constrained bounds and a matrix of parameters.
- </para>
- <programlisting role="example"><![CDATA[
+ ]]></programlisting>
+ <!-- Due to random input data, it is preferable to have a fixed image instead of automated.
+ <scilab:image>
+ // ===============
+ gcf().axes_size = [520 550];
+ </scilab:image>
+ -->
+ <inlinemediaobject>
+ <imageobject>
+ <imagedata fileref="../../images/datafit_B.png"/>
+ </imageobject>
+ </inlinemediaobject>
+ <para/>
+ <para>
+ <emphasis role="bold">Extraction of 3 overlapping gaussian curves (9 parameters)</emphasis>:
+ example with constrained bounds and a matrix of parameters.
+ </para>
+ <programlisting role="example"><![CDATA[
// Basic gaussian component:
function Y = gauss(X, average, stDev, ymax)
Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
// Gap function:
function r = gap(P,XY)
- r = XY(2,:) - model(P,XY(1,:))
+ r = XY(2,:) - model(P, XY(1,:))
endfunction
-// Actual parameters :
+// True parameters :
Pa = [ 1.1 4 7 // averages
1 0.8 1.5 // standard deviations
3 2 5.5 // ymax
- ];
+ ];
Nc = size(Pa,2); // Number of summed curves
Np = 200; // Number of data points
xmax = 10;
X = linspace(xmin-2, xmax+2, Np);
Y = model(Pa, X);
-Data = [X ; Y]; // Curve to analyze / fit
+Data = [X ; Y]; // Curve to analyze / fit
// FITTING
// -------
-// Fitting parameters:
-// Initial values & bounds (amplitudes must be > 0)
-Po = [ linspace(xmin,xmax,Nc)
- ones(1,Nc)*0.5
- ones(1,Nc)*max(Y)/2 ];
+// Fitting parameters: initial values and bounds (amplitudes must be > 0)
+Po = [ linspace(xmin,xmax,Nc) ; ones(1,Nc)*0.5 ; ones(1,Nc)*max(Y)/2];
Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3 ; ones(1,Nc)*max(Y)*1.2];
// Fitting
-[P, dmin] = datafit(gap,Data,'b',Pinf,Psup,Po,'qn','ar',200,200);
+[P, dmin] = datafit(gap, Data, 'b',Pinf,Psup, Po, 'qn','ar',200,200);
// Display
// -------
plot(X, Y, "-b", X, model(P,X), "-r");
gca().children.children(1:2).thickness = 2;
for i = 1:Nc
- // True gaussian components
- plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c");
- // Extracted gaussian components
- plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g");
+plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c"); // True gaussian components
+plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g"); // Extracted gaussian components
end
-legend("Data", "Parametric Model", ..
- "True components", "Extracted components", 2);
+legend("Data", "Parametric Model","True components", "Extracted components", 2);
xsetech([0 0.75 1 0.25])
plot(X, model(P,X)-Y);
gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
ylabel("Model - Data")
xgrid(color("grey70"))
- ]]></programlisting>
- </td>
- <td>
- <scilab:image>
- // Basic gaussian component:
- function Y = gauss(X, average, stDev, ymax)
- Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
- endfunction
-
- // Model: cumulated gaussian laws
- function r = model(P, X)
- // P(1,:): averages P(2,:): standard deviation P(3,:): ymax
- r = 0;
- for i = 1:size(P,2)
- r = r + gauss(X, P(1,i), P(2,i), P(3,i));
- end
- endfunction
-
- // Gap function:
- function r = gap(P,XY)
- r = XY(2,:) - model(P,XY(1,:))
- endfunction
-
- // True parameters :
- Pa = [ 1.1 4 7 // averages
- 1 0.8 1.5 // standard deviations
- 3 2 5.5 // ymax
- ];
- Nc = size(Pa,2); // Number of summed curves
- Np = 200; // Number of data points
-
- // Data generation (without noise)
- xmin = 0;
- xmax = 10;
- X = linspace(xmin-2, xmax+2, Np);
- Y = model(Pa, X);
- Data = [X ; Y]; // Curve to analyze / fit
-
- // FITTING
- // -------
- // Fitting parameters: initial values and bounds (amplitudes must be > 0)
- Po = [ linspace(xmin,xmax,Nc) ; ones(1,Nc)*0.5 ; ones(1,Nc)*max(Y)/2];
- Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
- Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3 ; ones(1,Nc)*max(Y)*1.2];
- // Fitting
- [P, dmin] = datafit(gap, Data, 'b',Pinf,Psup, Po, 'qn','ar',200,200);
-
- // Display
- // -------
- clf()
- xsetech([0 0 1 0.8]);
- plot(X, Y, "-b", X, model(P,X), "-r");
- gca().children.children(1:2).thickness = 2;
- for i = 1:Nc
- plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c"); // True gaussian components
- plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g"); // Extracted gaussian components
- end
- legend("Data", "Parametric Model","True components", "Extracted components", 2);
-
- xsetech([0 0.75 1 0.25])
- plot(X, model(P,X)-Y);
- gca().x_location = "top";
- gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
- ylabel("Model - Data")
- xgrid(color("grey70"))
- // ===============
- gcf().axes_size = [500 500];
- </scilab:image>
- </td>
- </tr>
- <tr>
- <td>
- <para>
- <emphasis role="bold">With a parametric curve: Fitting a off-centered tilted 2D elliptical arc (5 parameters)</emphasis>.
- </para>
- <programlisting role="example"><![CDATA[
+ ]]></programlisting>
+ <scilab:image>
+ // Basic gaussian component:
+ function Y = gauss(X, average, stDev, ymax)
+ Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
+ endfunction
+
+ // Model: cumulated gaussian laws
+ function r = model(P, X)
+ // P(1,:): averages P(2,:): standard deviation P(3,:): ymax
+ r = 0;
+ for i = 1:size(P,2)
+ r = r + gauss(X, P(1,i), P(2,i), P(3,i));
+ end
+ endfunction
+
+ // Gap function:
+ function r = gap(P,XY)
+ r = XY(2,:) - model(P,XY(1,:))
+ endfunction
+
+ // True parameters :
+ Pa = [ 1.1 4 7 // averages
+ 1 0.8 1.5 // standard deviations
+ 3 2 5.5 // ymax
+ ];
+ Nc = size(Pa,2); // Number of summed curves
+ Np = 200; // Number of data points
+
+ // Data generation (without noise)
+ xmin = 0;
+ xmax = 10;
+ X = linspace(xmin-2, xmax+2, Np);
+ Y = model(Pa, X);
+ Data = [X ; Y]; // Curve to analyze / fit
+
+ // FITTING
+ // -------
+ // Fitting parameters: initial values and bounds (amplitudes must be > 0)
+ Po = [ linspace(xmin,xmax,Nc) ; ones(1,Nc)*0.5 ; ones(1,Nc)*max(Y)/2];
+ Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
+ Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3 ; ones(1,Nc)*max(Y)*1.2];
+ // Fitting
+ [P, dmin] = datafit(gap, Data, 'b',Pinf,Psup, Po, 'qn','ar',200,200);
+
+ // Display
+ // -------
+ clf()
+ xsetech([0 0 1 0.8]);
+ plot(X, Y, "-b", X, model(P,X), "-r");
+ gca().children.children(1:2).thickness = 2;
+ for i = 1:Nc
+ plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c"); // True gaussian components
+ plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g"); // Extracted gaussian components
+ end
+ legend("Data", "Parametric Model","True components", "Extracted components", 2);
+
+ xsetech([0 0.75 1 0.25])
+ plot(X, model(P,X)-Y);
+ gca().x_location = "top";
+ gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
+ ylabel("Model - Data")
+ xgrid(color("grey70"))
+ // ===============
+ gcf().axes_size = [520 400];
+ </scilab:image>
+ <para/>
+ <para>
+ <emphasis role="bold">With a parametric curve: Fitting a off-centered tilted 2D elliptical arc (5 parameters)</emphasis>.
+ </para>
+ <programlisting role="example"><![CDATA[
function g = Gap(p, Data)
// p: a, b, center xc, yc, alpha tilt °
C = cosd(p(5)); S = -sind(p(5));
"Actual:" Patxt
"Fit:" Ptxt ])
xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
- ]]></programlisting>
- </td>
- <td>
- <scilab:image><![CDATA[
-function g = Gap(p, Data)
- // p: a, b, center xc, yc, alpha tilt °
- C = cosd(p(5)); S = -sind(p(5));
- x = Data(1,:) - p(3);
- y = Data(2,:) - p(4);
- g = ((x*C - y*S )/p(1)).^2 + ((x*S + y*C)/p(2)).^2 - 1;
-endfunction
-
-// Generating the data
-// -------------------
-// Actual parameters :
-Pa = [3, 1.3, -2, 1.5, 30];
-Np = 30; // Number of data points
-Theta = grand(1,Np, "unf",-100, 170);
-// untilted centered noisy arc
-x = Pa(1)*(cosd(Theta) + grand(Theta, "unf",-0.07, 0.07));
-y = Pa(2)*(sind(Theta) + grand(Theta, "unf",-0.07, 0.07));
-// Tilting and off-centering the arc
-A = Pa(5);
-C = cosd(A); S = sind(A);
-xe = C*x + y*S + Pa(3);
-ye = -S*x + y*C + Pa(4);
-Data = [xe ; ye];
-
-// Retrieving parameters from noisy data
-// -------------------------------------
-// Initial estimation
-ab0 = (strange(xe)+strange(ye))/2;
-xc0 = mean(xe);
-yc0 = mean(ye);
-A0 = -atand(reglin(xe,ye));
-P0 = [ab0 ab0/2 xc0 yc0 A0];
-// Parameters bounds
-Pinf = [ 0 0 xc0-2*ab0, yc0-2*ab0 -90];
-Psup = [3*ab0 3*ab0 xc0+2*ab0, yc0+2*ab0 90];// Fitting
-
-// FITTING
-[P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
-
-// P(1) is the longest axis:
-if P(1)<P(2) then
- P([1 2]) = P([2 1]);
- P(5) = 90 - P(5);
-end
-
-// DISPLAY
-// -------
-clf
-plot(xe, ye, "+") // Data
-// Model
-Theta = 0:2:360;
-x = P(1) * cosd(Theta);
-y = P(2) * sind(Theta);
-A = P(5);
-[x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
-plot(x, y, "r")
-isoview
-// Parameters values:
-Patxt = msprintf("%5.2f %5.2f %5.2f %5.2f %5.1f", Pa);
-Ptxt = msprintf("%5.2f %5.2f %5.2f %5.2f %5.1f", P);
-xstring(-3.7, 1.2, ["" " a b xc yc A°"
- "Actual:" Patxt
- "Fit:" Ptxt ]);
-xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
-title("datafit() example #4", "fontsize", 3)
-// ===============
-gcf().axes_size = [500 500];
-gcf().children(1).margins(1) = 0.05;
-]]> </scilab:image>
- </td>
- </tr>
- </table>
+ ]]></programlisting>
+ <!-- Due to random input data, it is preferable to have a fixed image instead of automated.
+ <scilab:image>
+ // ===============
+ gcf().axes_size = [520 550];
+ </scilab:image>
+ -->
+ <inlinemediaobject>
+ <imageobject>
+ <imagedata fileref="../../images/datafit_C.png"/>
+ </imageobject>
+ </inlinemediaobject>
</refsection>
<refsection role="see also">
<title>See also</title>
<revision>
<revnumber>6.1</revnumber>
<revdescription>
- <para>
- New option <varname>Wd</varname> to weight measured data <varname>Data</varname>.
- </para>
- <para>
- The gap function G() may now be vectorized over <literal>Data</literal> points.
- </para>
- <para>
- Initial parameters <varname>p0</varname> and their bounds may now be
- provided as a matrix instead of a mandatory column.
- </para>
- <para>
- The <literal>err</literal> output argument is replaced with
- <literal>dmin</literal> = average Model-to-data distance.
- </para>
- <para>
- New output argument <literal>status</literal> added (for the "qn" algorithm).
- </para>
+ <itemizedlist>
+ <listitem>
+ New option <varname>Wd</varname> to weight measured data
+ <varname>Data</varname>.
+ </listitem>
+ <listitem>
+ The gap function G() may now be vectorized over <literal>Data</literal>
+ points.
+ </listitem>
+ <listitem>
+ Initial parameters <varname>p0</varname> and their bounds may now be
+ provided as a matrix instead of a mandatory column.
+ </listitem>
+ <listitem>
+ The <literal>err</literal> output argument is replaced with
+ <literal>dmin</literal> = average Model-to-data distance.
+ </listitem>
+ <listitem>
+ New output argument <literal>status</literal> added
+ (for the "qn" algorithm).
+ </listitem>
+ </itemizedlist>
</revdescription>
</revision>
</revhistory>
projects / scilab.git / commitdiff
