2012-02-02 3 views
8

У меня есть длинные цифровые данные временных рядов приблизительно 200 000 строк (давайте назовем это Z).R: Эффективное размещение сегментов временных рядов с максимальной взаимной корреляцией с входным сегментом?

В цикле я подмножество х (около 30) последовательных строк из Z в то время, и рассматривать их в качестве точки запроса д.

Я хочу, чтобы определить местонахождение в пределах Zу (~ 300) наиболее коррелировали сегментов временных рядов длины х (наиболее коррелирует с ц).

Что такое эффективный способ достижения этого?

ответ

5

ниже код находит 300 сегментов вы ищете и работает через 8 секунд на моем не слишком мощном ноутбуке Windows, поэтому он должен быть достаточно быстрым для ваших целей.

Во-первых, он строит матрицу 30 на-199971 (Zmat), столбцы которой содержат все сегменты временных интервалов длиной 30 ", которые вы хотите изучить. Один вызов cor(), работающий на векторе q и матрице Zmat, затем вычисляет все требуемые коэффициенты корреляции. Наконец, результирующий вектор исследуется для идентификации 300 последовательностей, имеющих самые высокие коэффициенты корреляции.

# Simulate data 
nZ <- 200000 
nq <- 30 
Z <- rnorm(nZ) 
q <- seq_len(nq) 

# From Z, construct a 30 by 199971 matrix, in which each column is a 
# "time series segment". Column 1 contains observations 1:30, column 2 
# contains observations 2:31, and so on through the end of the series. 
Zmat <- sapply(seq_len(nZ - nq + 1), 
       FUN = function(X) Z[seq(from = X, length.out = nq)]) 

# Calculate the correlation of q with every column/"time series segment. 
Cors <- cor(q, Zmat) 

# Extract the starting position of the 300 most highly correlated segments  
ids <- order(Cors, decreasing=TRUE)[1:300] 

# Maybe try something like the following to confirm that you have 
# selected the most highly correlated segments. 
hist(Cors, breaks=100) 
hist(Cors[ids], col="red", add=TRUE) 
+0

Отличное решение. '' '' '' '' '' '' '' '' и, возможно, 'q' должен быть назначен' Z [seq_len (qlen)] '. В соответствии с терминологией Майка' qlen' = 'x'. – jbaums

+0

Спасибо для того, чтобы поймать это.Я теперь пересмотрел код, изменив имена переменных и добавив еще несколько комментариев/объяснений. Не уверен, что 'q' обязательно будет первым 30 элементами' Z', поэтому я не изменил это бит. –

3

Наивный решение действительно очень медленно (по крайней мере, несколько минут - я не достаточно терпеливы):

library(zoo) 
n <- 2e5 
k <- 30 
z <- rnorm(n) 
x <- rnorm(k) # We do not use the fact that x is a part of z 
rollapply(z, k, function(u) cor(u,x), align="left") 

Вы можете вычислить корреляцию вручную, с первых минут и comoments, но все еще занимает несколько минут.

y <- zoo(rnorm(n), 1:n) 
x <- rnorm(k) 
exy <- exx <- eyy <- ex <- ey <- zoo(rep(0,n), 1:n) 
for(i in 1:k) { 
    cat(i, "\n") 
    exy <- exy + lag(y,i-1) * x[i] 
    ey <- ey + lag(y,i-1) 
    eyy <- eyy + lag(y,i-1)^2 
    ex <- ex + x[i] # Constant time series 
    exx <- exx + x[i]^2 # Constant time series 
} 
exy <- exy/k 
ex <- ex/k 
ey <- ey/k 
exx <- exx/k 
eyy <- eyy/k 
covxy <- exy - ex * ey 
vx <- exx - ex^2 
vy <- eyy - ey^2 
corxy <- covxy/sqrt(vx * vy) 

После того, как у вас есть временной ряд корреляций, легко извлечь положение верхнего 300.

i <- order(corxy, decreasing=TRUE)[1:300] 
corxy[i]