Как мне подобрать синусоидальную кривую к моим данным с помощью pylab и numpy?

Question

Для школьного проекта я пытаюсь показать, что экономики следуют относительно синусоидальной схеме роста. Помимо экономики, которая, по общему признанию, является хитрой, я строю симуляцию питона, чтобы показать, что даже если мы допустим некоторую степень случайности, мы все равно можем произвести что-то относительно синусоидальное. Я доволен своими данными, которые я создаю, но теперь я хотел бы найти способ получить синусоидальный график, который очень близко соответствует данным. Я знаю, что вы можете сделать полиномиальную подгонку, но вы можете сделать синусную подгонку?

Заранее благодарны за Вашу помощь. Дайте мне знать, есть ли какие-то части кода, которые вы хотите увидеть.

Answer 1

Вы можете использовать функцию оптимизации наименьших квадратов в scipy, чтобы подогнать любую произвольную функцию к другой. В случае подбора синусоидальной функции для подгонки используются три параметра: смещение ("a"), амплитуда ("b") и фаза ("c").

Пока вы даете разумное первое предположение о параметрах, оптимизация должна хорошо сходиться. К счастью для синусоидальной функции, первые оценки 2 из них просты: смещение можно оценить, взяв среднее значение данных и амплитуду через RMS (3* стандартное отклонение / кв.м. (2)).

Примечание: в качестве более позднего редактирования также была добавлена ​​подгонка частоты. Это не очень хорошо работает (может привести к очень плохим припадкам). Таким образом, используйте по своему усмотрению, я бы посоветовал не использовать подгонку по частоте, если ошибка по частоте не меньше, чем несколько процентов.

Это приводит к следующему коду:

import numpy as np
from scipy.optimize import leastsq
import pylab as plt

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
f = 1.15247 # Optional!! Advised not to use
data = 3.0*np.sin(f*t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise

guess_mean = np.mean(data)
guess_std = 3*np.std(data)/(2**0.5)/(2**0.5)
guess_phase = 0
guess_freq = 1
guess_amp = 1

# we'll use this to plot our first estimate. This might already be good enough for you
data_first_guess = guess_std*np.sin(t+guess_phase) + guess_mean

# Define the function to optimize, in this case, we want to minimize the difference
# between the actual data and our "guessed" parameters
optimize_func = lambda x: x[0]*np.sin(x[1]*t+x[2]) + x[3] - data
est_amp, est_freq, est_phase, est_mean = leastsq(optimize_func, [guess_amp, guess_freq, guess_phase, guess_mean])[0]

# recreate the fitted curve using the optimized parameters
data_fit = est_amp*np.sin(est_freq*t+est_phase) + est_mean

# recreate the fitted curve using the optimized parameters

fine_t = np.arange(0,max(t),0.1)
data_fit=est_amp*np.sin(est_freq*fine_t+est_phase)+est_mean

plt.plot(t, data, '.')
plt.plot(t, data_first_guess, label='first guess')
plt.plot(fine_t, data_fit, label='after fitting')
plt.legend()
plt.show()

Изменить: я предположил, что вы знаете количество периодов в синусоиде. Если вы этого не сделаете, это немного сложнее, чтобы соответствовать. Вы можете попытаться угадать количество периодов вручную и постараться оптимизировать его как шестой параметр.

Answer 2

Вот функция подгонки без параметров fit_sin() это не требует ручного предположения частоты:

import numpy, scipy.optimize

def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)}

Начальное предположение частоты задается пиковой частотой в частотной области с использованием БПФ. Результат подгонки почти идеален, если предположить, что существует только одна доминирующая частота (кроме пика нулевой частоты).

import pylab as plt

N, amp, omega, phase, offset, noise = 500, 1., 2., .5, 4., 3
#N, amp, omega, phase, offset, noise = 50, 1., .4, .5, 4., .2
#N, amp, omega, phase, offset, noise = 200, 1., 20, .5, 4., 1
tt = numpy.linspace(0, 10, N)
tt2 = numpy.linspace(0, 10, 10*N)
yy = amp*numpy.sin(omega*tt + phase) + offset
yynoise = yy + noise*(numpy.random.random(len(tt))-0.5)

