2012-02-27 30 views
12

yと名付けられた2396x34 double matrixを有し、各行(2396)は、34の連続する時間セグメントからなる別個の状況を表す。重み付きピアソン相関?

numeric[34]xという名前で、34の連続した時間セグメントの単一の状況を表しています。現在、私はこのようyxの各行の間の相関を計算しています

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

私が今必要なものは加重相関と上記のステートメントでcor機能を交換することです。ウエイトベクトルxy.wtは34エレメント長であり、34の連続する時間セグメントのそれぞれに異なる重みを割り当てることができる。

Weighted Covariance Matrixファンクションcov.wtが見つかりました。最初にscaleのデータを入力すると、corのように機能するはずです。実際には、相関行列も返す関数を指定することができます。残念ながら、2つの変数(xy)を別々に供給することができないため、同じ方法で使用できるようには見えません。

多くのスピードを犠牲にすることなく記述した方法で加重相関を得る方法は誰にも分かりますか?

編集:おそらくいくつかの数学関数は、私が探している同じ結果を得るために前cor関数にyにも適用することができる。たぶん私は各要素をxy.wt/sum(xy.wt)で掛けていますか?

編集#2bootパッケージに別の機能corrが見つかりました。

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. 

これは私が必要とするものではありませんが、それはより近くです。あなたが戻って相関の定義に行くことができます

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 

答えて

4

:ここ

編集#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) 
+0

優秀!それがそれでした。再度、感謝します! Rで書かれた関数は、Rに組み込まれた関数よりもかなり遅くなると思っていましたが、そうは思いませんか? –

22

残念ながら、yが複数の行のマトリックスである場合、受け入れられた答えは間違っています。エラーは、ライン我々はwによってyの列を掛けたい

vy <- rowSums(w * y * y) 

であるが、これは、必要に応じてリサイクル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) 
} 

を次のように機能を修正し、bootパッケージからcorrによって生成されるものに対する結果を確認することができます:

> 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 

自体に問題があることができることを別の方法を提供します解決される。ここ

+0

@vincentzoonekyndおそらくあなたはこれを見てコメントする必要がありますか? – Andrie

+0

本当に私の答えにバグがあります(私はそれを削除したかったが、受け入れられた回答を削除することはできません)。私は通常、間違った次元でオブジェクトを掛けるときに警告を期待しますが、この場合は何もありませんでした... –

+0

その後、コメントを追加して、あなたの答えを編集させていただきました。少なくともバグは今すぐ上がっていますが、あなたはまだ仕事の大部分を行っていると信じています! –

2

は、(代わりに、ベクトル及び行列の、元の質問のような)は、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))) 
} 

ヘザーから上記の実施例との相関関数を使用して、我々はそれを確認することができます。

> 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 
012:構文を呼び出すという点で

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

が、これはcor重み付けされていないに似ています

関連する問題