Wie passe ich eine Sinuskurve an meine Daten mit pylab und numpy?

Question

Für ein Schulprojekt versuche ich zu zeigen, dass Volkswirtschaften einem relativ sinusförmigen Wachstumsmuster folgen. Jenseits der Ökonomie, die zugegebenermaßen zweifelhaft sind, baue ich eine Python-Simulation, um zu zeigen, dass selbst wenn wir ein gewisses Maß an Zufälligkeit zulassen, wir immer noch etwas relativ sinusförmiges produzieren können. Ich bin glücklich mit meinen Daten, die ich produziere, aber jetzt möchte ich einen Weg finden, ein Sinus-Diagramm zu erhalten, das den Daten ziemlich genau entspricht. Ich weiß, dass du eine polynomische Anpassung machen kannst, aber kannst du sine fit machen?

Vielen Dank für Ihre Hilfe im Voraus. Lassen Sie mich wissen, ob Teile des Codes angezeigt werden sollen.

Answer 1

Sie können die Funktion least-square optimization in scipy verwenden, um eine beliebige Funktion an eine andere anzupassen. Im Falle der Anpassung einer Sinusfunktion sind die 3 Parameter, die passen, der Offset ('a'), die Amplitude ('b') und die Phase ('c').

Solange Sie eine vernünftige erste Schätzung der Parameter liefern, sollte die Optimierung gut konvergieren. Zum Glück für eine Sinusfunktion sind erste Schätzungen von 2 einfach: Der Offset kann geschätzt werden, indem der Mittelwert der Daten genommen wird und die Amplitude über den RMS (3 * Standardabweichung/sqrt (2)).

Dies führt zu dem folgenden Code:

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: Ich nahm an, dass Sie die Anzahl der Perioden in den Sinuswellen kennen. Wenn Sie es nicht tun, ist es etwas schwieriger zu passen. Sie können die Anzahl der Perioden durch manuelles Plotten erraten und versuchen, sie als vierten Parameter zu optimieren.

Answer 2

Die aktuellen Methoden zum Anpassen einer Sinuskurve an einen gegebenen Datensatz erfordern eine erste Schätzung der Parameter, gefolgt von einem interativen Prozess. Dies ist ein nichtlineares Regressionsproblem.

Eine andere Methode besteht darin, die nichtlineare Regression durch eine einfache Integralgleichung in eine lineare Regression zu transformieren. Dann besteht keine Notwendigkeit für eine anfängliche Schätzung und keine Notwendigkeit für einen iterativen Prozess: Die Anpassung wird direkt erhalten.

Bei der Funktion y = a + r*sin(w*x+phi) oder y=a+b*sin(w*x)+c*cos(w*x), Seiten 35-36 des Papiers sehen "Régression sinusoidale" veröffentlicht am Scribd

Bei der Funktion y = a + p*x + r*sin(w*x+phi): Seiten 49-51 des Kapitels „Mixed linear und sinus Regressionen“ .

Bei komplizierteren Funktionen wird der allgemeine Prozess "Generalized sinusoidal regression" Seiten 54-61, gefolgt von einem numerischen Beispiel y = r*sin(w*x+phi)+(b/x)+c*ln(x), Seiten 62-63

Answer 3

Weiteren anwenderfreundliche wir im Kapitel ist die Funktion CurveFit. Hier ein Beispiel:

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

Hier ist eine parameterfreie Anpassungsfunktion fit_sin(), die keine manuelle Vermutung der Frequenz erfordert:

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

Die Anfangsfrequenz Vermutung durch die Spitzenfrequenz im Frequenzbereich gegeben ist mit FFT. Das Anpassungsergebnis ist nahezu perfekt, wenn angenommen wird, dass nur eine dominante Frequenz existiert (anders als die Nullfrequenzspitze).

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

Das Ergebnis ist gut, auch mit hohem Lärm:

Amplitude = 1,00660540618, Angular freq = 2,03370472482, Phase = ,360276844224, offset = 3,95747467506, Max.. Cov. = 0,0122923578658