python - 为什么我的线拟合在对数对数刻度中拟合时会导致 np.nan 的梯度？

Question

我试图找到在双对数刻度线拟合时绘制图表的梯度。因此我写了下面的函数。

def calc_coefficients_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    wnss = white_noise_signal_shift(n, N, num, shift, operations)

    whitenoise_filtered = np.abs(wnss[1][0:n // 2, :]).copy()
    shifted_whitenoise = np.abs(wnss[4][0:n // 2, :]).copy()
    x = np.linspace(1, n, n // 2)

    original_coefficients = []
    shifted_coefficients = []

    for i in range(num):
        original_coefficients.append(np.polyfit(np.log10(x), np.log10(whitenoise_filtered[:, i]), 1))
        shifted_coefficients.append(np.polyfit(np.log10(x), np.log10(shifted_whitenoise[:, i]), 1))

    original_coefficients, shifted_coefficients = \
        np.asarray((original_coefficients, shifted_coefficients))

    fig, (ax1, ax2) = plt.subplots(2, 1, sharey=True, figsize=(10, 7.5))

    ax1.plot(whitenoise_filtered)
    ax1.loglog(10 ** (original_coefficients.mean() * np.log10(x) + 1), 'r')
    ax1.set_title('Original noise')

    ax2.loglog(shifted_whitenoise)
    ax2.loglog(10 ** (shifted_coefficients.mean() * np.log10(x) + 1), 'r')
    ax2.set_title('Shifted noise')

    print(original_coefficients.mean(), shifted_coefficients.mean())

---

该函数的目标calc_coefficients_signal_shift是找出在移动信号后图形的梯度是否发生变化。我希望它在某个地方，-5/3.因为这是我在white_noise使用变量导入函数后应用的值slope_loglog（在所述斜率系数下过滤噪声）。当为执行移位0的次数输入a 时operations，它应该为两个系数生成相同的数组。但是，它会导致nan原始噪声和移位噪声的实际值（在这种情况下根本没有移位，因此是原始噪声）。

有人可以告诉我我做错了什么吗？

注意：您可能会假设移位操作正常工作，因为我已经调试了一段时间，并且它的行为正常。这个问题纯粹是关于我的线拟合功能。

import matplotlib.pyplot as plt
import numpy as np
import numpy.fft as fft
import numpy.random as rnd

# slope to use for every function
grad = -5 / 3.

rnd.seed(10)

def white_noise(n: int, N: int, slope: int = grad):
    x = np.linspace(1, n, n // 2)
    slope_loglog = 10 ** (slope * np.log10(x) + 1)

    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 = np.concatenate((whitenoise, whitenoise[0:1, :], np.conj(whitenoise[-1:0:-1, :])), axis=0)
    whitenoise_filtered = 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 random_arrays(N: int, num: int):
    res = np.asarray(range(N))
    rnd.shuffle(res)

    return res[:num]


def modified_roll(inp, shift: int, operations: int):
    count = 0
    array = inp[:]
    array_rolled = array.copy()

    for k in range(operations):
        count += shift
        array = np.roll(array, shift, axis=0)
        array[:count] = 0
        array_rolled += array

    out = array_rolled / operations

    return out


def white_noise_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    whitenoise, whitenoise_filtered, whitenoise_signal = white_noise(n, N)[:3]

    # only showing the selected arrays
    arrays_to_select = random_arrays(N, num)
    selected_whitenoise = whitenoise[:, arrays_to_select].copy()
    selected_whitenoise_filtered = whitenoise_filtered[:, arrays_to_select].copy()
    selected_whitenoise_signal = whitenoise_signal[:, arrays_to_select].copy()

    # shifting the signal as a field of different refractive index would do
    if operations == 0:
        shifted_signal = selected_whitenoise_signal
    else:
        shifted_signal = modified_roll(selected_whitenoise_signal.copy(), shift, operations)

    # fourier transform back to the power frequency domain
    shifted_whitenoise = fft.fft(shifted_signal, axis=0)

    return selected_whitenoise, selected_whitenoise_filtered, selected_whitenoise_signal, shifted_signal, \
           shifted_whitenoise

def calc_coefficients_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    wnss = white_noise_signal_shift(n, N, num, shift, operations)

    whitenoise_filtered = np.abs(wnss[1][0:n // 2, :]).copy()
    shifted_whitenoise = np.abs(wnss[4][0:n // 2, :]).copy()
    x = np.linspace(1, n, n // 2)

    original_coefficients = []
    shifted_coefficients = []

    for i in range(num):
        original_coefficients.append(np.polyfit(np.log10(x), np.log10(whitenoise_filtered[:, i]), 1))
        shifted_coefficients.append(np.polyfit(np.log10(x), np.log10(shifted_whitenoise[:, i]), 1))

    original_coefficients, shifted_coefficients = \
        np.asarray((original_coefficients, shifted_coefficients))

    fig, (ax1, ax2) = plt.subplots(2, 1, sharey=True, figsize=(10, 7.5))

    ax1.loglog(whitenoise_filtered)
    ax1.loglog(10 ** (original_coefficients.mean() * np.log10(x) + 1), 'r')
    ax1.set_title('Original noise')

    ax2.loglog(shifted_whitenoise)
    ax2.loglog(10 ** (shifted_coefficients.mean() * np.log10(x) + 1), 'r')
    ax2.set_title('Shifted noise')

    print(original_coefficients.mean(), shifted_coefficients.mean())

calc_coefficients_signal_shift(200, 1, 1, 0, 0)

标签： pythonnumpysignal-processingcurve-fitting

Answer 1

经过一番研究，我发现slope_loglog定义不正确。它的定义方式对于在双对数图中绘制直线是正确的，但是由于我正在研究幂律行为，因此我需要使用幂律公式。因此，slope_loglog成为c * x ** m, with c = 10（或我想使用的任何值，因为这不会对结果产生太大影响）并且m是我想要计算的梯度。

用函数实现这个scipy.optimize.curve_fit给了我最终的结果。

def calc_coefficients_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    wnss = white_noise_signal_shift(n, N, num, shift, operations)

    whitenoise_filtered = np.abs(wnss[1][0:n // 2, :]).copy()
    shifted_whitenoise = np.abs(wnss[4][0:n // 2, :]).copy()
    x = np.linspace(1, n, n // 2)

    original_coefficients = []
    shifted_coefficients = []

    for i in range(num):
        original_coefficients.append(curve_fit(lambda x, m: c * x ** m, x, whitenoise_filtered[:, i], p0=(-5/3.))[0][0])
        shifted_coefficients.append(curve_fit(lambda x, m: c * x ** m, x, shifted_whitenoise[:, i])[0][0])

    original_coefficients, shifted_coefficients = \
        np.asarray((original_coefficients, shifted_coefficients))

    fig, (ax1, ax2) = plt.subplots(2, 1, sharey=True, figsize=(10, 7.5))

    ax1.loglog(whitenoise_filtered)
    ax1.loglog(c * x ** original_coefficients.mean(), 'k:', label=f'Slope: {np.round(original_coefficients.mean(), 3)}')
    ax1.set_xlabel('Log(f)')
    ax1.set_ylabel('Log(p)')
    ax1.set_title('Original noise')
    ax1.legend(loc=0)

    ax2.loglog(shifted_whitenoise)
    ax2.loglog(c * x ** shifted_coefficients.mean(), 'k:', label=f'Slope: {np.round(shifted_coefficients.mean(), 3)}')
    ax2.set_ylabel('Log(p)')
    ax2.set_xlabel('Log(f)')
    ax2.set_title('Shifted noise')
    ax2.legend(loc=0)
    plt.tight_layout()

    print(f'The mean of the original coefficients is: {stats.describe(original_coefficients)}')
    print(f'The mean of the shifted coefficients is: {stats.describe(shifted_coefficients)}')