0

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')
Mohamed
  • 1
  • 1

1 Answers1

0

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.

marwan
  • 504
  • 4
  • 14
  • I literally cited the same package in the question... It only uses a single array and does MLE to find the coefficients of its PDF which is not what I'm looking for. I'm looking for MLE on power spectrum which can be regarded as it's own probability distribution. Please read my question again. – Mohamed Feb 17 '23 at 13:14
  • apologies - did you check this article by any chance? http://keatonb.github.io/archivers/powerspectrumfits – marwan Feb 17 '23 at 18:59
  • Thanks for the reply, the article is really interesting. I already changed their implementation from Lorentzian into power law distribution and added the code to the question. Do you think the to-be-minimized equation " output = np.sum(np.log(model) + y/model) " should be the same even if you change the probability distribution? I tried reading their references and I'm still not sure if it's the correct assumption, nevertheless, it gives promising results already. – Mohamed Feb 17 '23 at 23:07