Взвешенная корреляция Пирсона?

У меня есть 2396x34 double matrix с именем y, где каждая строка (2396) представляет отдельную ситуацию, состоящую из 34 последовательных временных сегментов.

У меня также есть numeric[34] с именем x, который представляет собой единую ситуацию из 34 последовательных временных сегментов.

В настоящее время я вычисляю корреляцию между каждой строкой в y и x следующим образом:

crs[,2] <- cor(t(y),x)

Теперь мне нужно заменить функцию cor в приведенном выше выражении на взвешенную корреляцию. Весовой вектор xy.wt составляет 34 элемента, поэтому для каждого из 34 последовательных временных сегментов может быть назначен другой вес.

Я нашел функцию Weighted Covariance Matrix cov.wt и подумал, что если я сначала scale данные должны работать так же, как cor. Фактически вы можете указать, чтобы функция возвращала корреляционную матрицу. К сожалению, похоже, что я не могу использовать его таким же образом, потому что я не могу предоставить две переменные (x и y) отдельно.

Кто-нибудь знает, как я могу получить взвешенную корреляцию в том, как я описал, не жертвуя большой скоростью?

Изменить: возможно, какая-то математическая функция может быть применена к y до функции cor, чтобы получить те же результаты, что и я. Может быть, если я умножу каждый элемент на xy.wt/sum(xy.wt)?

Изменить # 2 Я нашел еще одну функцию corr в пакете boot.

corr(d, w = rep(1, nrow(d))/nrow(d))

d   
A matrix with two columns corresponding to the two variables whose correlation we wish to calculate.

w   
A vector of weights to be applied to each pair of observations. The default is equal weights for each pair. Normalization takes place within the function so sum(w) need not equal 1.

Это тоже не то, что мне нужно, но оно ближе.

Изменить # 3 Вот несколько кодов для создания типа данных, с которыми я работаю:

x<-cumsum(rnorm(34))
y<- t(sapply(1:2396,function(u) cumsum(rnorm(34))))
xy.wt<-1/(34:1)

crs<-cor(t(y),x) #this works but I want to use xy.wt as weight

Ответ 1

Вы можете вернуться к определению корреляции.

f <- function( x, y, w = rep(1,length(x))) {
  stopifnot( length(x) == dim(y)[2] )
  w <- w / sum(w)
  # Center x and y, using the weighted means
  x <- x - sum(x*w)
  y <- y - apply( t(y) * w, 2, sum )
  # Compute the variance
  vx <- sum( w * x * x )
  vy <- rowSums( w * y * y ) # Incorrect: see Heather remark, in the other answer
  # Compute the covariance
  vxy <- colSums( t(y) * x * w )
  # Compute the correlation
  vxy / sqrt(vx * vy)
}
f(x,y)[1]
cor(x,y[1,]) # Identical
f(x, y, xy.wt)

Ответ 2

К сожалению, принятый ответ неверен, если y - это матрица из нескольких строк. Ошибка находится в строке

vy <- rowSums( w * y * y )

Мы хотим умножить столбцы y на w, но это умножит строки на элементы w, при необходимости переработанные. Таким образом,

> f(x, y[1, , drop = FALSE], xy.wt)
[1] 0.103021

является правильным, так как в этом случае умножение выполняется по-разному, что эквивалентно умножению по столбцам здесь, но

> f(x, y, xy.wt)[1]
[1] 0.05463575

дает неправильный ответ из-за умножения строки.

Мы можем исправить функцию следующим образом:

f2 <- function( x, y, w = rep(1,length(x))) {
  stopifnot(length(x) == dim(y)[2] )
  w <- w / sum(w)
  # Center x and y, using the weighted means
  x <- x - sum(x * w)
  ty <- t(y - colSums(t(y) * w))
  # Compute the variance
  vx <- sum(w * x * x)
  vy <- colSums(w * ty * ty)
  # Compute the covariance
  vxy <- colSums(ty * x * w)
  # Compute the correlation
  vxy / sqrt(vx * vy)
}

и проверьте результаты, полученные от corr, из пакета boot:

> res1 <- f2(x, y, xy.wt)
> res2 <- sapply(1:nrow(y), 
+                function(i, x, y, w) corr(cbind(x, y[i,]), w = w),
+                x = x, y = y, w = xy.wt)
> all.equal(res1, res2)
[1] TRUE

который сам по себе дает другой способ решить эту проблему.

Ответ 3

Вот обобщение для вычисления взвешенной корреляции Пирсона между двумя матрицами (вместо вектора и матрицы, как в исходном вопросе):

matrix.corr <- function (a, b, w = rep(1, nrow(a))/nrow(a)) 
{
    # normalize weights
    w <- w / sum(w)

    # center matrices
    a <- sweep(a, 2, colSums(a * w))
    b <- sweep(b, 2, colSums(b * w))

    # compute weighted correlation
    t(w*a) %*% b / sqrt( colSums(w * a**2) %*% t(colSums(w * b**2)) )
}

Используя приведенный выше пример и функцию корреляции от Хизер, мы можем проверить это:

> sum(matrix.corr(as.matrix(x, nrow=34),t(y),xy.wt) - f2(x,y,xy.wt))
[1] 1.537507e-15

В терминах синтаксиса вызова это напоминает невзвешенный cor:

> a <- matrix( c(1,2,3,1,3,2), nrow=3)
> b <- matrix( c(2,3,1,1,7,3,5,2,8,1,10,12), nrow=3)
> matrix.corr(a,b)
     [,1]      [,2] [,3]      [,4]
[1,] -0.5 0.3273268  0.5 0.9386522
[2,]  0.5 0.9819805 -0.5 0.7679882
> cor(a, b)
     [,1]      [,2] [,3]      [,4]
[1,] -0.5 0.3273268  0.5 0.9386522
[2,]  0.5 0.9819805 -0.5 0.7679882