To find exponent in power law equation of the form y = ax^m + b

I have X and Y points for a curve to be of the form Y = ax^m + b.

I want to find the exponent m, lets just say that m could be inbetween 1.2 - 2.5.

How can I find exact value for m?

Answer 1

[Matt J]

fminspleas downloadable from

https://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

is especially appropriate for power law fits.

a = 0.55;
m = 1.3;
b = -0.78;
% dummy data
x = (1:25)';
y = a*x.^m + b  +  randn(size(x));


m=fminspleas( {@(m,x)x.^m  , 1},  2,x,y, 1.2,2.5 )

m = 1.2933

[Faisal]

@Matt J

thank you for your simple code.

I have found the value of m.

Could you please explain that the following code will give the similar answer or it would be different?

file='03_L.xlsx';

[Fs] = xlsread(file,1,'B2:B1216');

[dep] = xlsread(file,1,'A2:A1216');

fcn = @(b,dep) b(1).*dep.^b(2) + b(3); % Objective Function

B = fminsearch(@(b) norm(Fs - fcn(b,dep)), rand(3,1)*2) % Estimate Parameters

figure

plot(dep, Fs, '.')

hold on

plot(dep, fcn(B,dep), '-r')

hold off

grid

[Matt J]

Probably similar, but with 3 unknowns fminsearch is not guaranteed to converge, so no rigorous predictions are possible.

Answer 2

[Mathieu NOE]

hello

try this

may need some refinement for the initial guess for the parameters depending of your data

a = 0.55;
m = 1.3;
b = -0.78;

% dummy data
x = (1:25);
y = a*x.^m + b  +  randn(size(x));

% equation model y = a*x^m + b
f = @(a,m,b,x) (a*x.^m + b);
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);

% IC guessed 
sol = fminsearch(obj_fun, rand(3,1));
a_sol = sol(1)

a_sol = 0.4997

m_sol = sol(2)

m_sol = 1.3288

b_sol = sol(3)

b_sol = -0.4661

y_fit = f(a_sol, m_sol, b_sol, x);

Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

figure(1)
plot(x,y,'rd',x,y_fit,'b-');
title(['Power Fit / R² = ' num2str(Rsquared) ], 'FontSize', 15)
ylabel('Intensity (arb. unit)', 'FontSize', 14)
xlabel('x(nm)', 'FontSize', 14)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R² correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

Answer 3

[Matt J]

If you have the Curve Fitting Toolbox,

a = 0.55;
m = 1.3;
b = -0.78;
% dummy data
x = (1:25)';
y = a*x.^m + b  +  randn(size(x));


fobj=fit(x,y,'power2','Lower',[-inf,1.2,-inf],'Upper',[+inf,2.5,+inf])

fobj = 
     General model Power2:
     fobj(x) = a*x^b+c
     Coefficients (with 95% confidence bounds):
       a =      0.5425  (0.2922, 0.7927)
       b =       1.295  (1.157, 1.434)
       c =     -0.3202  (-1.66, 1.02)

plot(fobj,x,y)