Curve fitter limitation for equations in the form f(x,y)=a*x^n/(y+b)^m

Dear Community Members,

When I use curve fitter for equation in the form f(x,y)=a*x^n/(y+b)^m, the results are very bad (R-square is very, very low).

After removing "b" or changing it with a random number (for example 8, then f(x,y)=a*x^n/(y+8)^m), R-square is very good.

The question is: Is this a limitation of the software or there is a way to overcome this problem and to use "b", which value wiil be calculated by Curve fitter.

[Torsten]

In order to avoid negative values for y+b or division by zero during the fitting process, you should use

f(x,y)=a*x^n*((y+b)^2)^m

as equation.

[Sam Chak]

Hi @Boyan,

Nonlinear surface fitting is generally challenging. Could you show us how the surface looks like?

[John D'Errico]

This is often the result of poor starting values. Of course, we don't know your data. We don't know what starting values were used. The mind reading toolbo has been tried, but it failed. Please concentrate, and we will try it again. Or provide the information we need.

[Alex Sha]

@Boyan show your data please.

[Boyan]

Dear all,

The surface looks like fig. 2. In Figure 1 are given just the points. These are three exponential curves, laying in vertical planes.

fig. 1

fig. 2

The data is given below in figure 3. I'm sorry for giving it in picture, not in table, but I don't know how to insert table.

fig. 3

Answer 1

[Sam Chak]

Hi @Boyan

Given that there are only data points, employing a bivariate polynomial of degree two in x and degree four in y will ensure a perfect fit, resulting in an R-squared value of 1.

% Data
x   = [  1   4   5   1   4   5   1   4   5   1   4   5   1   4   5];
y   = [  5   5   5  10  10  10  15  15  15  20  20  20  25  25  25];
z   = [257 410 455 185 306 322 151 246 261 125 191 205 109 159 169];

% p1*x^2 + p2*x + p3
[p21, gof] = fit([1 4 5]', [257 410 455]', 'poly2');
[p22, gof] = fit([1 4 5]', [185 306 322]', 'poly2');
[p23, gof] = fit([1 4 5]', [151 246 261]', 'poly2');
[p24, gof] = fit([1 4 5]', [125 191 205]', 'poly2');
[p25, gof] = fit([1 4 5]', [109 159 169]', 'poly2');

% p1*y^4 + p2*y^3 + p3*y^2 + p4*y + p5
[q14, gof] = fit([5 10 15 20 25]', [p21.p1 p22.p1 p23.p1 p24.p1 p25.p1]', 'poly4');
[q24, gof] = fit([5 10 15 20 25]', [p21.p2 p22.p2 p23.p2 p24.p2 p25.p2]', 'poly4');
[q34, gof] = fit([5 10 15 20 25]', [p21.p3 p22.p3 p23.p3 p24.p3 p25.p3]', 'poly4');

plot3(x, y, z, 'rp'), grid on, hold on
xlabel('x'), ylabel('y'), zlabel('z')

% Surface 
XX      = linspace(x(1), x(end), 101);
YY      = linspace(y(1), y(end), 101);
[X, Y]  = meshgrid(XX, YY);

% z(x, y) = f(y)*x^2 + g(y)*x + h(y)
f   = q14.p1*Y.^4 + q14.p2*Y.^3 + q14.p3*Y.^2 + q14.p4*Y + q14.p5;
g   = q24.p1*Y.^4 + q24.p2*Y.^3 + q24.p3*Y.^2 + q24.p4*Y + q24.p5;
h   = q34.p1*Y.^4 + q34.p2*Y.^3 + q34.p3*Y.^2 + q34.p4*Y + q34.p5;
Z   = f.*X.^2 + g.*X + h;

s  = surf(X, Y, Z, 'FaceAlpha', 0.25); s.EdgeColor = 'none';
title('Surface fitting with bivariate polynomial, poly24(x,y)')
view(145, 30)
hold off

[Alex Sha]

The original fitting function "f(x,y)=a*x^n/(y+b)^m" should produced a very good result like below:

Sum Squared Error (SSE): 388.7688306047
Root of Mean Square Error (RMSE): 5.09096474553825
Correlation Coef. (R): 0.998655962108673
R-Square: 0.9973137306552

Parameter	Best Estimate    
---------	-------------    
a        	11629.1948844926 
n        	0.339456613416582
b        	13.7213095385062 
m        	1.29844412685042 

[Boyan]

Hi Sam,

You are rigtt that bivariate polynomial of degree 2 for x and degree 4 for y gives good results.

But the problem is that I need the formula exact with this structure "f(x,y)=a*x^n/(y+b)^m" (not polynomial). This is because of traditional considerations in a specific engineering field.

Also thank you for the code (I will use it in other situations), but as fas as I have to use "f(x,y)=a*x^n/(y+b)^m" many times, it's better to past trough Curve Fitter App in Matlab.

So the problem with the precense of "b" like a coeffcient in the denominator still exist.

[Boyan]

Hello @Alex Sha,

How did you derive the parameters a, b, n, m in the "f(x,y)=a*x^n/(y+b)^m".

With Curve Fitter with some additional action? This solves my problem, but how did you do that?

[Sam Chak]

Hi @Boyan

For relatively few data points, the bivariate polynomial will guarantee a perfect fit, which is generally effective for prediction within the enclosed boundaries.

If a custom fit model is preferred, you should provide initial guesses for the parameters {a, b, m, n}. If the starting points are relatively close to the solution, they should yield the highest R-squared value.

If you find the codes helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Your support is greatly appreciated!

% your data
X   = [  1   4   5   1   4   5   1   4   5   1   4   5   1   4   5];
Y   = [  5   5   5  10  10  10  15  15  15  20  20  20  25  25  25];
Z   = [257 410 455 185 306 322 151 246 261 125 191 205 109 159 169];

% flattening the matrices for fitting
xData = X(:);
yData = Y(:);
zData = Z(:);

% proposed fit model with parameters a, b, m, n
ft = fittype('(a*x^n)/(y + b)^m', 'independent', {'x', 'y'}, 'dependent', 'z');

% Initial guesses for a, b, m, n
%               a      b   m  n
%               ↓      ↓   ↓  ↓
initialGuess = [12000, 14, 1, 0.5];

% Enter the fit command and view the Goodness of Fit
[fitResult, gof] = fit([xData, yData], zData, ft, 'StartPoint', initialGuess)

fitResult = 
     General model:
     fitResult(x,y) = (a*x^n)/(y + b)^m
     Coefficients (with 95% confidence bounds):
       a =   1.163e+04  (-1.567e+04, 3.893e+04)
       b =       13.72  (3.139, 24.3)
       m =       1.298  (0.7447, 1.852)
       n =      0.3395  (0.3155, 0.3634)

gof = struct with fields:
           sse: 388.7688
       rsquare: 0.9973
           dfe: 11
    adjrsquare: 0.9965
          rmse: 5.9450

% Plot the original data and the fit
figure;
plot(fitResult, [xData, yData], zData);
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
title('Surface Fitting');
view(145, 30)

[Alex Sha]

@Boyan other than Matlab, the result I provided above was obtained by a package named "1stOpt", guessing of initial start-values is no longer required completely, very friendly for all kinds of users.

Not sure how Sam Chak given a perfect guessing of initial start-values like "initialGuess = [12000, 14, 1, 0.5];", which are actually very close to the correct answers I provided before, hardly possible for most of users to give such idea guess.

[Boyan]

Hi @Alex Sha,

I'm sorry but can't see your result in our discussion above. I see only your request to show my data. Maybe it's somehow hidden in the comments above. But anyway the problem is solved by the code of @Sam Chak.

Anyway, because you are mention it, could I ask what is this package "1stOpt"?

[Sam Chak]

Hi @Boyan

I am not affiliated with the product "1stOpt," but it is a very powerful third-party tool, as demonstrated by @Alex Sha in multiple curve-fitting problems. The initial guesses were made after he published his 1stOpt solution. If you are interested, please take a look at the product page:

http://www.7d-soft.com/en/products.htm

I typically guess the starting point parameter values to be close to 0 (to avoid singularity), 1, or the midpoint between 0 and 1, and then gradually increase the magnitudes. In your case, since the z-axis data range is in the hundreds, an educated guess would be to set . See example below:

% your data
X   = [  1   4   5   1   4   5   1   4   5   1   4   5   1   4   5];
Y   = [  5   5   5  10  10  10  15  15  15  20  20  20  25  25  25];
Z   = [257 410 455 185 306 322 151 246 261 125 191 205 109 159 169];

% flattening the matrices for fitting
xData = X(:);
yData = Y(:);
zData = Z(:);

% proposed fit model with parameters a, b, m, n
ft = fittype('(a*x^n)/(y + b)^m', 'independent', {'x', 'y'}, 'dependent', 'z');

% Initial guesses for a, b, m, n
%               a    b  m  n
%               ↓    ↓  ↓  ↓
initialGuess = [100, 1, 1, 1];

% Enter the fit command and view the Goodness of Fit
[fitResult, gof] = fit([xData, yData], zData, ft, 'StartPoint', initialGuess)

fitResult = 
     General model:
     fitResult(x,y) = (a*x^n)/(y + b)^m
     Coefficients (with 95% confidence bounds):
       a =   1.163e+04  (-1.567e+04, 3.893e+04)
       b =       13.72  (3.139, 24.3)
       m =       1.298  (0.7447, 1.852)
       n =      0.3395  (0.3155, 0.3634)

gof = struct with fields:
           sse: 388.7688
       rsquare: 0.9973
           dfe: 11
    adjrsquare: 0.9965
          rmse: 5.9450

% Plot the original data and the fit
figure;
plot(fitResult, [xData, yData], zData);
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
title('Surface Fitting');
view(145, 30)