Кривая подгонка законодательства о мощности для двойных данных ⇐ Python

[Anonymous]

Я пытаюсь приспособиться к мощности данных, которые находятся в двойном логарифме. Поэтому я использовал функцию curve_fit (...) из пакета scipy.optimize . CODE COR_COEF = curve_fit (lambda x, m: c * x ** m, x, cor_ifg [: i]) [0] [0] , насколько мне известно, curve_fit ( ...) теперь должен правильно соответствовать составлению силы (будучи прямой) для моих данных.
Однако по какой-то причине я просто не получаю подходит правильно. См. Прикрепленное изображение для данных и его подгонки. закон." src = "https://i.sstatic.net/vmjut.png"/>
немного контекста в отношении минимального воспроизводимого примера (см. Ниже): br />
[*] Код генерирует случайный шум для моделирования, это делается в White_noise (...) < /code> < /li>
Это случайное шум, чем смещено (в a for -loop с различными долями смещения в соответствии с переменными fractions_to_shift , чтобы развитие мощного права можно было изучить) и вычитается из исходного шума, чтобы получить остаточный Сигнал
[*] Остаточный сигнал-это сигнал, который установлен с силой, в соответствии с
[*] curve_fit (...) применяется в SIM_POWERLAW_COEFFITH (...) function
Я знаю о том, что мой остаточный сигнал показывает некоторые артефакты, когда смещение становится больше, к сожалению, я не знаю, как это Избавьтесь от этих артефактов.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)



Подробнее здесь: https://stackoverflow.com/questions/616 ... e-log-data