improve fit compared to R

I struggle to fit a curve through 8 dots in Matlab. in R, exact same values of x,y,and x_values_finer, and 3 lines of code:

fitt <- smooth.spline(x, y, spar = 0.5)

predicted <- predict(fitt, x_values_finer)

frst_drvt <- predict(fitt, x_values_finer, deriv = 1)$y

gives me a nice u-ish, v-ish shape curve, and i get smooth first derivative.

in Matlab, nothing comes close. the line below gives me a straight line, where even visually i can see that is wrong.

fitt = csaps(x, y, 0.5);

x values =

0.858734025700000

0.975248753121875

1.02525125312812

1.13140821289062

1.15927407430000

1.18768630564687

1.21665290240000

1.27628156250000

y values =

0.0255926963287944

0.0206749804955220

0.0216695267388959

0.0221177913199036

0.0248476224632830

0.0277388242602486

0.0305936964241147

0.0362101199563917

x_values_finer is a 1x5684 double that goes from 0.5904900 to 2.0112400, like this:

0.590490000000000 0.590740000000000 0.590990000000000 0.591240000000000

to be precise, x_values_finer in R is defined like this:

x_values_finer=seq((1+(-0.1))^5, (1+0.15)^5, (0.15-(-0.1))/1000)

can you shed some light why it is so difficult for Matlab, or can offer solutions?

this is the best i can get, even in the x values area it is not very smooth, but the worst is that it drops to zero bewyond bounds of original x values:

[Renate]

this nice thing i get in R

Answer 1

