Power function fitting with my data
显示 更早的评论
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 个评论
jonas
2018-7-1
Can you provide the data (or some sample data) and show us what you have tried so far (code)?
John D'Errico
2018-7-1
Please don't tag your posts with the name of a person.
Ahmed Abdalazeez
2018-7-2
John D'Errico
2018-7-2
But what you did not attach is the data for the plot of y versus x. That was the curve where there is some question.
The data you did attach is trivial to fit. Just use polyfit, as I said. Or, if all you want is the curve drawn through your data, then just use the basic fitting tools on the plot itself.
Here you see a linear and quadratic curve fit through the data. It does appear that a quadratic removes the lack of fit in the linear curve. But there is no gain to be seen from a higher order fit.
If you want the coefficients of the polynomial itself, then use polyfit.
But be careful, do NOT just take a 5 significant approximation to those coefficients. You will get an inaccurate model.
polyfit(LL,eta_norm,1)
ans =
-5.8113e-07 0.028593
format long g
polyfit(LL,eta_norm,1)
ans =
-5.81128927834012e-07 0.0285929369824806
The quadratic coefficients:
polyfit(LL,eta_norm,2)
ans =
-2.19797821094579e-10 -3.99697534907371e-07 0.0285661662671785
With coefficients that vary by so many orders of magnitude, you will often need more than just 5 significant digits to adequately evaluate that model.
Again though, the curve that was of interest was the other one. We have seen no data on that. Nor have you explained why it is you think that a power function model is correct. (I doubt that it is correct here, as I argued in my answer.)
Image Analyst
2018-7-2
编辑：Image Analyst
2018-7-2
This data looks more like the second/lower plot rather than like the first/top one that John and I were working on.
Which one do you want? And do you have a model that you want to fit this to? What is the exact form of the fit that you want to use? Do you want a power law y=a^x, or a polynomial y=x^n? Is there a theoretical equation for this data, like one from a textbook?
Ahmed Abdalazeez
2018-7-2
Ahmed Abdalazeez
2018-7-2
John D'Errico
2018-7-2
The "simplest" constant slope is what polyfit(LL,eta_norm,1) returns. The slope is the first element of that vector. But again, I would not that a linear curve seems to show moderately lack of fit there.
Image Analyst
2018-7-2
You say "it is a theoretical equation " but you forgot to share what that equation is. Please give the form of the equation so we can determine the parameters to make it fit best.
Ahmed Abdalazeez
2018-7-2
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.
Ahmed Abdalazeez
2018-7-2
采纳的回答
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
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.
Ahmed Abdalazeez
2018-7-2
更多回答（1 个）
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')
类别
在 帮助中心 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
