Azzera filtri
Azzera filtri

improve fit compared to R

97 visualizzazioni (ultimi 30 giorni)
Renate
Renate il 11 Ott 2024 alle 13:11
Modificato: Pavl M. il 11 Ott 2024 alle 17:59
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:
  1 Commento
Renate
Renate il 11 Ott 2024 alle 14:27
this nice thing i get in R

Accedi per commentare.

Risposte (3)

John D'Errico
John D'Errico il 11 Ott 2024 alle 14:11
Modificato: John D'Errico il 11 Ott 2024 alle 14:21
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.
  2 Commenti
Renate
Renate il 11 Ott 2024 alle 15:04
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
John D'Errico il 11 Ott 2024 alle 15:29
Modificato: John D'Errico il 11 Ott 2024 alle 15:29
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.

Accedi per commentare.


Mathieu NOE
Mathieu NOE il 11 Ott 2024 alle 13:51
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')
  1 Commento
Renate
Renate il 11 Ott 2024 alle 15:00
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.

Accedi per commentare.


Pavl M.
Pavl M. il 11 Ott 2024 alle 14:36
Modificato: Walter Roberson il 11 Ott 2024 alle 14:59

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
%+380990535261
%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.

  2 Commenti
Renate
Renate il 11 Ott 2024 alle 15:01
thank you very much, i also tried all of these options and got these (to me unfortunately unsatisfying ) results. :-/
Pavl M.
Pavl M. il 11 Ott 2024 alle 16:43
Modificato: Pavl M. il 11 Ott 2024 alle 17:59

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
%+380990535261
%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.

Accedi per commentare.

Categorie

Scopri di più su Spline Postprocessing in Help Center e File Exchange

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by