Pasos claros para calcular coherencia entre dos series de tiempo

Question

Originalmente publicado en este stackoverflow.com y la elimina y se trasladó aquí

Mi pregunta es similar a La similitud de los dos discreta de fourier tranforms (específicamente la respuesta seleccionada). También he recogido alguna información útil a partir de este R ayuda hilo.

Estoy atascado en la aplicación real de este proceso, sin embargo.

Tengo dos series de Fourier; mi objetivo final es calcular la suma de la ponderación de la coherencia. (En la actualidad tengo millones de series correspondientes a los píxeles, pero, vamos a llamarlo dos por ahora).

    TS1 <- c(5051.29, 5355.31, 5602.18, 5784.4, 5896.45, 5934.9, 5898.6, 
5788.64, 5608.37, 5363.27, 5060.78, 4710.09, 4321.86, 3907.89, 
3480.75, 3053.42, 2638.89, 2249.75, 1897.82, 1593.81, 1346.93, 
1164.71, 1052.67, 1014.21, 1050.52, 1160.47, 1340.74, 1585.84, 
1888.33, 2239.02, 2627.25, 3041.22, 3468.36, 3895.69, 4310.22, 
4699.36, 5051.29)
    TS2 <- c(4192.83, 4532.62, 4836.41, 5094.96, 5300.41, 5446.53, 5528.88, 
5544.95, 5494.25, 5378.33, 5200.71, 4966.78, 4683.65, 4359.93, 
4005.45, 3630.98, 3247.9, 2867.85, 2502.38, 2162.59, 1858.8, 
1600.25, 1394.8, 1248.67, 1166.33, 1150.26, 1200.95, 1316.87, 
1494.5, 1728.43, 2011.55, 2335.28, 2689.76, 3064.23, 3447.31, 
3827.36, 4192.83)

De los enlaces arriba veo que necesito para calcular la coherencia de la cruz espectro en cada frecuencia, y el peso de las frecuencias que tienen una mayor densidad espectral de potencia.

Sé que se puede extraer de la densidad espectral de

spectrum(data.frame(TS1, TS2))$spec

Sé que se puede extraer de coherencia en el uso de:

spectrum(data.frame(TS1, TS2))$coh

Puedo ver que el valor de la coherencia extraído es siempre igual a 1. Parece que la "solución" a este problema es el lag de la serie de tiempo de alguna manera.

No estoy seguro de cómo 1) calcular el óptimo lag, sobre todo porque me gustaría que sean comparables a través de todos los cálculos de coherencia, y 2) cómo especificar lag en el espectro(). spec.pgram toma "abarca", pero los que parecen tener más que ver con el suavizado de compensación.

He intentado un máximo de ccf función definida aquí, pero en muchos casos (no en el ejemplo que me has proporcionado, tendrás que confiar en mí) el "óptimo lag" devuelto es 0, lo cual tiene sentido porque la ccf es 1, pero no es realmente lo que estoy buscando.

Debo ser el uso de acf() o ccf() a calcular manualmente la coherencia? ¿Qué pasos (R funciones sería realmente útil) necesito realizar para llegar desde donde yo estoy, para tener una sola ponderado coherencia valor?

Answer 1

Los siguientes son los temas que voy a resumir de manera sucinta en la secuencia, ya que cada uno de ellos son un problema interesante en sí mismo, y luego seguir con mis soluciones:

i) Cómo unlag dos series de tiempo para perder el mínimo de datos debido a la zaga, y al mismo tiempo para maximizar la correlación cruzada o coherencia espectral.

ii) ¿cómo podría esto ser hecho cuando se trata de lidiar con multivariante de series de tiempo. Como es mayor que el de los dos, es como resolver un cubo de rubiks de forma óptima, mientras que quedando varios de series de tiempo en el conjunto de datos con el fin de preservar un resumió medida de (i).

