Power function fitting with my data

17 views (last 30 days)
I have a question regarding the power function fitting.... I have attached my two pictures. They may help you to understand what I want to do.
1. I am looking for the best fitting with my data (blue curve).
using a power function respect to x. I tried one, but it's not as good as I want.
2. For linear fitting with (blue 'o') respect to x.
Thank you so much in advance.
Best regards,
Image Analyst
Image Analyst on 2 Jul 2018
No problem. It looks like John assumed a model and it works for you since you accepted his answer. So I guess we're done then.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 1 Jul 2018
Edited: John D'Errico on 1 Jul 2018
Since you have not posted your data, it is difficult to help you very much. But...
You can use polyfit to fit the linear curve. However, I see that it is not really a straight line. And you seem to have connected only two points as you drew the line. But just use polyfit to fit a straight line. A quadratic polynomial would be a much better fit however.
Next, you described the first curve as a power function. But if you chose to use the WRONG power function, it will fit terribly poorly. In fact, the curve as it has been drawn does not look like a true power law curve. I'd guess it to be a hyperbolic curve, just looking at the shape of the tails. They look to be asymptotically linear in both tails, and there simply is no way a power fit will be linear in both tails. Anyway, if you think that a pure power fit, thus something like
will fit that first curve, you are simply wrong. I might consider a model of the more general form
a + b*(x-c)^d
Because both a constant offset in x and y will be necessary. It still will not have the correct shape in the upper tail from what I see.
Really, the model I would try fitting to that curve would be something like this:
y = a - b/(x-c)
It won't fit that well, because that curve has a vertical asymptote (so vertical slope on the right end) at x==c. But it would be a start. So then you would next consider a more general hyperbolic curve. I'd suggest the model, if and when you post some data.
So there are several reasons that may explain why you got a poor fit.
So what model were you trying to fit there exactly? Why do you think that curve follows a power law (in any form)?
AND POST YOUR DATA, if you want serious help.
Edit: just for kicks, with no data at all, here is a curve that is starting to look much like your curve.
yhyp = @(x,a,b,x0,y0) y0 + (a*b*(x-x0 + sqrt(- a^2 + b^2 + (x-x0).^2)))/(b^2 - a^2);
ezplot(@(x) yhyp(x,.001,10,265,0.028),[220 290])
It is an arc of a hyperbola, with linear asymptotes as x goes in both directions. The curve fitting toolbox would fit it pretty easily, although I would suggest providing starting values similar to those I used.
John D'Errico
John D'Errico on 2 Jul 2018
Edited: Image Analyst on 2 Jul 2018
Now that you have provided data for the curves y(x), though I note that you still called the variables LL and eta_norm. Sigh. You are asking for trouble in your code by doing things like that.
Using the hyperbolic model I suggested in my answer...
ft = fittype('y0 + (a*b*(x-x0 + sqrt(- a^2 + b^2 + (x-x0).^2)))/(b^2 - a^2)','independent','x','coeff',{'a','b','x0','y0'})
ft =
General model:
ft(a,b,x0,y0,x) = y0 + (a*b*(x-x0 + sqrt(- a^2 + b^2 + (x-x0).^2)))/(b^2 - a^2)
mdl = fit(LL,eta_norm,ft,'start',[.01,10,950,0.028])
mdl =
General model:
mdl(x) = y0 + (a*b*(x-x0 + sqrt(- a^2 + b^2 + (x-x0).^2)))/(b^2 - a^2)
Coefficients (with 95% confidence bounds):
a = 0.002316 (0.002088, 0.002545)
b = 25.34 (23.75, 26.92)
x0 = 945.8 (943.8, 947.9)
y0 = 0.02786 (0.0278, 0.02792)
hold on
I get a very reasonable fit. The horizontal asymptote is at y0 = 0.02786, which seems reasonable from the plot.
You can also compute the linear asymptote on the right side simply enough.
Thus, for large x, we have the limiting behavior:
sqrt(- a^2 + b^2 + (x-x0).^2) --> (x-x0)
So that model reduces, for large x, to:
y0 + 2*a*b*(x-x0)/(b^2 - a^2)
Thus the slope of the linear asymptote is 2*a*b/(b^2-a^2)
In fact, as I predicted, a hyperbolic fit does do quite nicely here, because your data really seems to approach linearity in each wing.
Note that a power function does not fit this data terribly well. In fact, a power curve fit will be more difficult to do, because of numerical issues.

Sign in to comment.

More Answers (1)

Image Analyst
Image Analyst on 1 Jul 2018
Since you're not providing us with data, I had to make code to create some. It's at the top of this script. But it seems to do a pretty good job.
% Uses fitnlm() to fit a non-linear model (an power law curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Create the X coordinates: 30 points from 0.01 to 20, inclusive.
X = linspace(230, 290, 50);
% Define function that the X values obey.
a = .000005;
b = 2 % Arbitrary sample values I picked.
c = .2
d = 225
e = 0.028
Y = a * b .^ (c*( X - d)) + e; % Get a vector. No noise in this Y yet.
% Add noise to Y.
Y = Y + 0.001 * randn(1, length(Y));
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 8);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X', Y');
% Define the model as Y = a * (x .^ b) + c
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * b(2) .^ (b(3)*(x(:, 1) - b(4))) + b(5);
beta0 = [.0001, 2, .2, 225, .028]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * coefficients(2) .^ (coefficients(3)*(X - coefficients(4))) + coefficients(5);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Power Law Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
message = sprintf('Coefficients for Y = a * b .^ (c*(X ^ d)) + e:\n a = %8.5f\n b = %8.5f\n c = %8.5f\n d = %8.5f\n e = %8.5f',...
coefficients(1), coefficients(2), coefficients(3), coefficients(5), coefficients(5));
text(235, .055, message, 'FontSize', 23, 'Color', 'r', 'FontWeight', 'bold', 'Interpreter', 'none');
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

Community Treasure Hunt

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

Start Hunting!