powered by

DATA FITTING IN SCILAB
In this tutorial the reader can learn about data fitting, interpolation andapproximation in Scilab. Interpolation is very important in industrial applications
for data visualization and metamodeling.
Level

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License.

www.openeering.com

Step 1: Purpose of this tutorial
Many industrial applications require the computation of a fitting function inorder to construct a model of the data.
Two main data fitting categories are available:
Interpolation which is devoted to the development of numericalmethods with the constraint that the fitting function fits exactly all
the interpolation points (measured data);
Approximation which is devoted to the development of numericalmethods where the type of function is selected and then allparameters are obtained minimizing a certain error indicator toobtain the best possible approximation.

The first category is useful when data does not present noise, while thesecond one is used when data are affected by error and we want toremove error and smooth our model.

Step 2: Roadmap
In the first part we present some examples of polynomial interpolation andapproximation. After we propose exercises and remarks.

Data Fitting

www.openeering.com

Descriptions
Interpolation
Approximation or curve fitting
Notes
ExerciseConclusion and remarks

Steps
3-16
17-20
21-22
23
24-25

page 2/16

Step 3: Interpolation
The idea of approximation is to replace a function
with a function
selected from a given class of approximation function .

Two main cases exist:
Continuous function: In this case the function
is known
analytically and we want to replace it with an easier function, forexample we may want to replace a complex function with apolynomial for which integration or differentiation are easy;
Discrete function: In this case only some values of the function
are known, i.e.
and we want to make a mathematicalmodel
which is close to the unknown function
such that
it is possible to establish the value ofoutside the known
points.

Example of approximation of a continuous function using a piecewise
approximation

Step 4: Main class of interpolation function
Several families

of interpolation functions exist. The most common are:Polynomial interpolation of degree: In this case, weapproximate data with a polynomialof degree of the form

Piecewise polynomial: In this case the interval is subdivided intosubintervals in which we define a polynomial approximation of low
degree with or without continuity on the derivatives between eachsubinterval.

Data Fitting

Example piecewise interpolation with first derivative continuous anddiscontinuous connections

www.openeering.com

page 3/16

Step 5: Test case: Runge function
The Runge function is defined as

Typically is considered in the interval [-5,5].
In our examples we consider 7 interpolation points, denoted withequally distributed in the interval [-5,5].

,The following code implements the Runge function and perform
visualization.

// Define Runge function
deff('[y]=f(x)','y = 1 ./(1+x.^2)');// Interpolation points
xi = linspace(-5,5,7)'; yi = f(xi);// Dataxrunge = linspace(-5,5,101)'; yrunge = f(xrunge);

Interpolation function// Plot Runge function
scf(1); clf(1);plot(xrunge,yrunge,'b-');plot(xi,yi,'or');xlabel("x");ylabel("y");title("Runge function");

Data Fitting

www.openeering.com

page 4/16

Step 6: Piecewise constant interpolation
Piecewise constant interpolation is the simplest way to interpolate data. It
consists on locating the nearest data value and assigning the same value
to the unknown point.

Piecewise constant interpolation

Step 7: Piecewise constant interpolation in Scilab
The Scilab command used to perform piecewise interpolation is"interp1" where the third argument is "nearest". The fourth argument
specifies if an extrapolation method should be used when the evaluation
points are outside the interval of the interpolation points.
// Interpolation// Evaluation points
xval = linspace(-6,6,101)';
xx_c = xval;
yy_c = interp1(xi,yi,xx_c,'nearest','extrap');// Plotscf(3);
clf(3);
plot(xrunge,yrunge,'k-');plot(xx_c,yy_c,'b-');plot(xi,yi,'or');xlabel("x");ylabel("y");title("Piecewise interpolation");legend(["Runge func";"Interp.";"Interp. val"]);Data Fitting

Piecewise constant interpolationwww.openeering.com

page 5/16

Step 8: Piecewise linear interpolation
Linear interpolation is a polynomial of degree 1 that connects two points,, andthe interpolant is given by

