Вопрос:
Для школьного проекта я пытаюсь показать, что экономики следуют относительно синусоидальной модели роста. Помимо экономической теории, которые, по общему признанию, изворотливые, я строю моделирование на питоне, чтобы показать, что даже когда мы допустим некоторую степень случайности, мы все равно можем создать нечто относительно синусоидальное. Я доволен своими данными, которые я создаю, но теперь Id хотел найти какой-то способ получить синус-граф, который очень близко соответствует данным. Я знаю, что вы можете делать полиномиальную посадку, но можете ли вы сделать синус?
Спасибо за вашу помощь заранее. Дайте мне знать, есть ли какие-либо части кода, которые вы хотите увидеть.
Лучший ответ:
Вы можете использовать функцию оптимизации наименьших квадратов в scipy для подгонки любой произвольной функции к другой. В случае подбора функции sin, 3 подходящих параметра – это смещение (“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()
Изменение: я предположил, что вы знаете количество периодов в синусоиде. Если вы этого не сделаете, это несколько сложнее, чтобы соответствовать. Вы можете попытаться угадать количество периодов вручную и постараться оптимизировать его как шестой параметр.
Ответ №1
Вот функция установки без параметров 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, Angular частота = 2,03370472482, фаза = 0,360276844224, смещение = 3,95747467506, макс. Cov. = +0,0122923578658
Ответ №2
Другим дружественным для нас способом является функция 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() Ответ №3
Существующие методы, чтобы соответствовать кривой sin для данного набора данных, требуют первой догадки параметров, за которыми следует интерактивный процесс. Это проблема нелинейной регрессии.
Другой метод заключается в преобразовании нелинейной регрессии в линейную регрессию благодаря удобному интегральному уравнению. Тогда нет необходимости в первоначальной догадки и нет необходимости в итеративном процессе: фитинг получен непосредственно.
В случае функции 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
Ответ №4
Следующий скрипт Python выполняет подгонку синусоид методом наименьших квадратов, как описано в Bloomfield (2000), “Анализ Фурье временных рядов”, Wiley. Ключевыми шагами являются следующие:
- Определите диапазон различных возможных частот.
- Для каждой из частот, указанных в точке 1 выше, оцените все параметры синусоиды (среднее значение, амплитуду и фазу) по обычным наименьшим квадратам (OLS).
- Среди всех различных наборов параметров, оцененных в точке 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()
