Power law curve help

line([0 3 6],[0 4244 6000],'color','m')

I have to make a power law curve on a graph from points (0,0) to (6,6000) and it has to go through (3,4244).

How would I write this as code?

What i have above isnt right...

Thanks

[J. Alex Lee]

forget matlab and coding for now. how would you define this power law curve based on the points it has to go through? do you have enough information? do you have too much?

[Zoe Jackson]

I have no idea. I think I have the needed info since I'm expected to be able to do this though.

[J. Alex Lee]

What is the form of the power law you want (how many parameters)?

[Zoe Jackson]

I have the starting point, ending point, and a point it must pass through.

its the basic y=bx^m

Though I think I've figured out some of it by going through the exam review answer key (the one i need this for is a different assignment but...)

Apparently you use polyfit command to find m and b... plus log10? (I missed the day this was explained in class and no one helped fill me in and the notes are fairly useless so I come here for help)

[Image Analyst]

You can get a fairly good answer that way, but probably not the best, most accurate solution as if you'd used fitnlm (Fit NonLinear Model). Did you even see or try my solution below? You didn't respond at all to it so I'm not sure.

[Zoe Jackson]

I did look at it, it's a bit hard for me to follow (might just be cause I'm sleep deprived) but i don't think we're suppoesd to use the Fit NonLinear Model.

It's harder for me to say cause, again, i missed that day in class but we're supposed to do it in the method taught.

But less accurate is fine for now as long as i get the right answer.

Thanks for the help though!

[Image Analyst]

You need more points because there is no curve that is a power law y = b * x .^ m that can go through (0,0). The only way to get that is if b = 0, and then you essentially have no equation at all. Do you have more data?

[Zoe Jackson]

I have starting x point, ending x point, starting y, ending y, and a middle x and y point

oh, the 0,0 is the starting x and y, where the curve begins

[Image Analyst]

Yes, that may be, but my comment stands. No power law of the formula you gave will go through (0,0).

[J. Alex Lee]

As long as the power law is as simple as you say,

Then you can indeed fit a linear law to the resultant equation of taking log of both sides

I'm not sure why Image Analyst says that will not give "as good" an answer as doing a nonlinear fit on the original variable, although I guess the procedures are not equivalent, i.e., errors would not be the same in the respective "best"-ness of fit?

But for sure such an expression cannot go through (0,0).

I don't know if there's a strict definition of "power law", but consider

Now what would you do?

[Zoe Jackson]

At this point I've had to turn in the assignment, but the study guide for the exam showed me the method to use for future reference.

so id use polyfit, then split the answer into the needed numbers, then plug it into the equation.

Thanks for helping though.

[Image Analyst]

In the future, put a "tag" on your post if it's homework. We don't want anyone to get into trouble with their instructors should they turn in someone else's online code as their own. Tagging it as homework, lets us know to give hints, not complete solutions so you can't get into trouble.

[Zoe Jackson]

Okay, didn't think of that, thanks for the advice.

Answer 1

[Image Analyst]

Here is the method using polyfit on the transformed equation, followed by the method of using fitnlm(). I had to make the x(1) and y(1) slightly non-zero to prevent the process from giving a bunch of nan. But as you can see the polyfit method fit is not as good as with fitnlm():

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;

% I have to make a power law curve on a graph from points (0,0) to (6,6000) and it has to go through (3,4244). 
x = [0.0001 3 6]
y = [0.0001 4244 6000]
plot(x, y, 'b.', 'MarkerSize', 50);
grid on;
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);

% Use polyfit to estimate coefficients.
% Residuals using this method may be higher than if you'd used fitnlm().
% y=bx^m
% log(y) =  m * log(x) + log(b)
logy = log(y);
logx = log(x);
coefficients = polyfit(logx, logy, 1)
% So now, coefficients(1) = m, and coefficients(2) = log(b)

% Now get the fit using the computed coefficients.
% Method 1: directly using formula:
b = exp(coefficients(2));
m = coefficients(1);
x1 = min(x)
x2 = max(x)
xFit = linspace(x1, x2, 500);
y2 = b * xFit .^ m;
hold on;
plot(xFit, y2, 'r-', 'LineWidth', 7);

% Method 2: using polyval
% Now make a fit from left to right with 1000 points.
yLogFit = polyval(coefficients, log(xFit))
% This gives an estimate of log(y), not y itself.  Exponentiate to get y
yFit = exp(yLogFit);
hold on;
plot(xFit, yFit, 'g-', 'LineWidth', 2);

%=======================================================================================================================
% Uses fitnlm() to fit a non-linear model (a power law curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.

% 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) * x(:, 1) .^ b(2) + b(3);  
beta0 = [2000, .4, -20]; % Guess values to start with.  Just make your best guess.  They don't have to match the [a,b,c] values from above because normally you would not know those.
% 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) * x .^ coefficients(2) + coefficients(3);
% Do another fit but for a lot more points, including points in between the training points.
X1000 = linspace(x(1), x(end), 1000);
yFitted1000 = coefficients(1) * X1000 .^ coefficients(2) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
% Plot fitted values at all the 1000 X values with a red line.
plot(X1000, yFitted1000, 'c-', 'LineWidth', 4);
grid on;
title('Power Law Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Training Data', 'y2 = b * x .^ m', 'yFit = exp(yLogFit)', 'fitnlm', 'Location', 'northwest');

