How to use Maximum likelihood estimator to find a slope of power law distribution

Question

I have a discrete power spectrum f(x) which is a power law with a specific exponent b: f(x) = A*x^(b). I want to estimate that exponent when I have f(x) and x. I'm interested in using a Maximum Likelihood Estimator (MLE) to find it. Why? Well, there's this paper that says it's the best method for such distributions with least bias: Clauset et al. 2007 .

I already tried many different methods, for instance, I took the commonly known route and took the log for both power and frequency then applied a least squares fit to the loglog spectrum but the slope (sought exponent) of such process is sometimes biased. There's already a Python implementation of the mentioned paper package but it does MLE based on PDF or CDF not the probability distribution you provide. I tried editing it but it's so complicated. I also tried using pytorch to do optimization with no success since I'm not experienced with it.

How can I implement such MLE in python for my problem?

Edit: I have found someone ask a similar question but in R. Someone answered them here but again it's in R not python. Edit: added code

from astroML.time_series import generate_power_law
import numpy as np
import scipy
N=72 # number of points
dt=5*60 #  time resolution
exponent= -2 # power law exponent
rand_seed= 1
time_max = 6*6*60 # observation time
x =  np.linspace(0,  time_max, N ) # time array
y =  generate_power_law(N, dt, -exponent, random_state= rand_seed) 
# amplitudes with power law array
# Get the spectrum by FFT for example
sample_freq = scipy.fft.rfftfreq(1*N, d=dt)[1:] # sampling frequency
sig_fft = scipy.fft.rfft(y,1*N)[1:] # FFT amplitudes
psd = (np.abs(sig_fft)**2)*2*dt/(N) # power spectral density
# do least squares fitting
logA = np.log(sample_freq ) 
logB = np.log(psd )
slope, intercept = np.polyfit(logA, logB, 1)

Another attempt for MLE based on the comment which links article :

import numpy as np
from scipy.optimize import minimize
def func_powerlaw(x, params):
    slope = params[0]
    coefficient = params[1]
    return  (x**slope)*coefficient

params = [1,1]
#For MLE, minimize the negative log likelihood
def neglnlike(params, x, y):
    model = func_powerlaw(x, params)
    output = np.sum(np.log(model) + y/model)
    #Check that this is valid, returning large number if not
    if not np.isfinite(output):
        return 1.0e30
    return output

res = minimize(neglnlike, params, args=(sample_freq , psd),
              method='Nelder-Mead')

Answer 1

See this library named powerlaw

https://github.com/jeffalstott/powerlaw

powerlaw is a toolbox using the statistical methods developed in Clauset et al. 2007 and Klaus et al. 2011 to determine if a probability distribution fits a power law.