Curve Fitting with Scilab

Neither Scilab nor Scicoslab have a function for straight curve fitting, such as the polyfit function that we can find in Matlab. However, it’s not that difficult to develop (or find?) a custom made function for that purpose. 

In [1] we can find a suggestion for this task. This is an edited version of the transcription of the 15-line code.
 

function p = polyfit(x, y, n)
  if length(x) ~= length(y)
    error('x and y vectors must be the same size')
  end 

  x = x( : );
  y = y( : );
  V = ones(length(x), n+1); 

  for j = n : -1 : 1
    V(:, j) = x .* V (:, j+1);
  end 

  [Q, R] = qr(V);
  QTy = Q' * y;
  p = R(1 : n+1, 1 : n+1) \ QTy(1 : n+1);
  p = p.';
endfunction
 

There are three inputs: vectors x and y, and a scalar n. x and y are the vectors defining our points to be fitted, and n is the order of the polynomial to be found.
 

On the other hand, we also need an equivalent of the polyval function given in Matlab. Naturally there are many possibilities here. The following code is much simpler than the previous one and we’re showing our suggestion.
 

function y = polyval(p, x)
  y = 0*x;
  p = mtlb_fliplr(p);
  for ix = 1 : length(p)
    y =  y + p(ix) * x.^(ix-1);
  end
endfunction
 

The goal of this function is to use the polynomial found by polyfit, receive an x-value and release the corresponding y-value produced by that polynomial. There are two inputs: we need the coefficients of the polynomial, in p, and the x-value, in x. We just do the required math in a simple and straight way. We use the function mtlb_fliplr (provided in the M2SCI toolbox) just to sort the order of the exponent of the variable. 

Now, let’s see in action both functions.
 

// Prepare our environment
clear; clc 

// Declare our two functions
getf('polyfit.sci')
getf('polyval.sci') 

// Define (arbitrarily) the number of points to take into account
np = 10; 

// Define the x-vector and two functions, the second function
// is a noised version of the first one
x = linspace(0, 1, np);
y1 = x.^3 - 5*x.^2 - 3*x - 7;
y2 = y1 + .1*rand(1, np); 

// Enter the x and y vectors, and the order of the polynomial
// that we want to obtain
p = polyfit(x, y2, 3) 

// Define other x-values and find the original function
x = 0 : .4 : 1
y1 = x.^3 - 5*x.^2 - 3*x - 7 

// Use polyval to find the equivalent values in the
// noised function
y = polyval(p, x) 

// Divide the values just for comparison purposes
ratio = y1./y
 

Now, the (commented) results shown by Scilab are: 

// This is our 3rd. order polynomial from the noised function
p  =  1.1103879  - 5.0400213  - 3.052195  - 6.9291697 

// Another set of points generated
x  =  0.    0.4    0.8
y1 =  - 7.  - 8.936  - 12.088 

// Find the y-values of the noised function
y  =  - 6.9291697  - 8.8853863  - 12.028021 

// Compare the original vs the noised results
ratio = 1.0102221    1.0056963    1.0049866 
 

We can see that the polyfit and polyval functions coded in Scilab are working pretty fine.  

Reference:
[1] Medrano Sanchez, C., Plaza Garcia, I.; Software libre para Calculo Numerico; Alfaomega; Mexico, DF, 2009.