iii) Cómo calcular un promedio ponderado de la coherencia?

Voy a omitir ii) y la respuesta i) y iii) debido a razones de conflicto de interés y colaboración en la investigación/confidencial acuerdos sobre mi trabajo en ii). Entiendo que iii) pasa a ser la cuestión central para empezar, seguido por i).

Mi solución para la iii)

Pre-Procesamiento Pasos: Decir que había dos series de tiempo t1 y t2. Determinar un CCF maximización de lag, decir optLag y unlag t1 y t2 como unlaggedTS = lag(t2,optLag) y realizar una series de tiempo de la unión del conjunto de datos como unionTS = ts.union(t1,unlaggedTS). Siguiente que, como el de series de tiempo unlagging produciría un par de NA, como bajar de entradas por desfasado, usted necesita para establecer estos valores a cero en unionTS según sea necesario, en el sentido de que nada se observó en el NA puntos.

Paso 2: Ahora, obtener un abinded de series de tiempo, con el lapso de los parámetros según sea necesario.

bindTSPair= abind(((spectrum(unionTS,span=c(16,16),plot=F)$spec)[,1]),((spectrum(unionTS,span=c(16,16),plot=F)$spec)[,2])) 

Siguiente, a normalizar los datos en bindTSPair como decir:

normalizedTS=round(data.Normalization(c(bindTSPair),type="n4"),6)

Luego se obtienen los siguientes dos mitades:

TSHalf1=normalizedTS[1:(length(bindTSPair)/2)]

TSHalf2=normalizedTS[((length(bindTSPair)/2)+1):(length(bindTSPair))]

y, finalmente, obtener el 'peso' de la versión de coherencia espectral de la siguiente manera:

 weightedCoherence=sum((10*(TSHalf1 + TSHalf2)) * (spectrum(unionTS,span=c(10,10),plot=F)$coh))

Usted puede cambiar la ponderación requerida o el uso de 'antes' info.

Mi solución para i) para obtener un CCF maximización de lag:

 for(i in 2:ncol(myTS))
{
    for(j in 2:ncol(myTS))
    {
r=ccf(myTS[,i],myTS[,j],lag.max=1000,plot=TRUE)
tp=sort((r$acf),index.return=TRUE)$ix
tn=sort((r$acf),index.return=TRUE,decreasing=T)$ix
cn=r$acf[tn[length(tn)]]
    cp=r$acf[tp[length(tp)]]
    if(abs(cp) > abs(cn))
    {
    temp=r$lag[tp[length(tp)]]
    	 tcc=cp
    	}
    		if(abs(cn) > abs(cp))
    		{
    		temp=r$lag[tn[length(tn)]]
    		tp=tn
    		 tcc=cn
    		}
    if(abs(cn) > abs(cp))
    {
    temp=r$lag[tp[length(tp)]]
}
show(i)
if(length(tp) != 0)
{

mlag[i,j]=temp
mcc[i,j]=tcc
}
    }

}

mlag=mlag[2:ncol(myTS),2:ncol(myTS)]
mlag=data.frame(mlag)

Answer 2

a menos que su tarea es poner en práctica la coherencia de cálculo, yo no lo haría, porque no tiene que ser un paquete. por ejemplo, veo coh paquete, lo he usado yo. esta función parece hacer todo que usted necesita para una determinada frecuencia f, por lo que puede pasar de una llamada a través del rango de frecuencia.

En matlab existe una función llamada mscohere, que produce todo un gráfico como este, es el cuadrado de la magnitud de la coherencia

Yo apolologise para el uso de Matlab ejemplo, pero yo no hice ningún análisis espectral en R mí todavía. Los conceptos son los mismos, aunque.

Answer 3

Uso `` para extraer los valores para el punto de máxima coherencia coherencia y lag:

La `` función está diseñada para trabajar con objetos de ts. Esta es una manera simple de hacer el cálculo Praneeth describe aquí.