curve fit a custom polynomial

8 次查看（过去 30 天）
Randy Chen
Randy Chen 2021-10-26
评论： Star Strider 2021-10-27
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
Star Strider 2021-10-27
Ran in:
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 个评论
Randy Chen
Randy Chen 2021-10-27
Thank you! I was trying out the code myself, but this error occured:
syms A B C D r z ;
sf = A*r^2 + B*z + C*r + D == 0;
sfiso = isolate(sf,z);
zfcn = matlabFunction(rhs(sfiso), 'Vars',{[A B C D],r});
z = [12.96 46.6 13.5 46.5188]
r = [2.143 102 2.41814 101.5];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0);
Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'});
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')
Unrecognized function or variable 'charcmd'.
Error in sym/isolate (line 108)
throw(CaughtMException);
Error in hw10 (line 4)
sfiso = isolate(sf,z);
Star Strider
Star Strider 2021-10-27
Ran in:
As always, my pleasure!
The isolate function was introduced in R2017a.
Use solve instead —
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
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 个评论
Randy Chen
Randy Chen 2021-10-27
I'm not sure how to rewrite that expression to be pure quadratic. Is it Bz = -Ar^2-Cr-D? But why do I have to rewrite it before using polyfit?
Rik
Rik 2021-10-27
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by