python - 曲线拟合幂律到双对数数据

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循环中，根据变量具有不同的未对齐分数，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)

标签： pythonnumpycurve-fittingscipy-optimizepower-law