Piecewise linear interpolation (green) and extrapolation (red)Step 9: Linear interpolation in Scilab
The Scilab command used to perform linear interpolation is again"interp1" but now the third argument is "linear". We can note that thecode is similar to the previous one for piecewise interpolation.
// Interpolation
xx_l = xval;
yy_l = interp1(xi,yi,xx_c,'linear','extrap');// Plotscf(4);
clf(4);
plot(xrunge,yrunge,'k-');plot(xx_l,yy_l,'b-');plot(xi,yi,'or');xlabel("x");ylabel("y");title("Linear interpolation");legend(["Runge func";"Interp.";"Interp. val"]);Data Fitting

Piecewise linear interpolationwww.openeering.com

page 6/16

Step 10: Polynomial interpolation
Given a set ofdata points, where all
are different, we findthe polynomial of degree which exactly passes through these points.
Polynomial interpolations may exhibits oscillatory effects at the end of thepoints. This is known as Runge phenomenon. For example, the Runge
function has this phenomenon in the interval [-5,+5] while in the interval [1,1] this effect is not present.
Polynomial interpolation

Step 11: Polynomial interpolation in Scilab
The Scilab command used to perform polynomial interpolation is"polyfit" which is included in the zip file (see references). The syntax
requires the interpolation points
and the degree of the polynomialinterpolation which is equal to the number of point minus one.
Note that the command "horner" evaluates a polynomial in a given set ofdata.
// Import function
exec("polyfit.sci",-1);// Interpolation
xx_p = xval;[Pn] = polyfit(xi, yi, length(xi)-1);
yy_p = horner(Pn,xx_p);// Plotscf(5);
clf(5);
plot(xrunge,yrunge,'k-');plot(xx_p,yy_p,'b-');plot(xi,yi,'or');xlabel("x");ylabel("y");title("Polynomial interpolation");legend(["Runge func";"Interp.";"Interp. val"]);Data Fitting

Polynomial interpolationwww.openeering.com

page 7/16

Step 12: Cubic spline interpolation
Cubic spline interpolation uses cubic polynomials on each interval. Theneach polynomial is connected to the next imposing further continuity
equations for the first and second derivatives.
Several kinds of splines are available and depend on how degrees offreedom are treated. The most well-known splines are: natural, periodic,not-a-knot and clamped.
Spline interpolation

Step 13: Cubic spline in Scilab
The Scilab command used to perform cubic spline interpolation is"splin". Several types of splines exist and can be specified by setting thethird argument of the function. The evaluation of the spline is done using
the command "interp" where it is possible to specify the extrapolation
strategy.
// Splines examples
d = splin(xi, yi,"not_a_knot");// d = splin(xi, yi,"natural");// d = splin(xi, yi,"periodic");
xx_s = xval;
yy_s = interp(xx_s, xi, yi, d, "linear");// Plotscf(6);
clf(6);
plot(xrunge,yrunge,'k-');plot(xx_s,yy_s,'b-');plot(xi,yi,'or');xlabel("x");ylabel("y");title("Spline interpolation");legend(["Runge func";"Interp.";"Interp. val"]);Data Fitting

Cubic spline interpolation with linear extrapolationwww.openeering.com

page 8/16

Step 14: Radial Basis Interpolation (RBF)
Radial basis function modeling consists of writing the interpolation function
as a linear combination of basis functions
that depends only on thedistanceof the interpolation points . This is equal to:

Many radial basis functions exist, the most famous are:Gaussian:
Multiquadratic:
Inverse multiquadratic:

Using the interpolation conditions
following linear system

we have to solve theThin plate spline:

where the element.

Step 15: Gaussian RBF in Scilab
In our example we use the Gaussian RBF.
// Gaussian RBF
deff('[y]=rbf_gauss(r,sigma)','y = exp(-r.^2 ./(2*sigma))');// Plotscf(7);
clf(7);
r = linspace(0,3);
y1 = rbf_gauss(r,0.1);
y2 = rbf_gauss(r,1.0);
y3 = rbf_gauss(r,2.0);
plot(r,y1,'k-');plot(r,y2,'b-');plot(r,y3,'r-');xlabel("$r$");ylabel("$\phi(r)$");
title("Gaussian rbf");legend(["$\sigma = 0.1$";"$\sigma = 1.0$";"$\sigma = 2.0$"]);

