[Python] Curvefitting powerlaw to double-log data

I'm trying to fit a power-law to data which is in the double log scale. Therefore I've used the curve_fit(...) function from the scipy.optimize package.

To run the function I've implemented the following piece of code COR_coef = curve_fit(lambda x, m: c * x ** m, x, COR_IFG[:, i])[0][0], to the best of my knowledge the curve_fit(...) should now correctly fit a power-law (being a straight line) to my data.

However, for some reason, I just do not seem to get the fit right. See the attached picture for the data and its fit.



Some more context with regards to the minimum reproducible example (see below):

• The code generates random noise for simulation purposes, this is done in the white_noise(...)
• This random noise is than misaligned (in a for-loop with different fractions of misalignment according to the variable fractions_to_shift so the development of the power-law can be studied) and subtracted from the original noise to gain a residual signal
• The residual signal is the signal the power-law is fitted to
• The curve_fit(...) is applied in the sim_powerlaw_coefficient(...) function
• I am aware of the fact that my residual signal shows some artifacts when the misalignment gets larger, unfortunately I don't know how to get rid of these artifacts.

MINIMUM REPRODUCIBLE EXAMPLE

import matplotlib.pyplot as plt
import numpy as np
import numpy.fft as fft
import numpy.random as rnd
from scipy.optimize import curve_fit

plt.style.use('seaborn-darkgrid')

rnd.seed(100) # to select a random seed for creating the "random" noise
grad = -5 / 3. # slope to use for every function
c = 1 # base parameter for the powerlaw
ylim = [1e-7, 30] # range for the double log plots of the powerfrequency domains
values_to_shift = [0, 2**-11, 2**-10, 2**-9, 2**-8, 2**-7, 2**-6, 2**-5, 2**-4, 2**-3, 2**-2, 2**-1, 2**0] # fractions of misalignment

def white_noise(n: int, N: int):
"""
- Creates a data set of white noise with size n, N;
- Filters this dataset with the corresponding slope;
This slope is usually equal to -5/3 or -2/3
- Makes sure the slope is equal to the requested slope in the double log scale.

@param n: size of random array
@param N: number of random arrays
@param slope: slope of the gradient
@return: white_noise, filtered white_noise and the original signal
"""

m = grad
x = np.linspace(1, n, n // 2)
slope_loglog = c * x ** m

whitenoise = rnd.randn(n // 2, N) + 1j * rnd.randn(n // 2, N)
whitenoise[0, :] = 0 # zero-mean noise
whitenoise_filtered = whitenoise * slope_loglog[:, np.newaxis]

whitenoise = 2 * np.pi * np.concatenate((whitenoise, whitenoise[0:1, :], np.conj(whitenoise[-1:0:-1, :])), axis=0)
whitenoise_filtered = 2 * np.pi * np.concatenate(
(whitenoise_filtered, whitenoise_filtered[0:1, :], np.conj(whitenoise_filtered[-1:0:-1, :])), axis=0)

whitenoise_signal = fft.ifft(whitenoise_filtered, axis=0)
whitenoise_signal = np.real_if_close(whitenoise_signal)
if np.iscomplex(whitenoise_signal).any():
print('Warning! whitenoise_signal is complex-valued!')
whitenoise_retransformed = fft.fft(whitenoise_signal, axis=0)

return whitenoise, whitenoise_filtered, whitenoise_signal, whitenoise_retransformed, slope_loglog


def sim_powerlaw_coefficient(n: int, N: int, show_powerlaw=0):
"""
@param n: Number of values in the IFG
@param N: Number of IFG's
@return: Returns the coefficient after subtraction of two IFG's
"""

master = white_noise(n, N)
slave = white_noise(n, N)
x = np.linspace(1, n, n // 2)

signal_IFG = master[2] - slave[2]
noise_IFG = np.abs(fft.fft(signal_IFG, axis=0))[0:n // 2, :]

for k in range(len(values_to_shift)):
shift = np.int(np.round(values_to_shift[k] * n, 0))
inp = signal_IFG.copy()

# the weather model is a shifted copy of the actual signal, to better understand the errors that are introduced.
weather_model = np.roll(inp, shift, axis=0)
WM_IFG = np.abs(fft.fft(weather_model, axis=0)[0:n // 2, :])

signal_corrected = signal_IFG - weather_model
COR_IFG = np.abs(fft.fft(signal_corrected, axis=0)[0:n // 2, :])

COR_coef = np.zeros(N)

for i in range(N):
COR_coef = curve_fit(lambda x, m: c * x ** m, x, COR_IFG[:, i])[0][0]

plt.figure(figsize=(15, 10))

plt.title('Corrected IFG (combined - weather model)')
plt.loglog(COR_IFG, label='Corrected IFG')
plt.ylim(ylim)
plt.xlabel('log(k)')
plt.ylabel('log(P)')

plt.loglog(c * x ** COR_coef.mean(), '-.', label=f'COR powerlaw coef:{COR_coef.mean()}')

plt.legend(loc=0)
plt.tight_layout()

sim_powerlaw_coefficient(8192, 1, show_powerlaw=1)

Continue reading...