2010-10-27 6 views
3

J'ai 13 matrices de différentes dimensions que je voudrais utiliser dans les corrélations matricielles par paires avec une fonction personnalisée (qui calcule le coefficient Rv). La fonction prend deux arguments (matrix1, matrix2) et produit un scalaire (fondamentalement une valeur r multivariée). Je voudrais exécuter la fonction sur toutes les paires de matrices possibles (donc 78 corrélations au total) et produire une matrice 13 par 13 des valeurs Rv résultantes avec les noms des 13 matrices dans les lignes et les colonnes. J'ai pensé essayer de faire cela en plaçant les matrices dans une liste et en utilisant un double pour passer en revue les éléments de la liste, mais cela semble très complexe. J'ai donné un exemple réduit avec des données factices ci-dessous. Quelqu'un at-il des suggestions sur la façon dont cela pourrait être abordé? Merci d'avance.Corrélations matricielles par paires dans R - comment parcourir toutes les paires?

# Rv function 
Rv <- function(M1, M2) { 
    tr <- function(x) sum(diag(x)) 
    psd <- function(x) x %*% t(x) 
    AA <- psd(M1) 
    BB <- psd(M2) 
    num <- tr(AA %*% BB) 
    den <- sqrt(tr(AA %*% AA) * tr(BB %*% BB)) 
    Rv <- num/den 
    list(Rv=Rv, "Rv^2"=Rv^2) 
} 

# data in separate matricies 
matrix1 <- matrix(rnorm(100), 10, 10) 
matrix2 <- matrix(rnorm(100), 10, 10) 
# ... etc. up to matrix 13 

# or, in a list 
matrix1 <- list(matrix(rnorm(100), 10, 10)) 
rep(matrix1, 13) # note, the matrices are identical in this example 

# call Rv function 
Rv1 <- Rv(matrix1, matrix2) 
Rv1$Rv^2 

# loop through all 78 combinations? 
# store results in 13 by 13 matrix with matrix rownames and colnames? 
+0

Votre fonction 'psd' est déjà implémentée en tant que' tcrossprod' (et est plus efficace). – Marek

Répondre

3

Ce que j'utilisé dans le passé est expand.grid() suivie apply(). Voici un exemple plus simple utilisant seulement 1: 3 plutôt que 1:13.

R> work <- expand.grid(1:3,1:3) 
R> work 
    Var1 Var2 
1 1 1 
2 2 1 
3 3 1 
4 1 2 
5 2 2 
6 3 2 
7 1 3 
8 2 3 
9 3 3 
R> apply(work, 1, function(z) prod(z)) 
[1] 1 2 3 2 4 6 3 6 9 
R> 

Vous souhaitez évidemment une autre fonction de travail.

Questions connexes