Using maximum likelihood estimation for power law fitting in Python

powerlaw_fitting_MLE.py

from __future__ import print_function

import numpy as np
import scipy.stats
import matplotlib.pyplot as plt


# Exponent
a = 3.2

# Number of samples
n_samples = 1000


# Generate powerlaw data
data = scipy.stats.powerlaw.rvs(a, loc=0, scale=1, size=n_samples)

# Introduce some gaussian noise
data_noise = data + np.random.normal(0, 0.01, size=n_samples)


### Fit a powerlaw to given data

# Initial estimate of the exponent
exp_est = 3.0

# Initial estimate of x0
x0_est = 0

# Initial estimate of the scale
scale_est = 1


# Perform the fit
pl_fit = scipy.stats.powerlaw.fit(data_noise, exp_est, loc=x0_est, scale=scale_est)


print("Fit:", pl_fit)
###


x_arr = np.linspace(0, 1, 100)

# Plot CDF of the original
plt.plot(x_arr, scipy.stats.powerlaw.cdf(x_arr, a), color='k', linestyle='--', label='Original')
plt.hist(data, cumulative=1, normed=True, histtype='step', color='b', linestyle='--', label='Original')

# Plot CDF of the noisy data and the fit
plt.plot(x_arr, scipy.stats.powerlaw.cdf(x_arr, *pl_fit), color='g', label='Fit')
plt.hist(data_noise, cumulative=1, normed=True, histtype='step', color='r', label='Noise added')



plt.legend()

plt.show()