Data Fitting

www.openeering.com

Gaussian RBF for different value of sigma

page 9/16

Step 16: RBF in Scilab
On the right we report the optimal Gaussian RBF obtained for theinterpolation of the Runge function where we have optimized theparameters. The full code is reported in the Openeering web site.

Generally, each radial basis function depends on a parameter and thisparameter is known as modeling parameter. This parameter has effect on
the oscillation behavior of the function and the optimal choice is not aneasy task. Many techniques are available for finding the best modelingparameter, the most famous is probably the “leave one out”.

Behaviour of the error versus sigma

The leave-one-out cross-validation (LOOCV) consists in using a single
point from the original set of data as a validation data. The validation of
the model is given by that point. This process can be repeated for all thepoints in the data set such that, in the end, all the points are used once asvalidation point. This method can produce a mean value of all these leaveone-out errors and gives a global estimate of the model.
Since a value that estimates the model is available, we can use anoptimization solver for finding the best parameter in order to minimize theerror of the model.

Optimal RBF with

Data Fitting

www.openeering.com(left), non optimal solution(right)page 10/16

Step 17: Approximation or curve fitting
When data is affected by errors, polynomial interpolation cannot beappropriate since the approximation function is constrained to be throughthe interpolation points. So it makes sense to fit the data starting from agiven class of functions and minimizing the difference between the data
and the class of functions, i.e.

The "polyfit" function computes the best least square polynomial
approximation of data. If we choosein the "polyfit" function, we
approximate data with linear function of the form
, i.e. wecompute the linear least squares fitting.

Schematic example of min interpretation

Step 18: 1D approximation
In this example we add noise to the function
make polynomial approximation of order 1 and 2.

and then we

The critical code is reported here (only case 1):
np = 100; noise = 0.7*(rand(np,1)-0.5);
x = linspace(0,2,np)';
yexact = x.^2 + x;
ynoise = yexact + noise;// degree 1 approximation
p1 = polyfit(x, ynoise, 1);
p1val = horner(p1,x);
scf(10); clf(10);plot(x,yexact,'k-'); plot(x,ynoise,'b-'); plot(x,p1val,'r-');Comparison of best fitting for degree 1 and 2

For details download the zip file with the source codes.

Data Fitting

www.openeering.com

page 11/16

Step 19: 2D approximation
In this example we want to approximate scattered data with a linear leastsquare fitting in two dimensions.
Given the polynomial approximationand thescattered interpolation points, the optimal computation of theunknown parametersrequires to solve the following
overdetermined linear system

The matrix is known as the Vandermonde matrix and arises in polynomialinterpolation. The i-th row of the matrix corresponds to the polynomial
evaluated at point i-th. In our case, the i-th row corresponds to:This is done using the Singular Value Decomposition (SVD) or
equivalently using the backslash command. Here, we report only thesolution stage. The full code can be downloaded from our web page.
// Generating random points along a plane
np = 30;
noise = 0.5*(rand(np,1)-0.5);

Example of least square approximation in 2D// Extract datax = rand(np,1);
y = rand(np,1);
znoise = -x+2*y+noise;// Vandermonde matrix for P(x,y) = a+b*x+c*y
V = [ones(np,1),x,y];// Find coefficient i.e. minimize error norm
coeff = V\znoise;

Data Fitting

www.openeering.com

page 12/16

Step 20: nD linear approximation
The previous example can be easily extended to an n-dimensional
problem, i.e. we search for an approximation of the formn = 10;// Problem dimensions
neval = 100;// Number of evaluation points
pcoeff = 1:n+1; // Problem coefficient ([a0, a1, ..., an])pcoeff([2,5,8]) = 0; // Some zero coefficients// Evaluation points[xdata,ydata] = generatedata(pcoeff, n, neval);.
Here, we report the code for the problem where some coefficients may beequal to zero and could be deleted from the estimation problem.

The idea here implemented consists of the following strategy:// Estimate coefficients
pstar = estimatecoeff(xdata, ydata);// Find "zero" coefficient (define tolerance)
tol = 0.1;
zeroindex = find(abs(pstar)<=tol);// re-estimate coefficientspzero = estimatecoeffzero(xdata, ydata, zeroindex);

