2012-02-27 15 views
12

Ho un 2396x34 double matrix chiamato y cui ciascuna riga (2396) rappresenta una situazione separata costituita da 34 segmenti di tempo consecutivi.Correlazione di Pearson ponderata?

Ho anche un numeric[34] denominato x che rappresenta una singola situazione di 34 segmenti di tempo consecutivi.

Attualmente sto calcolando la correlazione tra ciascuna riga y e x simili:

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

cosa devo ora è quello di sostituire la funzione cor nell'istruzione sopra con un ponderato correlazione . Il peso vettore xy.wt è lungo 34 elementi in modo che sia possibile assegnare un peso diverso a ciascuno dei 34 segmenti di tempo consecutivi.

Ho trovato la funzione Weighted Covariance Matrixcov.wt e ho pensato che se prima i dati di scale dovessero funzionare come la funzione cor. In effetti è possibile specificare che la funzione restituisca anche una matrice di correlazione. Sfortunatamente non sembra che io possa usarlo allo stesso modo perché non posso fornire separatamente le mie due variabili (x e).

Qualcuno sa di un modo per ottenere una correlazione ponderata nel modo che ho descritto senza sacrificare la velocità?

Edit: Forse qualche funzione matematica potrebbe essere applicato a y prima della funzione di cor al fine di ottenere gli stessi risultati che sto cercando. Forse se moltiplico ciascun elemento per xy.wt/sum(xy.wt)?

Modifica n. 2 Ho trovato un'altra funzione corr nel pacchetto 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. 

Anche questo non è quello che mi serve ma è più vicino.

Modifica # 3 Ecco il codice per generare il tipo di dati che sto lavorando con:

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 

risposta

4

si può tornare alla definizione della correlazione.

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) 
+0

Eccellente! Questo è stato. Grazie ancora! Pensavo che le funzioni scritte in R sarebbero state sostanzialmente più lente di quelle incorporate in R ... ma immagino di no? –

22

Purtroppo la risposta accettata è sbagliato quando y è una matrice di più di una riga. L'errore è nella linea

vy <- rowSums(w * y * y) 

vogliamo moltiplicare le colonne di y da w, ma questo moltiplicherà le righe dagli elementi di w, riciclato come necessario.Così

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

è corretto, perché in questo caso la moltiplicazione viene eseguita elemento saggio, che equivale a colonna-saggio moltiplicazione qui, ma

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

dà una risposta sbagliata a causa della row- saggia moltiplicazione

Possiamo correggere la funzione come segue

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) 
} 

e verificare i risultati contro quelli prodotti dalla corr dal pacchetto 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 

che di per sé dà un altro modo in cui questo problema potrebbe essere risolto.

+0

@vincentzoonekynd Forse dovresti dare un'occhiata a questo e commentare? – Andrie

+0

Nella mia risposta c'è un bug (volevo cancellarlo, ma non è possibile cancellare le risposte accettate). Di solito mi aspetto un avvertimento quando moltiplico gli oggetti con dimensioni errate, ma in questo caso non ce n'erano ... –

+0

Ho pensato che sarebbe stato meglio aggiungere un commento e modificare la tua risposta, mi dispiace per quello. Almeno adesso il bug è stato contrassegnato e tu hai ancora il merito di fare la maggior parte del lavoro! –

2

Ecco una generalizzazione per calcolare la correlazione di Pearson ponderata tra due matrici (invece di un vettore e una matrice, come nella domanda iniziale):

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))) 
} 

Utilizzando l'esempio di cui sopra e la funzione di correlazione da Heather , possiamo verificare che:

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

In termini di chiamare sintassi, questo assomiglia alla non ponderata 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