res = fit_sin(tt, yynoise)
print( "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )

plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt, yynoise, "ok", label="y with noise")
plt.plot(tt2, res["fitfunc"](tt2), "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

Результат хорош даже при высоком уровне шума:

Амплитуда =1.00660540618, угловая частота =2.03370472482, фаза =0.360276844224, смещение =3.95747467506, макс. Ков.= +0,0122923578658

Answer 3

Более удобной для нас является функция Curvefit. Вот пример:

import numpy as np
from scipy.optimize import curve_fit
import pylab as plt

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
data = 3.0*np.sin(t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise

guess_freq = 1
guess_amplitude = 3*np.std(data)/(2**0.5)
guess_phase = 0
guess_offset = np.mean(data)

p0=[guess_freq, guess_amplitude,
    guess_phase, guess_offset]

# create the function we want to fit
def my_sin(x, freq, amplitude, phase, offset):
    return np.sin(x * freq + phase) * amplitude + offset

# now do the fit
fit = curve_fit(my_sin, t, data, p0=p0)

# we'll use this to plot our first estimate. This might already be good enough for you
data_first_guess = my_sin(t, *p0)

# recreate the fitted curve using the optimized parameters
data_fit = my_sin(t, *fit[0])

plt.plot(data, '.')
plt.plot(data_fit, label='after fitting')
plt.plot(data_first_guess, label='first guess')
plt.legend()
plt.show()

Answer 4

Текущие методы для подгонки синусоидальной кривой к заданному набору данных требуют первого угадывания параметров с последующим интерактивным процессом. Это проблема нелинейной регрессии.

Другой метод состоит в преобразовании нелинейной регрессии в линейную регрессию благодаря удобному интегральному уравнению. Тогда нет необходимости в первоначальном предположении и в итеративном процессе: подгонка получается напрямую.

В случае функции y = a + r*sin(w*x+phi) или же y=a+b*sin(w*x)+c*cos(w*x)см. стр. 35-36 статьи "Régression sinusoidale" опубликовано на Scribd

В случае функции y = a + p*x + r*sin(w*x+phi): стр. 49-51 главы "Смешанные линейные и синусоидальные регрессии".

В случае более сложных функций общий процесс объясняется в главе "Generalized sinusoidal regression" страницы 54-61, сопровождаемые числовым примером y = r*sin(w*x+phi)+(b/x)+c*ln(x)стр. 62-63

Answer 5

Все приведенные выше ответы основаны на подборе кривой, и в большинстве из них используется итеративный метод - все они работают очень хорошо, но я хотел добавить другой подход с использованием БПФ. Здесь мы преобразуем данные, устанавливаем все частоты, кроме пика, равными нулю, а затем выполняем обратное преобразование. Обратите внимание, что вы, вероятно, захотите удалить среднее значение данных (и устранение тренда) перед выполнением БПФ, а затем вы можете добавить их обратно после.

      import numpy as np
import pylab as plt

# fake data

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
f = 1.05 # Optional!! Advised not to use
data = 3.0*np.sin(f*t+0.001) + np.random.randn(N) # create artificial data with noise

# FFT...
mfft=np.fft.fft(data)
imax=np.argmax(np.absolute(mfft))
mask=np.zeros_like(mfft)
mask[[imax]]=1
mfft*=mask
fdata=np.fft.ifft(mfft)


plt.plot(t, data, '.')
plt.plot(t, fdata,'.', label='FFT')
plt.legend()
plt.show()

Answer 6

Если вы уже знаете частоту, вы можете выполнить линейную аппроксимацию, которая в вычислительном отношении более эффективна, чем методы нелинейной аппроксимации. Как отмечает @JJacquelin, никаких первоначальных предположений не требуется.

Явно тебе подойдетy=a+b*sin(x)+c*sin(x)к данным. Обратите внимание: это эквивалентноA*sin(x+phi), по тригонометрическим тождествам. Однако это выражается линейным образом в подходящих параметрах (но не в x,y). Таким образом, мы можем использовать линейную подгонку в Python.

Просто притворись, чтоx1 = sin(x)иx2 = cos(x)являются входными данными, используйте функцию линейной аппроксимацииy = a + b* x1 + c* x2

Для этого используйте

      from sklearn.linear_model import LinearRegression
reg = LinearRegression()

x= # your x values list
y = # your y values list
X = np.column_stack((np.sin(x), np.cos(x)))

reg.fit(X, y)

Вы можете получить доступ к параметрам подгонки:

      a = reg.intercept_
b = reg.coef_[0]
c = reg.coef_[1]

Answer 7

Следующий скрипт Python выполняет аппроксимацию синусоид методом наименьших квадратов, как описано в Bloomfield (2000), "Анализ Фурье временных рядов", Wiley. Ключевые шаги следующие:

1. Определите диапазон различных возможных частот.
2. Для каждой из частот, указанных в пункте 1 выше, оцените все параметры синусоиды (среднее значение, амплитуда и фаза) обычным методом наименьших квадратов (OLS).
3. Среди всех различных наборов параметров, оцененных в пункте 2 выше, выберите тот, который минимизирует остаточную сумму квадратов (RSS).

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
######################################################################################
# (1) generate the data
######################################################################################
omega=4.5  # frequency
theta0=0.5 # mean
theta1=1.5 # amplitude
phi=-0.25  # phase

n=1000 # number of observations
sigma=1.25 # error standard deviation
e=np.random.normal(0,sigma,n) # Gaussian error

t=np.linspace(1,n,n)/n # time index
y=theta0+theta1*np.cos(2*np.pi*(omega*t+phi))+e # time series
######################################################################################
# (2) estimate the parameters
######################################################################################
# define the range of different possible frequencies
f=np.linspace(1,12,100)

# create a data frame to store the estimated parameters and the associated
# residual sum of squares for each of the considered frequencies
coefs=pd.DataFrame(data=np.zeros((len(f),5)), columns=["omega","theta0","theta1","phi","RSS"])

for i in range(len(f)): # iterate across the different considered frequencies

    x1=np.cos(2*np.pi*f[i]*t) # define the first regressor
    x2=np.sin(2*np.pi*f[i]*t) # define the second regressor

    X=sm.add_constant(np.transpose(np.vstack((x1,x2)))) # add the intercept
    reg=sm.OLS(y,X).fit() # fit the regression by OLS

    A=reg.params[1] # estimated coefficient of first regressor
    B=reg.params[2] # estimated coefficient of second regressor

    # estimated mean
    mean=reg.params[0]

    # estimated amplitude
    amplitude=np.sqrt(A**2+B**2)

    # estimated phase
    if A>0:          phase=np.arctan(-B/A)/(2*np.pi)
    if A<0 and B>0:  phase=(np.arctan(-B/A)-np.pi)/(2*np.pi)
    if A<0 and B<=0: phase=(np.arctan(-B/A)+np.pi)/(2*np.pi)
    if A==0 and B>0: phase=-1/4
    if A==0 and B<0: phase=1/4

    # fitted sinusoid
    s=mean+amplitude*np.cos(2*np.pi*(f[i]*t+phase))

    # residual sum of squares
    RSS=np.sum((y-s)**2)

    # save the estimation results
    coefs["omega"][i]=f[i]
    coefs["theta0"][i]=mean
    coefs["theta1"][i]=amplitude
    coefs["phi"][i]=phase
    coefs["RSS"][i]=RSS

    del x1, x2, X, reg, A, B, mean, amplitude, phase, s, RSS
######################################################################################
# (3) analyze the results
######################################################################################
# extract the set of parameters that minimizes the residual sum of squares
coefs=coefs.loc[coefs["RSS"]==coefs["RSS"].min(),]

# calculate the fitted sinusoid
s=coefs["theta0"].values+coefs["theta1"].values*np.cos(2*np.pi*(coefs["omega"].values*t+coefs["phi"].values))

# plot the fitted sinusoid
plt.plot(y,color="black",linewidth=1,label="actual")
plt.plot(s,color="lightgreen",linewidth=3,label="fitted")
plt.ylim(ymax=np.max(y)*1.3)
plt.legend(loc=1)
plt.show()