Estimación del pequeño cambio de tiempo entre dos series de tiempo

Tengo dos series de tiempo, y sospecho que hay un cambio de tiempo entre ellas, y quiero estimar este cambio de tiempo.

Esta pregunta se ha formulado anteriormente en: Encuentre la diferencia de fase entre dos ondas (inarmónicas) y encuentre el cambio en el tiempo entre dos formas de onda similares, pero en mi caso, el cambio en el tiempo es más pequeño que la resolución de los datos. por ejemplo, los datos están disponibles a una resolución por hora, 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 unos pocos minutos de cambio en su tiempo.

¿Algún algoritmo por ahí que pueda estimar este cambio, preferiblemente sin usar interpolación?

Predicción de la irradiación solar y medición de la irradiación solar.

Este es un problema bastante interesante. Aquí hay un bash de una solución parcial utilizando transformadas de Fourier. Esto se basa en que los datos son moderadamente periódicos. No estoy seguro de 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 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 la mayor potencia (N = 1, N = 2) proporcionan estimaciones razonables, N = 3 también está bien si observa el valor absoluto (np.absolute), aunque I ‘ No puedo explicar por qué sería eso.

Tal vez alguien más familiarizado con las matemáticas pueda tomarlo desde aquí para dar una mejor respuesta …

Uno de los enlaces que proporcionaste tiene la idea correcta (de hecho, aquí estoy haciendo casi lo mismo)

 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 = ba 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 la siguiente salida:

 Preset shift: -3 Calculated shift: -2.99 

introduzca la descripción de la imagen aquí

Puede ser necesario comprobar

  1. Scipy Correlate
  2. Analaysis Delay De Tiempo

Tenga en cuenta que el argmax () de la correlación muestra la posición de la alineación, debe ser escalado por la longitud de ba = 10-0 = 10 y N para obtener el valor real.

Al verificar la fuente de la fuente correlacionada, no es del todo obvio cómo se comporta la función importada de sigtools. Para grandes conjuntos de datos, la correlación circular (a través de las 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 asegurarlo. Una búsqueda del archivo en mi carpeta python2.7 solo devolvió el archivo comstackdo en C pyd.

Este es un problema muy interesante. Originalmente, iba a sugerir una solución basada en correlación cruzada similar a la de user948652. Sin embargo, a partir de la descripción de su problema, hay dos problemas con esa solución:

  1. La resolución de los datos es mayor que el cambio de tiempo, y
  2. En algunos días, el valor predicho y los valores medidos tienen una correlación muy baja entre sí

Como resultado de estos dos problemas, creo que la aplicación directa de la solución de correlación cruzada probablemente boostá su turno de tiempo, especialmente en los días en que los valores predichos y medidos tienen una correlación muy baja entre sí.

En mi comentario anterior, le pregunté si tuvo algún evento que ocurriera en ambas series de tiempo y dijo que no. Sin embargo, según tu dominio, creo que en realidad tienes dos:

  1. amanecer
  2. Puesta de sol

Incluso si el rest de la señal está mal correlacionada, la salida y la puesta del sol deberían estar algo correlacionadas, ya que boostán / disminuirán monótonamente a la línea de base nocturna. Entonces, 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 la salida del sol aproximada / puesta del sol

Esto debería ser bastante fácil, simplemente tome los primeros y últimos puntos de datos que son más altos que la línea plana de la noche, y etiquételos como el amanecer y el atardecer aproximados. Luego, me concentraría en esos datos, 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. numpy.diff , lo numpy.diff en un valor de umbral específico, y numpy.diff el primer y último punto por encima de ese valor. También puede leer los datos nocturnos desde una gran cantidad de archivos, calcular la media y la desviación estándar, y buscar el primer y último punto de datos que exceda, por ejemplo, 0.5 * st_dev de los datos nocturnos. También podría hacer algún tipo de coincidencia de plantillas basadas en clústeres, en particular si diferentes clases de día (es decir, soleado frente a nublado o muy nublado) tienen eventos de amanecer / atardecer altamente estereotipados.

2. Volver a muestrear datos

No creo que haya ninguna manera de resolver este problema sin alguna interpolación. Yo usaría remuestrear los datos a una tasa de muestreo más alta que el cambio. Si el cambio está en la escala de minutos, aumente 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. La convolución gaussiana

Dado que hay cierta interpolación, entonces no sabemos cómo se predijeron precisamente el amanecer y el atardecer. Entonces, podemos convertir 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

Utilice el método de correlación cruzada en la respuesta del usuario948652 para obtener el cambio de hora.

Hay muchas preguntas sin respuesta en este método que requieren un examen y experimentación con los datos para concretar más específicamente, como cuál es el mejor método para identificar el amanecer / anochecer, qué tan amplia debe ser la ventana gaussiana, etc. Cómo empezaría a atacar el problema. ¡Buena suerte!

De hecho, problema interesante, pero ninguna respuesta satisfactoria todavía. Tratemos de cambiar eso …

Usted dice que prefiere no utilizar la interpolación, pero, según entiendo por su comentario, lo que realmente quiere decir es que le gustaría evitar el muestreo a una resolución más alta. Una solución básica utiliza un ajuste de mínimos cuadrados con una función de interpolación lineal, pero sin muestrear 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 da una solución bastante precisa:

 Unknown shift: 0.0695701123582 Found shift: 0.0696105501967 

Si uno incluye ruido en la Y desplazada:

 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 

He utilizado con éxito (en un canal AWGN) el enfoque de filtro adaptado, que proporciona la energía pico m [n] en el índice n; luego ajuste un polinomio de segundo grado f (n) a m [n-1], m [n], m [n + 1] y encuentre el mínimo estableciendo 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].

Optimizar para la mejor solución.

Para las restricciones dadas, a saber, que la solución se desplaza en fase por una pequeña cantidad menor que el método de muestreo, un simple algoritmo simplex de bajada funciona bien. He modificado el problema de muestra de @mgilson para mostrar cómo hacerlo. Tenga en cuenta que esta solución es robusta, ya que puede manejar el ruido.

Función de error : puede haber cosas más ó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