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

Question

У меня есть 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)?

Edit # 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.

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

Edit # 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

Answer 1

К сожалению, принятый ответ неверен, когда 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

, что само по себе дает другой способ решения этой проблемы.

Answer 2

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

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

Answer 3

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

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's 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)