python - montenegro - series de tiempo machine learning
Estimando el pequeño turno de tiempo entre dos series de tiempo (6)
Tengo dos series de tiempo, y sospecho que hay un turno de tiempo entre ellas, y quiero estimar este turno de tiempo.
Esta pregunta se ha planteado anteriormente en: Encuentre la diferencia de fase entre dos ondas (inarmónicas) y encuentre el cambio de tiempo entre dos formas de onda similares, pero en mi caso, el cambio de tiempo es menor que la resolución de los datos. por ejemplo, los datos están disponibles en resolución horaria, y el cambio de tiempo es de solo unos minutos (ver imagen).
La causa de esto es que el registrador de datos utilizado para medir una de las series tiene pocos minutos de cambio en su tiempo.
¿Algún algoritmo que pueda estimar este cambio, preferiblemente sin usar interpolación?
Este es un problema bastante interesante. Aquí hay un intento de una solución parcial usando transformadas de Fourier. Esto se basa en que los datos son moderadamente periódicos. No estoy seguro si funcionará con sus datos (donde los derivados en los puntos finales no parecen coincidir).
import numpy as np
X = np.linspace(0,2*np.pi,30) #some X values
def yvals(x):
return np.sin(x)+np.sin(2*x)+np.sin(3*x)
Y1 = yvals(X)
Y2 = yvals(X-0.1) #shifted y values
#fourier transform both series
FT1 = np.fft.fft(Y1)
FT2 = np.fft.fft(Y2)
#You can show that analyically, a phase shift in the coefficients leads to a
#multiplicative factor of `exp(-1.j * N * T_d)`
#can''t take the 0''th element because that''s a division by 0. Analytically,
#the division by 0 is OK by L''hopital''s<sp?> rule, but computers don''t know calculus :)
print np.log(FT2[1:]/FT1[1:])/(-1.j*np.arange(1,len(X)))
Una inspección rápida de la salida impresa muestra que las frecuencias con mayor potencia (N = 1, N = 2) dan estimaciones razonables, N = 3 también lo hace si se mira el valor absoluto (np.absolute), aunque I '' Estoy en una pérdida para explicar por qué sería eso.
Tal vez alguien más familiarizado con las matemáticas pueda sacarlo de aquí para dar una mejor respuesta ...
Uno de los enlaces que proporcionó tiene la idea correcta (de hecho, estoy haciendo más o menos lo mismo aquí)
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
a,b, N = 0, 10, 1000 #Boundaries, datapoints
shift = -3 #Shift, note 3/10 of L = b-a
x = np.linspace(a,b,N)
x1 = 1*x + shift
time = np.arange(1-N,N) #Theoritical definition, time is centered at 0
y1 = sum([np.sin(2*np.pi*i*x/b) for i in range(1,5)])
y2 = sum([np.sin(2*np.pi*i*x1/b) for i in range(1,5)])
#Really only helps with large irregular data, try it
# y1 -= y1.mean()
# y2 -= y2.mean()
# y1 /= y1.std()
# y2 /= y2.std()
cross_correlation = correlate(y1,y2)
shift_calculated = time[cross_correlation.argmax()] *1.0* b/N
y3 = sum([np.sin(2*np.pi*i*(x1-shift_calculated)/b) for i in range(1,5)])
print "Preset shift: ", shift, "/nCalculated shift: ", shift_calculated
plt.plot(x,y1)
plt.plot(x,y2)
plt.plot(x,y3)
plt.legend(("Regular", "Shifted", "Recovered"))
plt.savefig("SO_timeshift.png")
plt.show()
Esto tiene el siguiente resultado:
Preset shift: -3
Calculated shift: -2.99
Puede ser necesario verificar
Tenga en cuenta que la argmax () de la correlación muestra la posición de la alineación, tiene que ser escalada por la longitud de ba = 10-0 = 10 y N para obtener el valor real.
Comprobando el origen de la fuente de correlación no es del todo obvio cuál es el comportamiento de la función importada de sigtools. Para grandes conjuntos de datos, la correlación circular (a través de transformadas rápidas de Fourier) es mucho más rápida que el método directo. Sospecho que esto es lo que se implementa en sigtools, pero no puedo decirlo con certeza. Una búsqueda del archivo en mi carpeta python2.7 solo devolvió el archivo compilado C pyd.
Optimizar para la mejor solución
Para las restricciones dadas, es decir, que la solución está desplazada en fase por una pequeña cantidad menor que el método de muestreo, un simple algoritmo simplex cuesta abajo funciona bien. Modifiqué el problema de muestra de @mgilson para mostrar cómo hacer esto. Tenga en cuenta que esta solución es robusta, ya que puede manejar el ruido.
Función de error : puede haber más cosas óptimas para optimizar, pero esto funciona sorprendentemente bien:
np.sqrt((X1-X2+delta_x)**2+(Y1-Y2)**2).sum()
Es decir, minimice la distancia euclidiana entre las dos curvas ajustando solo el eje x (fase).
import numpy as np
def yvals(x):
return np.sin(x)+np.sin(2*x)+np.sin(3*x)
dx = .1
unknown_shift = .03 * np.random.random() * dx
X1 = np.arange(0,2*np.pi,dx) #some X values
X2 = X1 + unknown_shift
Y1 = yvals(X1)
Y2 = yvals(X2) # shifted Y
Y2 += .1*np.random.normal(size=X1.shape) # now with noise
def err_func(p):
return np.sqrt((X1-X2+p[0])**2+(Y1-Y2)**2).sum()
from scipy.optimize import fmin
p0 = [0,] # Inital guess of no shift
found_shift = fmin(err_func, p0)[0]
print "Unknown shift: ", unknown_shift
print "Found shift: ", found_shift
print "Percent error: ", abs((unknown_shift-found_shift)/unknown_shift)
Una ejecución de muestra da:
Optimization terminated successfully.
Current function value: 4.804268
Iterations: 6
Function evaluations: 12
Unknown shift: 0.00134765446268
Found shift: 0.001375
Percent error: -0.0202912082305
He utilizado con éxito (en el canal awgn) el enfoque de filtro combinado, que da la energía máxima m [n] en el índice n; luego ajustando un polinomio de segundo grado f (n) a m [n-1], m [n], m [n + 1] y encontrando el mínimo mediante el ajuste f ''(n) == 0.
La respuesta no es necesariamente absolutamente lineal, especialmente si la autocorrelación de la señal no desaparece en m [n-1], m [n + 1].
De hecho, es un problema interesante, pero todavía no hay una respuesta satisfactoria. Tratemos de cambiar eso ...
Usted dice que prefiere no utilizar la interpolación, pero, como entiendo de su comentario, lo que realmente quiere decir es que le gustaría evitar el muestreo superior a una resolución más alta. Una solución básica hace uso de un ajuste por mínimos cuadrados con una función de interpolación lineal, pero sin subir el muestreo a una resolución más alta:
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import leastsq
def yvals(x):
return np.sin(x)+np.sin(2*x)+np.sin(3*x)
dx = .1
X = np.arange(0,2*np.pi,dx)
Y = yvals(X)
unknown_shift = np.random.random() * dx
Y_shifted = yvals(X + unknown_shift)
def err_func(p):
return interp1d(X,Y)(X[1:-1]+p[0]) - Y_shifted[1:-1]
p0 = [0,] # Inital guess of no shift
found_shift = leastsq(err_func,p0)[0][0]
print "Unknown shift: ", unknown_shift
print "Found shift: ", found_shift
Una ejecución de muestra proporciona una solución bastante precisa:
Unknown shift: 0.0695701123582
Found shift: 0.0696105501967
Si uno incluye ruido en el Y desplazado:
Y_shifted += .1*np.random.normal(size=X.shape)
Uno obtiene resultados algo menos precisos:
Unknown shift: 0.0695701123582
Found shift: 0.0746643381744
La precisión en presencia de ruido mejora cuando hay más datos disponibles, por ejemplo con:
X = np.arange(0,200*np.pi,dx)
Un resultado típico es:
Unknown shift: 0.0695701123582
Found shift: 0.0698527939193
Este es un problema muy interesante. Originalmente, iba a sugerir una solución basada en la correlación cruzada similar a la de user948652. Sin embargo, de la descripción de su problema, hay dos problemas con esa solución:
- La resolución de los datos es mayor que el cambio de hora, y
- En algunos días, el valor predicho y los valores medidos tienen una muy baja correlación entre sí
Como resultado de estas dos cuestiones, creo que es probable que la aplicación directa de la solución de correlación cruzada aumente realmente su cambio de horario, particularmente en los días en que los valores predichos y medidos tienen una muy baja correlación entre sí.
En mi comentario anterior, le pregunté si tenía algún evento que ocurriera en ambas series de tiempo, y usted dijo que no. Sin embargo, según su dominio, creo que realmente tiene dos:
- amanecer
- Puesta de sol
Incluso si el resto de la señal está mal correlacionada, el amanecer y la puesta de sol deberían estar correlacionados de alguna manera, ya que monótonamente aumentarán / disminuirán a la línea de base de la noche. Así que aquí hay una solución potencial, basada en estos dos eventos, que debería minimizar la interpolación necesaria y no depender de la correlación cruzada de señales mal correlacionadas.
1. Encuentra el amanecer / ocaso aproximado
Esto debería ser lo suficientemente fácil, simplemente tome el primer y último punto de datos que sean más altos que la línea plana de la noche, y etiquete aquellos con el amanecer y el ocaso aproximados. Luego, me centraría en esa información, así como en los puntos inmediatamente a cada lado, es decir:
width=1
sunrise_index = get_sunrise()
sunset_index = get_sunset()
# set the data to zero, except for the sunrise/sunset events.
bitmap = zeros(data.shape)
bitmap[sunrise_index - width : sunrise_index + width] = 1
bitmap[sunset_index - width : sunset_index + width] = 1
sunrise_sunset = data * bitmap
Hay varias formas de implementar get_sunrise() y get_sunset() dependiendo de la cantidad de rigor que necesite en su análisis. Yo usaría numpy.diff , lo numpy.diff a un valor específico y tomaría el primer y último punto por encima de ese valor. También puede leer los datos de la noche a partir de una gran cantidad de archivos, calcular la media y la desviación estándar, y buscar los puntos de datos primero y último que excedan, digamos, 0.5 * st_dev de los datos de la noche. También podría hacer algún tipo de coincidencia de plantilla basada en clúster, en particular si las diferentes clases de días (es decir, soleado vs parcialmente nublado vs muy nublado) tienen eventos altamente estereotípicos de salida / puesta del sol.
2. Remuestrear datos
No creo que haya ninguna forma de resolver este problema sin alguna interpolación. Utilizaría remuestrear los datos a una frecuencia de muestreo más alta que el cambio. Si el cambio está en la escala de minutos, entonces la muestra ascendente a 1 minuto o 30 segundos.
num_samples = new_sample_rate * sunrise_sunset.shape[0]
sunrise_sunset = scipy.signal.resample(sunrise_sunset, num_samples)
Alternativamente, podríamos usar una spline cúbica para interpolar los datos (ver aquí ).
3. Convolución gaussiana
Como hay algo de interpolación, entonces no sabemos cómo se predijo exactamente el amanecer y el ocaso reales. Entonces, podemos convolucionar la señal con un gaussiano, para representar esta incertidumbre.
gaussian_window = scipy.signal.gaussian(M, std)
sunrise_sunset_g = scipy.signal.convolve(sunrise_sunset, gaussian_window)
4. Correlación cruzada
Use el método de correlación cruzada en la respuesta del usuario948652 para obtener el cambio de tiempo.
Hay muchas preguntas sin respuesta en este método que requerirían un examen y una experimentación con los datos para delimitar más específicamente, como cuál es el mejor método para identificar el amanecer / ocaso, qué tan amplia debe ser la ventana gaussiana, etc. Pero es cómo comenzaría a atacar el problema. ¡Buena suerte!