curve fit a custom polynomial

4 次查看（过去 30 天）
Randy Chen2021-10-26

I have the following 2nd order polynomial in the r-z coordinates: Right now I have four sets of coordinats (r,z), how should I do a curve fit such that I can get an expression for Z in terms of r?
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];

采纳的回答

Star Strider 2021-10-27
One approach —
syms A B C D r z
sf = A*r^2 + B*z + C*r + D == 0
sf = sfiso = isolate(sf, z)
sfiso = zfcn = matlabFunction(rhs(sfiso), 'Vars',{[A B C D],r})
zfcn = function_handle with value:
@(in1,r)-(in1(:,4)+in1(:,3).*r+in1(:,1).*r.^2)./in1(:,2)
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0)
Warning: The model is overparameterized, and model parameters are not identifiable. You will not be able to compute confidence or prediction intervals, and you should use caution in making predictions.
nlm =
Nonlinear regression model: y ~ F(in1,r) Estimated Coefficients: Estimate SE tStat pValue ________ ________ _______ __________ b1 -0.61998 0.070975 -8.7353 0.072563 b2 2.4979 0.87073 2.8688 0.21353 b3 29.339 1.4571 20.135 0.031591 b4 -275.65 0.16252 -1696.1 0.00037534 Number of observations: 4, Error degrees of freedom: 1 Root Mean Squared Error: 3.89 R-Squared: 0.998, Adjusted R-Squared 0.995 F-statistic vs. constant model: 290, p-value = 0.0415
Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'})
Out = 4×3 table
Parameter Value p-Value _________ ________ __________ {'A'} -0.61998 0.072563 {'B'} 2.4979 0.21353 {'C'} 29.339 0.031591 {'D'} -275.65 0.00037534
rv = linspace(min(r), max(r));
zv = predict(nlm, rv(:));
figure
plot(r, z, 'pg')
hold on
plot(rv, zv, '-r')
hold off
grid
xlabel('r')
ylabel('z')
legend('Data','Model Fit', 'Location','best') The Warning was thrown because the number of parameters are not less than the number of data pairs.
Experiment to get different results.
.
2 个评论显示隐藏 1更早的评论
Star Strider 2021-10-27
As always, my pleasure!
The isolate function was introduced in R2017a.
syms A B C D r z
sf = A*r^2 + B*z + C*r + D == 0
sf = sfiso = solve(sf, z)
sfiso = zfcn = matlabFunction(sfiso, 'Vars',{[A B C D],r})
zfcn = function_handle with value:
@(in1,r)-(in1(:,4)+in1(:,3).*r+in1(:,1).*r.^2)./in1(:,2)
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0)
Warning: The model is overparameterized, and model parameters are not identifiable. You will not be able to compute confidence or prediction intervals, and you should use caution in making predictions.
nlm =
Nonlinear regression model: y ~ F(in1,r) Estimated Coefficients: Estimate SE tStat pValue ________ ________ _______ __________ b1 -0.725 0.082997 -8.7353 0.072563 b2 2.9211 1.0182 2.8688 0.21353 b3 34.308 1.7039 20.135 0.031591 b4 -322.35 0.19005 -1696.1 0.00037534 Number of observations: 4, Error degrees of freedom: 1 Root Mean Squared Error: 3.89 R-Squared: 0.998, Adjusted R-Squared 0.995 F-statistic vs. constant model: 290, p-value = 0.0415
Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'})
Out = 4×3 table
Parameter Value p-Value _________ _______ __________ {'A'} -0.725 0.072563 {'B'} 2.9211 0.21353 {'C'} 34.308 0.031591 {'D'} -322.35 0.00037534
rv = linspace(min(r), max(r));
zv = predict(nlm, rv(:));
figure
plot(r, z, 'pg')
hold on
plot(rv, zv, '-r')
hold off
grid
xlabel('r')
ylabel('z')
legend('Data','Model Fit', 'Location','best') I like isolate because of the output format, and some of its other characteristics.
.

更多回答（1 个）

Rik 2021-10-27
You have two options: rewrite your equation to be a pure quadratic and use polyfit, or use a function like fit or fminsearch on this shape.
If you have trouble implementing either of these two, feel free to comment with what you tried.
2 个评论显示隐藏 1更早的评论
Rik 2021-10-27
Polyfit will fit a pure polynomial of the form
f(x)=p(1)*x^n +p(2)*x^(n-1) ... +p(n)*x +p(n+1)
That means you can determine the values of -A/B, -C/B, and -D/B with polyfit.
As you may conclude from this: there is no unique solution for your setup, unless you have other restrictions to the values you haven't told yet.
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
p=polyfit(r,z,2)
p = 1×3
0.2482 -11.7452 110.3524

Community Treasure Hunt

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

Start Hunting!