Problem with fit Model with odd order polynomial

4 views (last 30 days)
Hi, I have the data as below.
Its fits OK'ish to a 3rd order (odd polynomial), but then becomes terrible when trying a 5th order (odd) polynomial
3rd order fit below
and here is the 5th order
Here is my code:
ax=app.UIAxes3;
[x,y] = getDataFromGraph(app,ax,1); % My function which grabs the data from the plot
% 3rd order (odd)
ft=fittype('a*x^3+b*x'); % ft=fittype('a*x^5+b*x^3+c*x');
[fitobj,~,~,~]=fit(x,y,ft); %[fitobj, goodness, output, convmsg]=fit(x,N(:,i),ft)
coeff1=fitobj.a;
coeff2=fitobj.b;
yfit=coeff1*x.^3+coeff2*x;
r=y-yfit;
hold(ax,'on');
plot(ax,x,yfit,'r--');
%5th order (odd)
ft=fittype('a*x^5+b*x^3+c*x');
[fitobj,~,~,~]=fit(x,y,ft); %[fitobj, goodness, output, convmsg]=fit(x,N(:,i),ft)
coeff1=fitobj.a;
coeff2=fitobj.b;
coeff3=fitobj.c;
yfit=coeff1*x.^5+coeff2*x.^3+coeff3*x; %Remember dot notation
plot(ax,x,yfit,'g--');

Accepted Answer

Matt J
Matt J on 2 Jul 2025
Edited: Matt J on 2 Jul 2025
Don't use a custom fit type. Use the builtin poly5 fit type, with bounds on the even degree coefficients.
[x,y] = readvars('DistortionTableData.csv'); % My function which grabs the data from the plot
lb=[-inf,0,-inf,0,-inf,0];
% 3rd order (odd)
ft=fittype('poly3');
[fitobj,~,~,~]=fit(x,y,ft,'Lower',lb(1:4),'Upper',-lb(1:4),'Normalize','on')
fitobj =
Linear model Poly3: fitobj(x) = p1*x^3 + p2*x^2 + p3*x + p4 where x is normalized by mean 1.285e-14 and std 2024 Coefficients (with 95% confidence bounds): p1 = 0.286 (0.2638, 0.3082) p2 = 0 (fixed at bound) p3 = -0.1792 (-0.2227, -0.1357) p4 = 0 (fixed at bound)
plot(fitobj,x,y)
%5th order (odd)
ft=fittype('poly5');
[fitobj,~,~,~]=fit(x,y,ft,'Lower',lb,'Upper',-lb,'Normalize','on')
fitobj =
Linear model Poly5: fitobj(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6 where x is normalized by mean 1.285e-14 and std 2024 Coefficients (with 95% confidence bounds): p1 = 0.04478 (0.01573, 0.07382) p2 = 0 (fixed at bound) p3 = 0.1371 (0.03795, 0.2362) p4 = 0 (fixed at bound) p5 = -0.08362 (-0.1592, -0.008065) p6 = 0 (fixed at bound)
plot(fitobj,x,y)
  6 Comments
Matt J
Matt J on 2 Jul 2025
I changed my mind in light of the conversation in Torsten's answer thread. You should definitely use normalization. I've now edited my original answer accordingly.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 2 Jul 2025
In the case of the polynomial of degree 5, the design matrix is rank-deficient:
T = readmatrix("DistortionTableData.csv");
x = T(:,1);
y = T(:,2);
% 3rd order (odd)
A = [x.^3,x];
rank(A)
ans = 2
b = y;
sol = A\b;
a = sol(1);
b = sol(2);
figure(2)
hold on
plot(x,a*x.^3+b*x)
plot(x,y)
hold off
%5th order (odd)
A = [x.^5,x.^3,x];
rank(A)
ans = 2
b = y;
sol = A\b;
Warning: Rank deficient, rank = 2, tol = 4.322070e+05.
a = sol(1);
b = sol(2);
c = sol(3);
figure(4)
hold on
plot(x,a*x.^5+b*x.^3+c*x)
plot(x,y)
hold off
  4 Comments
Torsten
Torsten on 2 Jul 2025
Edited: Torsten on 2 Jul 2025
Did you look at cond(R) ?
According to the flow chart, "mldivide" uses the QR solver in case of a non-square matrix A:
Whats this method of finding the coefficients called so i can read further (i.e sol=A/b)
It's just solving the overdetermined linear system of equations for the unknown parameter vector in the least-squares sense:
a*x.^5 + b*x.^3 + c*x = y
Matt J
Matt J on 2 Jul 2025
Edited: Matt J on 3 Jul 2025
Did you look at cond(R) ?
We know without looking at cond( R ) that it will be the same as cond(A), but since R is triangular, the error amplification will not be as bad:
[x,~] = readvars('DistortionTableData.csv'); % My function which grabs the data from the plot
c=[1;2;3]; %ground truth coefficients
A = [x.^5,x.^3,x];
[errMLD, errQR]=compareSolvers(A,c)
Warning: Rank deficient, rank = 2, tol = 4.322070e+05.
errMLD = 3.0000
errQR = 0.0017
In any case, the best thing to do here would probably be to normalize the x-data, which improves cond(A) greatly:
x=x/4000;
A = [x.^5,x.^3,x];
cond(A)
ans = 38.6801
[errMLD, errQR]=compareSolvers(A,c)
errMLD = 9.9301e-16
errQR = 1.2608e-14
function [errMLD, errQR]=compareSolvers(A,c)
y=A*c;
[Q,R]=qr(A,0);
errMLD=norm(A\y-c); %error using direct mldivide
errQR=norm(R\(Q'*y)-c); %error using QR
end

Sign in to comment.

Categories

Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange

Tags

Products


Release

R2024b

Community Treasure Hunt

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

Start Hunting!