MATLAB: Power law curve help

graphlinepower law curveStatistics and Machine Learning Toolbox

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

Best Answer

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')
Related Question