Estimating the parameter of a power law distribution using maximum likelihood estimation (MLE)

I have a set of data points. I want to fit a power law distibution to these data points using the relation:

P(E) = CE^-m

where,

E = Set of data points

P(E) is the probability density function of E

C = a constant

m = power law exponent

I need to determine "C" and "m" using MLE and overlap the power law curve over P(E) versus E plot.

Please help.

Answer 1

[Animesh]

Hey @@Kashif Naukhez

To fit a power law distribution to your data points using Maximum Likelihood Estimation (MLE) in MATLAB, you can follow these steps:

1.Define the power law function powerLaw with parameters C and m:

powerLaw = @(C, m, E) C * E.^(-m);

2.Create a probability density function (PDF) of your data points E using the power law function:

pdf = @(C, m) powerLaw(C, m, E);

3.Define the negative log-likelihood function negLogLikelihood to be minimized:

negLogLikelihood = @(params) -sum(log(pdf(params(1), params(2))));

4.Use the fminsearch function to find the parameters C and m that minimize the negative log-likelihood:

initialGuess = [1, 1];  % Initial guess for C and m
params = fminsearch(negLogLikelihood, initialGuess);
C = params(1);
m = params(2);

5.Plot the data points E and the fitted power law curve:

scatter(E, P, 'b', 'filled');
hold on;
x = linspace(min(E), max(E), 100);
plot(x, powerLaw(C, m, x), 'r', 'LineWidth', 2);
xlabel('E');
ylabel('P(E)');
legend('Data', 'Power Law Fit');

Make sure to replace E and P with your actual data points and their corresponding probabilities.

[Kashif Naukhez]

I am having trouble creating the functions as I am new to MATLAB and getting some errors while running the code. Can u send the full code, if possible

Answer 2

[Shantanu Dixit]

Hi Kashif,

You can use fminsearch - MATLAB and obtain the C and m for the power law fit. See the below snippet for minimizing the negative log likelihood using fminsearch.

negLogLikelihood = @(params) -sum(log(params(1) * data.^(-params(2))));
initialGuess = [0.5,0.5];   % initial guess for c and m
options = optimset('MaxFunEvals', 1000); 
params = fminsearch(negLogLikelihood, initialGuess, options);
C = params(1);
m = params(2);