Пакет Wavemulcor - функция корреляции wave.multiple.cross.correlation - замена имеет нулевую длину

Я хотел бы использовать пакет wavemulcor и, в частности, функцию wave.multiple.cross.correlation для выполнения множественной взаимной корреляции вейвлетов на моих данных.

Я следую примеру в соответствии с руководством по пакету, но вместо этого использую свои данные. Функция работает с данными примера, но выдает ошибку, когда я пробую ее с моими. Ошибка относится к «замене имеет нулевую длину», но я не уверен, что именно это означает.

Я погуглил ошибку, но есть много примеров одной и той же проблемы для разных функций, и, как правило, все они имеют какое-то отношение к циклам в коде.

Затем я погуглил, как решить проблему, и прочитал об отладке. Я пробовал отлаживать код, но не могу понять, где он ломается, я все еще на ранних стадиях обучения программированию. Я думаю, что это может быть этот раздел кода в функции wave.multiple.cross.correlation, который вызывает проблему:

xy.cor.vec <- matrix(unlist(xy.cor), l, dd)
    xy.mulcor <- matrix(0, l, 2 * lm + 1)
    YmaxR <- vector("numeric", l)
    for (i in 1:l) {
        r <- xy.cor.vec[i, ]
        P <- diag(d)/2
        P[lower.tri(P)] <- r
        P <- P + t(P)
        Pidiag <- diag(solve(P))
        if (is.null(ymaxr)) {
            YmaxR[i] <- Pimax <- which.max(Pidiag)
        }

Есть ли другие способы определить, почему это не работает?

Фактическая ошибка выглядит следующим образом:

> Lst <- wave.multiple.cross.correlation(xx, lag.max = NULL, ymaxr = NULL)
Error in YmaxR[i] <- Pimax <- which.max(Pidiag) : 
  replacement has length zero

Я попытался отследить все переменные в исходном коде, но я просто не могу понять это.

Это код, который я пытаюсь использовать, и для полноты картины здесь - это ссылка на dput() xx, который представляет собой список переменных, которые я хочу использовать, см. Код ниже для получения подробной информации о xx

library(wavemulcor)
library(readxl)

rm(list = ls()) # clear objects
graphics.off() # close graphics windows

RNC20_30Hourly <- read_excel(RNC20_30Hourly.xlsx")

RNC20_30Hourly <- RNC20_30Hourly[-1]

RNC20_30TS <- ts(RNC20_30Hourly, start = 1, frequency = 23)

wf <- "d4"
J <- 6
lmax <- 36

n <- nrow(RNC20_30TS)

CK0158U09A3.modwt <- modwt(RNC20_30TS[,"CK0158U09A3"], wf, J)
CK0158U09A3.modwt.bw <- brick.wall(CK0158U09A3.modwt, wf)

CK0158U21A1.modwt <- modwt(RNC20_30TS[,"CK0158U21A1"], wf, J)
CK0158U21A1.modwt.bw <- brick.wall(CK0158U21A1.modwt, wf)

CK0158U21A2.modwt <- modwt(RNC20_30TS[,"CK0158U21A2"], wf, J)
CK0158U21A2.modwt.bw <- brick.wall(CK0158U21A2.modwt, wf)

CK0158U21A3.modwt <- modwt(RNC20_30TS[,"CK0158U21A3"], wf, J)
CK0158U21A3.modwt.bw <- brick.wall(CK0158U21A3.modwt, wf)

xx <- list(CK0158U09A3.modwt.bw, 
           CK0158U21A1.modwt.bw, 
           CK0158U21A2.modwt.bw,
           CK0158U21A3.modwt.bw)

Lst <- wave.multiple.cross.correlation(xx, lag.max = 13, ymaxr = NULL)

CK0158.RTWP.cross.cor <- as.matrix(Lst$xy.mulcor[1:J,])
YmaxR <- Lst$YmaxR

cell.names <- c("CK0158U09A3", "CK0158U21A1", "CK0158U21A2", "CK0158U21A3")
rownames(CK0158.RTWP.cross.cor)<-rownames(CK0158.RTWP.cross.cor,
                                          do.NULL = FALSE, prefix = "Level ")
lags <- length(-lmax:lmax)

lower.ci <- tanh(atanh(CK0158.RTWP.cross.cor) - qnorm(0.975) /
                   sqrt(matrix(trunc(n/2^(1:J)), nrow=J, ncol=lags)- 3))
upper.ci <- tanh(atanh(CK0158.RTWP.cross.cor) + qnorm(0.975) /
                   sqrt(matrix(trunc(n/2^(1:J)), nrow=J, ncol=lags)- 3))

par(mfrow=c(3,2), las=1, pty="m", mar=c(2,3,1,0)+.1, oma=c(1.2,1.2,0,0))
for(i in J:1) {
  matplot((1:(2*lmax+1)),CK0158.RTWP.cross.cor[i,], type="l", lty=1, ylim=c(-1,1),
          xaxt="n", xlab="", ylab="", main=rownames(CK0158.RTWP.cross.cor)[[i]][1])
  if(i<3) {axis(side=1, at=seq(1, 2*lmax+1, by=12),
                labels=seq(-lmax, lmax, by=12))}
  #axis(side=2, at=c(-.2, 0, .5, 1))
  lines(lower.ci[i,], lty=1, col=2) ##Add Connected Line Segments to a Plot
  lines(upper.ci[i,], lty=1, col=2)
  abline(h=0,v=lmax+1) ##Add Straight horiz and vert Lines to a Plot
  text(1,1, labels=cell.names[YmaxR[i]], adj=0.25, cex=.8)
}
par(las=0)
mtext('Lag (hours)', side=1, outer=TRUE, adj=0.5)
mtext('Wavelet Multiple Cross-Correlation', side=2, outer=TRUE, adj=0.5)

Будем очень признательны за любую помощь в решении/устранении этой проблемы.


r time-series cross-correlation wavelet-transform
person TheGoat    schedule 19.07.2017    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Просмотрев код гребешком с мелкими зубьями, я обнаружил ошибку в своем коде. Согласно мануалу использование следующее:

wave.multiple.cross.correlation(xx, lag.max = NULL, ymaxr = NULL)

однако в приведенном примере автор использует другую переменную с именем lmax, которая указывает используемую задержку.

Lst <- wave.multiple.cross.correlation(xx, lmax)

Как вы можете видеть из моего примера выше, я указал два разных аргумента для задержки: lmax = 36 и lag.max = 13.

Ну, это стоило мне всего 100 репутации!

person TheGoat    schedule 22.07.2017

arrow_upward
1
arrow_downward

Кажется, что, учитывая длину вашего временного ряда, вы установили максимальный уровень вейвлета J слишком высоким (строка 14 вашего кода):

14 J <- 6

Вместо этого попробуйте J ‹- 4 и запустите

Lst ‹- wave.multiple.cross.correlation.prueba(xx, lag.max = lmax, ymaxr = NULL)

где lmax — это максимальная задержка, которую вы ранее установили на 36 (или любое другое разумное значение по вашему выбору).

Как правило, рекомендуемое значение для максимального уровня вейвлета составляет J = trunc(log2(T))-3 (на всякий случай), где T = длина временного ряда.

Фактически, сообщение об ошибке было следствием вашего исходного списка xx, где уровни d6 и s6 заполнены NaN. Обратите внимание, что это не часть пакета wavemulcor, а предыдущие вычисления вейвлет-преобразования перед входом в выбранную вами функцию пакета wavemulcor.

Надеюсь это поможет.

person user9897491    schedule 05.06.2018