Problem with fit Model with odd order polynomial

2 次查看（过去 30 天）
Jason
Jason 2025-7-2
编辑： Matt J 2025-7-3
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--');

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2025-7-2
编辑：Matt J 2025-7-2
Ran in:
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 个评论
显示 4更早的评论
Matt J
Matt J 2025-7-2
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.

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2025-7-2
Ran in:
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 个评论
显示 2更早的评论
Torsten
Torsten 2025-7-2
编辑：Torsten 2025-7-2
@Matt J
Did you look at cond(R) ?
According to the flow chart, "mldivide" uses the QR solver in case of a non-square matrix A:
https://uk.mathworks.com/help/matlab/ref/double.mldivide.html
@Jason
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 2025-7-2
编辑：Matt J 2025-7-3
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by