legendHandle.FontSize = 25;

message = sprintf('Coefficients for Y = a * X ^ b + c:\n  a = %8.5f\n  b = %8.5f\n  c = %8.5f',...
coefficients(1), coefficients(2), coefficients(3));
text(0.5, 5500, 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') 

[J. Alex Lee]

The power law form is sufficient, as it is guaranteed to go through (0,0), so there would only be 2 remaining usable points for the power law fit, and only 2 parameters a and b. Thus:

clc;
close all;
clear;

xData = [0;6;3]
yData = [0;6000;4244]

% assume the law y = a x^b
% so log(y) = b log(x) + log(a)
% the point (0,0) is redundant, so we don't include in fit
z = log(yData(2:end));
w = log(xData(2:end))

% do the linear least squares on the log of variables
p = polyfit(w,z,1)

% the exponent b is the order 1 coefficient
b_num = p(1)

% the prefactor a is the exponential of the order 0 coefficient
a_num = exp(p(2))

% exact solution
log_a = (log(4244)/log(3) - log(6000)/log(6))/(1/log(3)-1/log(6))
a = exp( log_a )
b = log(4244)/log(3) - log_a/log(3)

% here is the fitted function
fitfun = @(x) a*x.^b

figure(1); cla; hold on
plot(xData,yData,'+','LineWidth',2)
fplot(fitfun,[-1,7],'LineWidth',2)

text(2,2000,["$y = a x^b - 1$";"a = "+sprintf("%g",a);"b = "+sprintf("%g",b)],"Interpreter","latex")

Where you can compare the numerical fit to the analytical exact expressions for the fit parameters as I have shown in the code.

[Image Analyst]

Yes, that -25.4 popped out of fitnlm() as the best offset to give the best fit. If I wanted to I could have put in -1 or any other number and the other coefficients would have adjusted for the new offset but I assume the residuals would have been higher than if I let the program just decide on its own what the optimal offset is. So you can see with your offset of -1, the values for a and b are very close to mine (because -1 or -25 is not much different when the other numbers are several thousand).

[J. Alex Lee]

My bad, Image Analyst you threw me off with the statement that a power law can't go through (0,0)...that sounded right, but that's for exponential forms. Not only can have go through (0,0), it must!

Editing my comment above to match.

Answer 2

[Image Analyst]

See attached demo.

[Faisal]

is there a way to write the code to find the exponent in power law function? assume the same function of y = a*x^b. How can I find b? if I have a data of about 1000 points?

[Image Analyst]

Yes. Did you run my demo? It fits y=a*x^b + c. If you want c to be zero, then just delete c everywhere you see c in my code. I've donke that below. If you have trouble, start your own, new question with your code and data attached.

% Demo by Image Analyst to fit data to a power law y = a * x ^ b
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 short g;
format compact;
fontSize = 22;
markerSize = 6;
% Uses fitnlm() to fit a non-linear model (a 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;

%-----------------------------------------------------------------------------------------------------------------------
% FIRST CREATE X AND Y DATA.  SKIP THIS IF YOU ALREADY HAVE DATA AND USE YOUR OWN INSTEAD OF THIS!
% Create the X coordinates: 30 points from 0.01 to 20, inclusive.
X = linspace(0.01, 20, 30);
% Define function that the X values obey.  Define Y as a function of X.
a = 10 % Arbitrary sample values I picked.
b = 0.4
Y = a * X .^ b; % Get a vector.  No noise in this Y yet.
% Add noise to Y so we don't have such a perfect fit.
Y = Y + 0.8 * randn(1, length(Y));

%-----------------------------------------------------------------------------------------------------------------------
% NOW WE HAVE OUR X AND Y DATA AND WE CAN BEGIN.
% Now we have noisy training Y data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 20);
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)
% 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) * x(:, 1) .^ b(2);  
beta0 = [10, .4]; % Guess values to start with.  Just make your best guess.  They don't have to match the [a,b,c] values from above because normally you would not know those.

% 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) * X .^ coefficients(2);
% Do another fit but for a lot more points, including points in between the training points.
X1000 = linspace(X(1), X(end), 1000);
yFitted1000 = coefficients(1) * X1000 .^ coefficients(2);

% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
% Plot red diamonds fitted values at the training X values.
plot(X, yFitted, 'rd', 'LineWidth', 2, 'MarkerSize', 10);
% Plot fitted values at all the 1000 X values with a red line.
plot(X1000, yFitted1000, 'r-', 'LineWidth', 2);
grid on;
title('Power Law Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Training Y', 'Fitted Y at training X', 'Fitted Y everywhere', 'Location', 'north');
legendHandle.FontSize = 25;
message = sprintf('Coefficients for Y = a * X ^ b:\n  a = %8.5f\n  b = %8.5f',...
    coefficients(1), coefficients(2));
text(8, 15, 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')