为什么在对数-对数刻度中进行线拟合时，我的线拟合会导致 np.nan 的梯度？
Why does my line fit result in a gradient of np.nan when line fitting in Log-Log scale?
我试图找到在双对数刻度线拟合时绘制图形的梯度。因此我写了下面的函数。
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 后应用的值（在所述斜率系数下过滤噪声）。当为 operations 的数量输入 0 时，将执行移位，这将导致两个系数的数组相同。但是，它会导致原始噪声的 nan 和移位噪声的实数值（在这种情况下根本没有移位，因此是原始噪声）。
谁能告诉我我做错了什么？
NOTE: you may assume that the shift operation is working properly since I've been debugging that one for a while now and it behaves as it should. This question is purely about my line fitting function.
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)
经过一番研究，我发现 slope_loglog 定义不正确。它的定义方式对于在双对数图中绘制直线是正确的，但由于我正在研究幂律行为，所以我需要使用幂律公式。因此，slope_loglog 变为 c * x ** m，其中 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)}')
我试图找到在双对数刻度线拟合时绘制图形的梯度。因此我写了下面的函数。
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 后应用的值（在所述斜率系数下过滤噪声）。当为 operations 的数量输入 0 时，将执行移位，这将导致两个系数的数组相同。但是，它会导致原始噪声的 nan 和移位噪声的实数值（在这种情况下根本没有移位，因此是原始噪声）。
谁能告诉我我做错了什么？
NOTE: you may assume that the shift operation is working properly since I've been debugging that one for a while now and it behaves as it should. This question is purely about my line fitting function.
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)
经过一番研究，我发现 slope_loglog 定义不正确。它的定义方式对于在双对数图中绘制直线是正确的，但由于我正在研究幂律行为，所以我需要使用幂律公式。因此，slope_loglog 变为 c * x ** m，其中 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)}')