1. Perform least square approximation of data;

The code use the following functions:2. Find the zero coefficients with some pre-fixed tolerance;
3. Re-perform least square approximation of the reduce problem
where we have deleted the columns of the Vandermonde matrix
corresponding to the zero coefficients.

The code reported on the right estimates these coefficients. The full code
can be downloaded from our website.

y=evaldata(pcoeff,x)that evaluates the polynomialon point
noise on the definition of the coefficients;

defined by coefficients
adding a uniform[xdata,ydata]=generatedata(pcoeff,n,neval)that generates the evaluation points and their values;
pstar=estimatecoeff(xdata,ydata)that estimates the polynomial coefficients using the least square
method;pzero=estimatecoeffzero(xdata,ydata,zeroindex)
that estimates the polynomial coefficients where some coefficients,
specified by the vector zeroindex , are equal to zero.

Data Fitting

www.openeering.com

page 13/16

Step 21: Another polyfit function
The numerical solution of the polynomial interpolation problem requires tosolve the linear system with the Vandermonde matrix. Numerically, thisproblem is ill-conditioned and requires an efficient strategy for a “correct”
solution. This is performed by using, for example, QR factorization.
The function “polyfit_fulldemo.sce”, provided by Konrad Kmieciak,
and contained in the zip file, is an alternative to the polyfit function.

function [u]=polyfit(x, y, n)// Vandermondefor k=n:-1:0
if k==n then w=0;
end
w=w+1;
Xu(:,w)=[x.^k];
end
// QR
[q r k]=qr(Xu,'0');s = inv(r) * (q' * y); // s = r \ (q' * y)for o=1:length(s)u(find(k(:,o)>0))=s(o);endendfunction

Step 22: Industrial applications of data fitting
Interpolation and approximation are two major techniques for constructing
mathematical models. Mathematical models can substitute complex modelin real-case applications.
When the original model is complex, or when it requires long and costly
evaluations (for example with finite element analysis – FEA), a simplified
model of the original model is required.
The model of the model is often called metamodel (a.k.a. responsesurfaces) and the metamodeling technique is widely used in industrial
applications.

Metamodeling+Response surface+
Self-Organizing Maps+
Neural Networks=Industrial applications

In real-case applications, when optimizing product or services, we need toevaluate several times the model. This means that having a metamodelthat can be evaluated faster is, most of the time, the only way for findingoptimal solutions.

Data Fitting

www.openeering.com

page 14/16

Step 23: Exercise on exponential decay
As an exercise, try to estimate the coefficient of a decay function of theform:

using a linear least square fitting of data.

Hits: Use a logarithmic change of variable to reduce the problem in thefollowing form

Exponential decay fitting of data

Step 24: Concluding remarks and References
1. Scilab Web Page: Available: www.scilab.org.
In this Scilab tutorial we have shown how to apply data fitting in Scilab
starting from piece-wise interpolation, presenting polynomial fitting andcubic spline, and ending up with radial basis functions (RBF).

2. Openeering: www.openeering.com.
3. Javier I. Carrero is the author of the polyfit function. The originalcode is available on the Scilab.org web pages and is included in thezip file.

In the case of polynomial interpolation we used the function "polyfit"
which can be downloaded from the Scilab webpage and is included in theprovided source codes.

Data Fitting

www.openeering.com

page 15/16

Step 25: Software content
To report bugs or suggest improvements please contact the Openeering
team. We thank Konrad Kmieciak for reporting us a new version of thepolyfit function.
-------------Main directory-------------ex0.sce
ex1.sceinterpolation.sce
polyfit.sci
polyfit_manual.pdf
polyfit_fulldemo.sce
estimate_lincoeff.sce
license.txt

www.openeering.com: Plotting of the first figure: Solution of exercise 1: Interpolation examples codes: Polyfit function: Polyfit manual: Another version of polyfit: Estimantion of nD linear model: The license fileThank you for your attention,
Manolo Venturin
Silvia Poles

Data Fitting

www.openeering.com

page 16/16