Power function fitting with my data

18 次查看（过去 30 天）
Ahmed Abdalazeez
Ahmed Abdalazeez 2018-6-30
评论： Ahmed Abdalazeez 2018-7-2
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,
Ahmed
  12 个评论
显示 10更早的评论
Image Analyst
Image Analyst 2018-7-2
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.

请先登录，再进行评论。

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

采纳的回答

John D'Errico
John D'Errico 2018-7-1
编辑：John D'Errico 2018-7-1
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
a*x^b
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.
  2 个评论
John D'Errico
John D'Errico 2018-7-2
编辑：Image Analyst 2018-7-2
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)
plot(LL,eta_norm,'o')
hold on
plot(mdl)
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.

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2018-7-1
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')

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by