python - ¿Cómo puedo ajustar una curva sinusoidal a mis datos con pylab y numpy?

sine economics

Question

Para un proyecto escolar, estoy tratando de mostrar que las economías siguen un patrón de crecimiento relativamente sinusoidal. Más allá de su economía, que ciertamente son poco fiables, estoy construyendo una simulación de python para mostrar que incluso cuando permitimos que se logre cierto grado de aleatoriedad, todavía podemos producir algo relativamente sinusoidal. Estoy contento con los datos que estoy produciendo, pero ahora me gustaría encontrar alguna forma de obtener un gráfico sinusoidal que coincida bastante con los datos. Sé que puedes hacer un ajuste polinomial, pero ¿puedes hacer un ajuste sinusoidal?

Gracias por su ayuda por adelantado. Déjame saber si hay alguna parte del código que quieres ver.

Answer 1

Aquí hay una función de ajuste sin parámetros fit_sin() que no requiere una estimación manual de la frecuencia:

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

El cálculo de frecuencia inicial viene dado por la frecuencia pico en el dominio de la frecuencia utilizando FFT. El resultado de ajuste es casi perfecto, suponiendo que solo hay una frecuencia dominante (distinta del pico de frecuencia cero).

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

El resultado es bueno incluso con alto ruido:

Amplitud = 1.00660540618, Frecuencia angular = 2.03370472482, fase = 0.360276844224, desplazamiento = 3.95747467506, Máx. Cov. = 0.0122923578658

Answer 2

Los métodos actuales para ajustar una curva de pecado a un conjunto de datos dado requieren una primera conjetura de los parámetros, seguidos de un proceso interactivo. Este es un problema de regresión no lineal.

Un método diferente consiste en transformar la regresión no lineal en una regresión lineal gracias a una ecuación integral conveniente. Entonces, no hay necesidad de una estimación inicial ni de un proceso iterativo: la adaptación se obtiene directamente.

En el caso de la función y = a + r*sin(w*x+phi) o y=a+b*sin(w*x)+c*cos(w*x) , consulte las páginas 35-36 del documento "Régression sinusoidale" publicado en Scribd

En el caso de la función y = a + p*x + r*sin(w*x+phi) : páginas 49-51 del capítulo "Regresiones lineales mixtas y sinusoidales".

En el caso de funciones más complicadas, el proceso general se explica en el capítulo "Generalized sinusoidal regression" páginas 54-61, seguido de un ejemplo numérico y = r*sin(w*x+phi)+(b/x)+c*ln(x) , páginas 62-63

Answer 3

Más amigable para nosotros es la función curvefit. Aquí un ejemplo:

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

Puede usar la función de optimización de mínimos cuadrados en scipy para ajustar cualquier función arbitraria a otra. En el caso de ajustar una función sin, los 3 parámetros que se ajustan son el desplazamiento (''a''), la amplitud (''b'') y la fase (''c'').

Siempre que proporcione una primera aproximación razonable de los parámetros, la optimización debería converger bien. Afortunadamente para una función sinusoidal, las primeras estimaciones de 2 de estos son fáciles: el desplazamiento puede estimarse tomando la media de los datos y la amplitud a través de el RMS (3 * desviación estándar / sqrt (2)).

Nota: como una edición posterior, también se ha agregado el ajuste de frecuencia. Esto no funciona muy bien (puede llevar a ajustes extremadamente pobres). Por lo tanto, use a su discreción, mi consejo sería no usar el ajuste de frecuencia a menos que el error de frecuencia sea menor que un pequeño porcentaje.

Esto conduce al siguiente código:

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

Edición: asumí que conoces la cantidad de períodos en la onda sinusoidal. Si no lo haces, es un poco más difícil de ajustar. Puede probar y adivinar el número de períodos mediante el trazado manual y tratar de optimizarlo como su sexto parámetro.