Подбор кривой для двойного журнала данных

Question

Я пытаюсь подобрать степенной закон к данным, которые находятся в двойном логарифмическом масштабе. Поэтому я использовал функцию curve_fit(...) из пакета scipy.optimize. Чтобы запустить функцию, я реализовал следующий фрагмент кода COR_coef[i] = curve_fit(lambda x, m: c * x ** m, x, COR_IFG[:, i])[0][0], насколько мне известно, curve_fit(...) теперь должен правильно соответствовать степенному закону (являющемуся прямой линией) моим данным. Однако по какой-то причине мне кажется, что я не совсем подхожу. См. Прилагаемое изображение для данных и их соответствия.

---

Еще немного контекста относительно минимально воспроизводимого примера (см. Ниже):

• Код генерирует случайный шум для целей моделирования, это делается в white_noise(...)
• Этот случайный шум не совмещен (в for -l oop с различные доли рассогласования в соответствии с переменной fractions_to_shift, чтобы можно было изучить развитие степенного закона) и вычесть из исходного шума, чтобы получить остаточный сигнал
• Остаточный сигнал - это сигнал степенного закона соответствует
• curve_fit(...) применяется в sim_powerlaw_coefficient(...) функции
• Я знаю, что мой остаточный сигнал показывает некоторые артефакты, когда смещение увеличивается, к сожалению, я не Я не знаю, как избавиться от этих артефактов.

---

МИНИМАЛЬНЫЙ ВОСПРОИЗВОДИМЫЙ ПРИМЕР

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 missalignment

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[i] = 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)