[John D'Errico]

x = [0.858734025700000
0.975248753121875
1.02525125312812
1.13140821289062
1.15927407430000
1.18768630564687
1.21665290240000
1.27628156250000];
y = [0.0255926963287944
0.0206749804955220
0.0216695267388959
0.0221177913199036
0.0248476224632830
0.0277388242602486
0.0305936964241147
0.0362101199563917];

fitt = csaps(x, y, 0.5)

fitt = struct with fields:
      form: 'pp'
    breaks: [0.8587 0.9752 1.0253 1.1314 1.1593 1.1877 1.2167 1.2763]
     coefs: [7x4 double]
    pieces: 7
     order: 4
       dim: 1

So fitt is a smoothing spline, in a pp-form.

fplot(@(X) ppval(fitt,X),[min(x),max(x)])
hold on
plot(x,y,'bo')

What is wrong? YOU USED THE WRONG VALUE FOR THE SMOOTHING PARAMETER! There is no presiumption that the smoothing parameter in R will be the same as in MATLAB, for two differenti codes, written by two different people.

If the smoothing parameter is too small, then the rsult will be a curve close to a straight line fit. (You would have learned this had you read the help for csaps.)

help csaps

 csaps - Cubic smoothing spline
    This MATLAB function returns the cubic smoothing spline interpolation to
    the given data (x,y) in ppform.

    Syntax
      pp = csaps(x,y)
      pp = csaps(x,y,p)
      pp = csaps(x,y,p,[],w)
      values = csaps(x,y,p,xx)
      values = csaps(x,y,p,xx,w)
      [___] = csaps({x1,...,xm},y,___)
      [___,P] = csaps(___)

    Input Arguments
      x - Data sites
        vector | cell array
      y - Data values to fit
        vector | matrix | array
      p - Smoothing parameter
        scalar in the range [0,1] | vector | cell array | empty array
      w - Error measure weights
        vector | cell array
      xx - Evaluation points
        vector | cell array

    Output Arguments
      pp - Spline in ppform
        spline structure
      values - Evaluated spline
        vector | matrix | array
      P - Smoothing parameter
        scalar | cell array

    Examples
      openExample('curvefit/FitSplinesWithDifferentSmoothingParametersExample')
      openExample('curvefit/AdjustSmoothingParametersAndWeightsExample')
      openExample('curvefit/SmoothingBivariateDataExample')

    See also spap2, spaps, tpaps

    Introduced in Curve Fitting Toolbox before R2006a
    Documentation for csaps
       doc csaps

Instead, you can allow csaps to choose a value that it thinks is better, for this set of data.. It chose a value very close to 1.

[fitt2,P] = csaps(x, y);
format long g
P

P = 
         0.999991231058044

figure
fplot(@(X) ppval(fitt2,X),[min(x),max(x)])
hold on
plot(x,y,'bo')

And here we see a resaonable fit, that passes smoothly through the data. If you felt that bump in the bottom is still too much, then you could fix that, by using a slightly smaller smoothing parameter.

fitt3 = csaps(x, y,0.99995);
figure
fplot(@(X) ppval(fitt3,X),[min(x),max(x)])
hold on
plot(x,y,'bo')

Finally, you complain that the curve drops to zero outside of the bounds.

fplot(@(X) ppval(fitt3,X),[0.75 1.25])
hold on
plot(x,y,'bo')

As you can see, it does not. ppval is able to extrapolate the curve, but extrapolating a spline is a truly dangerous thing to do.

[Renate]

thank you very much for your answer, i tried different smoothing params, up to 1.00, but also one has this downward sloping tail in the lowest values. i will check the ppval ability to extrapolate the curve, thank you very much again.

[John D'Errico]

That is a completely different question, how to extrapolate a spline well. Thankfully, MATLAB has the FNXTR utility. Lets see how this works on your problem.

help fnxtr

 fnxtr - Extrapolate spline
    This MATLAB function returns a spline of order order that extrapolates
    the spline f.

    Syntax
      pp = fnxtr(f,order)
      pp = fnxtr(f)

    Input Arguments
      f - Spline to extrapolate
        structure
      order - Order of extrapolating spline
        integer | vector of integers

    Output Arguments
      pp - Spline in ppform
        spline structure

    Examples
      openExample('curvefit/ExtrapolateCubicSmoothingSplineExample')
      openExample('curvefit/ExtrapolateBivariateBSplineExample')

    See also ppmak, spmak, fn2fm

    Introduced in Curve Fitting Toolbox in R2006a
    Documentation for fnxtr
       doc fnxtr

x = [0.858734025700000
    0.975248753121875
    1.02525125312812
    1.13140821289062
    1.15927407430000
    1.18768630564687
    1.21665290240000
    1.27628156250000];
y = [0.0255926963287944
    0.0206749804955220
    0.0216695267388959
    0.0221177913199036
    0.0248476224632830
    0.0277388242602486
    0.0305936964241147
    0.0362101199563917];
fitt3 = csaps(x, y, 0.99995);

I could specify a constant extrpolant. Sometimes even a linear extrapolant is too much. Be careful, as it is easy to just assume that 1 yields a linear extrapolant.

fittextr = fnxtr(fitt3,1); % A constant extrapolant
fplot(@(X) fnval(fittextr,X),[-1,2])
hold on
plot(x,y,'bo')
hold off
grid on

Here I'll use a linear extrapolant, which is what I presume you would prefer.

fittextr = fnxtr(fitt3,2) % A linear extrapolant

fittextr = struct with fields:
      form: 'pp'
    breaks: [-0.1413 0.8587 0.9752 1.0253 1.1314 1.1593 1.1877 1.2167 1.2763 2.2763]
     coefs: [9x4 double]
    pieces: 9
     order: 4
       dim: 1

fplot(@(X) fnval(fittextr,X),[-1,2])
hold on
plot(x,y,'bo')
grid on

If I just throw a spline directly into ppval or fnval without intervenign with fnextr, then it uses the end segments of the spline, and extrapolates them. And this is what you did not like.

Answer 2

[Mathieu NOE]

hello

I interpreted your question as how to fit a U shaped function to your data , so at the end I opted for a parabola , is that what you wanted ?

% Given Data
x =[0.858734025700000
0.975248753121875
1.02525125312812
1.13140821289062
1.15927407430000
1.18768630564687
1.21665290240000
1.27628156250000];

y = [0.0255926963287944
0.0206749804955220
0.0216695267388959
0.0221177913199036
0.0248476224632830
0.0277388242602486
0.0305936964241147
0.0362101199563917];

% interpolation or smoothing or fit a parabola ? 
p = polyfit(x,y,2);  % fit a parabola

% Evaluate the fitted polynomial p and plot:
x_values_finer = linspace(0.5904900,2.0112400,100);
f = polyval(p,x_values_finer);
plot(x,y,'o',x_values_finer,f,'-')
legend('data','fit')

[Renate]

Thank you very much, indeed the fit looks nice here, it is just i was after the same shape first derivative like in R (i attached a figure to my question). also, i think i need to keep the flexibility of splines for other input data. the first derivate of this parabola is a straight line.

Answer 3

[Pavl M.]

clc;
clear global;
tol = 1e-07;
format('long','g') %format of the numeric values: each number represented by about 15 digits
%format native-bit
format longg
rand('state',1)
function stasis_level = collate(params)
%ToDo: may be continued if you hire me
end
%% Measurement(Experiment, Lab) Input Datum:
x_val = [0.858734025700000 0.975248753121875 1.02525125312812 1.13140821289062 1.15927407430000 1.18768630564687 1.21665290240000 1.27628156250000];
y_val = [0.0255926963287944 0.0206749804955220 0.0216695267388959 0.0221177913199036 0.0248476224632830 0.0277388242602486 0.0305936964241147 0.0362101199563917];
%% Real world natural datum to extrapolate to:
x_values_finer = (1+(-0.1))^5:(0.15-(-0.1))/1000:(1+0.15)^5;
%% System/function identification by polynomial interpolation/fit:
PP = spline (x_val, y_val)
%% Interpolated system application on real world natural datum(Extrapolation):
YI = spline (x_val, y_val, x_values_finer);
%% Validation:
y_int = ppval (PP, x_val);
%% Resultants:
figure
plot(x_val,y_val)
title('yval vs xval')
figure
plot(x_values_finer)
title('X Values Finer')
figure
plot(x_val,y_int)
title('Y computed from Spline interp.')
figure
plot(x_val,y_val,'o',x_values_finer,YI,'-')
title('Pol1')
figure
plot(YI)
title('Interpolated with Spline values')
% Calculate derivative:
dfdx1 = fnder(PP,1)
dfdx2 = fnder(PP,2)
%Compute derivative:
der1 = ppval(dfdx1, x_values_finer);
der2 = ppval(dfdx2,x_values_finer);
%Plot derivative:
figure
plot(x_values_finer,der1,'-b',x_values_finer, der2,'-r');
title('1st approx. derivative order 1 und 2')
%Constructed by P.Mazniker:
%https://independent.academia.edu/PMazniker
%[phone]
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://join.skype.com/invite/iEiqnsfbpgF8
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
%The End.

[Renate]

thank you very much, i also tried all of these options and got these (to me unfortunately unsatisfying ) results. :-/

[Pavl M.]

Confidential Solution.
Of course it's true, correct, right, proven by facts and sci-tech processes. My results are satisfyying, yup ok. Data analysis: Spline is cubic and as accurate polyfit of order 2, because of the input data, more input will lead to the distinguishing what is more optimal 2 or 3 order of polynomial (spline or quadratic). I have updated, spline and polynomial derivative. It is just left to find derivative of polynomial by only coefficients and not like spline when you find derivative by the full structure(made already) of order in TCE with builtin or workaround functions ... m working on it ... I hope I answered your question.
Here it is the derivative values I found first (for spline and polynome):

dfdx1p = polyder(p2)
dfdx1p =
0.4571 -0.4629

here is the replete, outright answer from me. Please accept.

%Begin:

clc
clear all
close all

%% Measurement(Experiment, Lab) Input Datum:

x_val = [0.858734025700000 0.975248753121875 1.02525125312812 1.13140821289062 1.15927407430000 1.18768630564687 1.21665290240000 1.27628156250000];
y_val = [0.0255926963287944 0.0206749804955220 0.0216695267388959 0.0221177913199036 0.0248476224632830 0.0277388242602486 0.0305936964241147 0.0362101199563917];

%% Real world natural datum to extrapolate to:

x_values_finer = (1+(-0.1))^5:(0.15-(-0.1))/1000:(1+0.15)^5;
%equivalent to:
%x_values_finer = linspace(0.5904900,2.0112400,100);

%% System/function identification by polynomial interpolation/fit:

PP = spline (x_val, y_val)
PP2 = pchip(x_val, y_val)
p2 = polyfit(x_val,y_val,2); % fit a parabola
p22 = mkpp(p2,x_val);
%p22 = polyder(p2);
%p9 = polyfit(x_val,y_val,9);
%% Interpolated system application on real world natural datum(Extrapolation):
YI = spline (x_val, y_val, x_values_finer);
%% In order to compute derivative:
bk = PP.breaks;
cf = PP.coefs;
[m,n] = size(cf);
syms z
v = z.^((n-1:-1:0)');
eq = cf*v;
%% Calculate derivative:
dfdx1 = fnder(PP,1);
dfdx2 = fnder(PP,2);
dfdx1p = polyder(p2);
dfdx2p = polyder(dfdx1p);

%Compute derivative:
der1 = ppval(dfdx1, x_values_finer);
der2 = ppval(dfdx2,x_values_finer);
%Plot derivative:
figure
plot(x_values_finer,der1,'-b',x_values_finer, der2,'-r');
title(' Spline derivative order 1 & 2')

%Compute derivative:
der2p = dfdx2p*ones(1,length(x_values_finer));%= ppval(mkpp(dfdx2p,x_val), x_values_finer);
der1p = x_values_finer*dfdx1p(1) + dfdx1p(2); %ppval(mkpp(dfdx1p,x_val), x_values_finer);

%Plot derivative:
figure
plot(x_values_finer,der1p,'-b',x_values_finer, der2p,'-r');
title('Quadratic poly derivative order 1 & 2')

YI1 = pchip(x_val,y_val,x_values_finer);

f2 = polyval(p2,x_values_finer);
%f9 = polyval(p9,x_values_finer);

%% Validation:
y_int = ppval (PP, x_val);
y_int2 = ppval(PP2,x_val);
y2 = polyval(p2,x_val);
%y9 = polyval(p9,x_val);

%% Resultants:
figure
plot(x_val,y_val)
title('yval vs xval')

figure
plot(x_values_finer)
title('X Values Finer')

figure
plot(x_val,y_int)
title('Y computed from Spline interp.')

figure
plot(YI)
title('Interpolated with Spline values')

figure
plot(x_val,y_int2)
title('Y computed from Piecewise Cubic Hermite Interpolating Polynomial interp.')

figure
plot(YI1)
title('Interpolated with Piecewise Cubic Hermite Interpolating Polynomial values')

figure
plot(x_val,y_val,'o',x_values_finer,f2,'-')
title('Pol1')
%figure
%plot(x_val,y_val,'o',x_values_finer,f9,'*')
%title('Pol2')
figure
plot(x_val,y2,'.')
title('Pol3')
%figure
%plot(x_val,y9,'+')
%title('Pol4')

%Developed from needing help code by
%https://independent.academia.edu/PMazniker
%[phone]
%https://join.skype.com/invite/oXnJhbgys7oW
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w

%The End.