Jak dopasować krzywą sinusoidalną do moich danych za pomocą pylab i numpy?

Question

W przypadku projektu szkolnego staram się pokazać, że gospodarki wykazują względnie sinusoidalny wzorzec wzrostu. Poza ekonomią tego, co jest z pewnością sprytne, buduję symulację Pythona, aby pokazać, że nawet jeśli pozwolimy na pewien stopień losowości, nadal możemy wytworzyć coś stosunkowo sinusoidalnego. Cieszę się z moich danych, które produkuję, ale teraz chciałbym znaleźć jakiś sposób, by uzyskać wykres sinusoidalny, który dokładnie pasuje do danych. Wiem, że możesz zrobić dopasowanie wielomianowe, ale czy możesz zrobić sinus fit?

Dzięki za pomoc z góry. Daj mi znać, jeśli są jakieś fragmenty kodu, który chcesz zobaczyć.

Answer 1

Możesz użyć funkcji least-square optimization w scipy, aby dopasować dowolną funkcję do innej. W przypadku dopasowania funkcji sin, 3 parametry do dopasowania to przesunięcie ("a"), amplituda ("b") i faza ("c").

Dopóki podajesz rozsądne pierwsze przypuszczenie o parametrach, optymalizacja powinna zbiegać się dobrze. Na szczęście dla funkcji sinusowej, pierwsze oszacowania 2 z nich są łatwe: przesunięcie można oszacować, biorąc średnią danych oraz amplitudę za pomocą RMS (3 * odchylenie standardowe/sqrt (2)).

Prowadzi to do następującego kodu:

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) 
data = 3.0*np.sin(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) 
guess_phase = 0 

# 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(t+x[1]) + x[2] - data 
est_std, est_phase, est_mean = leastsq(optimize_func, [guess_std, guess_phase, guess_mean])[0] 

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

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

EDIT: Zakłada się, że wiesz liczbę okresów w sinusa. Jeśli nie, to nieco trudniejsze dopasowanie. Możesz spróbować zgadnąć liczbę okresów przez ręczne kreślenie i spróbuj zoptymalizować go jako swój czwarty parametr.

Answer 2

Obecne metody dopasowania krzywej sinus do danego zestawu danych wymagają najpierw odgadnięcia parametrów, a następnie procesu mieszania. Jest to problem regresji nieliniowej.

Inna metoda polega na przekształceniu regresji nieliniowej w regresję liniową dzięki dogodnemu równaniu całemu. Wtedy nie ma potrzeby wstępnego odgadywania i nie ma potrzeby powtarzania procesu: dopasowanie jest uzyskiwane bezpośrednio.

W przypadku funkcji y = a + r*sin(w*x+phi) lub y=a+b*sin(w*x)+c*cos(w*x) opisano na stronach 35-36 w dokumencie "Régression sinusoidale" opublikowanym Scribd

W przypadku funkcji y = a + p*x + r*sin(w*x+phi): strony 49-51 rozdziału „Mixed liniowej regresji sinusoidalnym” .

W przypadku bardziej skomplikowanych funkcji, ogólny proces jest wyjaśnione w rozdziale "Generalized sinusoidal regression" stronach 54-61, a następnie przez przykładzie liczbowym y = r*sin(w*x+phi)+(b/x)+c*ln(x), strony 62-63

Answer 3

Bardziej przyjazny dla użytkownika jest dla nas curvefit funkcja. Oto przykład:

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

tutaj jest funkcją montażu parametr wolne fit_sin() które nie wymagają ręcznego przypuszczenie częstotliwości:

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)} 

Początkowy przypuszczenie częstotliwości podanej przez częstotliwość piku w dziedzinie częstotliwości za pomocą FFT. Wynik dopasowania jest prawie idealny, zakładając, że istnieje tylko jedna dominująca częstotliwość (inna niż szczytowa częstotliwość zerowa).

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() 

Wynik jest dobry, nawet przy wysokim hałasem:

Amplituda = 1,00660540618, kątowe freq = 2,03370472482, faza = ,360276844224, offset = 3.95747467506, max.. Cov. = 0,0122923578658