max value of poynomial surface fit

9 次查看（过去 30 天）
Jason
Jason 2017-6-28
编辑： John D'Errico 2020-5-19
I'm trying to fit a 2D surface and extract the max value from the fit. The fit seems to be working:
x=[1 2 3 1 2 3 1 2 3]'
y=[1 1 1 2 2 2 3 3 3]'
z=[1 3 2 6 10 5 2 4 3]'
f=fit([x,y],z, 'poly22','Normalize','off')
figure
subplot(1,2,1);plot(f, [x,y],z)
So one way I thought I could find the max value was to use the fit at finer sampling, differentiate and find the smallest values in the differentials:
%Differentiate fit with higher x,y sampling:
[xx, yy] = meshgrid( 1:0.1:3, 1:0.1:3 ) %Finer sampling
[fx, fy] = differentiate( f, xx, yy )
then look at each fx, fy and find indicies of lowest values, then get the coordinates using these indices
X=0;
Y=0;
xm=min(abs(fx(:)))
ym=min(abs(fy(:)))
for i=1:numel(fx)
for j=1:numel(fy)
if (fx(i)==xm)&&(fy(j)==ym)
X=xx(i);
Y=yy(j);
disp('found')
[fx(i) fy(j)]
end
end
end
subplot(1,2,2); plot(fx,fy)
subplot(1,2,1); hold on; plot(X,Y,'ro')
But X & Y always come out as 0. I've notices that condition:
if (fx(i)==xm)&&(fy(j)==ym)
is never met and I can't see why. Thanks for any help

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

采纳的回答

Image Analyst
Image Analyst 2017-6-28
Since you're doing it numerically anyway, why not simply use the max() function:
[maxOfF, linearIndexOfMax] = max(f(:));
[row, col] = ind2sub(size(f), linearIndexOfMax)
  5 个评论
显示 3更早的评论
Image Analyst
Image Analyst 2017-6-28
If you want to plot a star in 3-D space, you need to use plot3(), not plot().
v = 1 : 0.1 : 3;
plot3(v(row), v(col), maxOfF, 'r*', 'MarkerSize', 10);

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2017-6-28
编辑：John D'Errico 2017-6-28
Simple enough. Using your data, I'll use my polyfitn tools, found on the file exchange for download. That toolbox has a utility in it that can translate the model into a sym.
P = polyfitn([x,y],z,2)
P =
struct with fields:
ModelTerms: [6×2 double]
Coefficients: [-2.5 1.8821e-15 10.167 -4.5 18.5 -20.667]
ParameterVar: [0.88889 0.44444 16.296 0.88889 16.296 29.432]
ParameterStd: [0.94281 0.66667 4.0369 0.94281 4.0369 5.4251]
DoF: 3
p: [0.076885 1 0.086294 0.017474 0.019509 0.0318]
R2: 0.91111
AdjustedR2: 0.76296
RMSE: 0.7698
VarNames: {'X1' 'X2'}
The same coefficients as you got, but in a different order.
Next, convert it to a symbolic polynomial.
Psym = polyn2sym(P)
Psym =
(61*X1)/6 + (37*X2)/2 + (4771727564350351*X1*X2)/2535301200456458802993406410752 - (5*X1^2)/2 - (9*X2^2)/2 - 62/3
Then differentiate with respect to each variable, and solve.
maxloc = solve(diff(Psym,X1) == 0, diff(Psym,X2) == 0)
maxloc =
struct with fields:
X1: [1×1 sym]
X2: [1×1 sym]
vpa(maxloc.X1)
ans =
2.0333333333333341070915836534138
vpa(maxloc.X2)
ans =
2.0555555555555559807740534792035
The value at that location is
vpa(subs(Psym,{X1,X2},[maxloc.X1,maxloc.X2]))
ans =
8.6833333333333411998755449208187
You should see that the max of the surface is lower than the max datapoint. if you plot the surface and the data (as you already did) you will see that one point stands quite a bit higher than the surface.
No need to do any search or approximation with meshgrid. I'm not sure if the curvefitting TB has a similar utility in it to convert to a sym. If not, I should probably write one someday.
  4 个评论
显示 2更早的评论
Florian Marco Lipp
Florian Marco Lipp 2020-5-19
I converted my model (higher dimension) to a symbolic polynomial. It's not possible to solve it because Operator '==' is not supported for operands of type 'sympoly'. Is there another opportunity to calculate the maximum value of polynomial surface fit?
John D'Errico
John D'Errico 2020-5-19
编辑：John D'Errico 2020-5-19
No, of course you cannot use solve with the sympoly toolbox. I did not write a solve function. Nor does it use the same syntax. HOWEVER, if you read the code I wrote, I did show how to solve the problem.
In what I wrote, I converted the polynomial to a sym, NOT to a sympoly.
Psym = polyn2sym(P)
If you do something sompletely different you must expect it to fail.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by