J'ai besoin de calculer la distance de Mahalanobis échantillon dans R entre chaque paire d'observations dans une matrice de covariables. J'ai besoin d'une solution efficace, c'est-à-dire que seules distances sont calculées et de préférence implémentées dans C / RCpp / Fortran etc. Je suppose que , la matrice de covariance de population, est inconnue et utilise la covariance d'échantillon matrice à sa place.
Je m'intéresse particulièrement à cette question car il ne semble pas y avoir de méthode "consensus" pour calculer les distances de Mahalanobis par paires dans R, c'est-à-dire qu'elle n'est pas implémentée dans la dist
fonction ni dans la cluster::daisy
fonction. La mahalanobis
fonction ne calcule pas les distances par paire sans travail supplémentaire du programmeur.
Cela a déjà été demandé ici la distance de Mahalanobis Pairwise en R , mais les solutions semblent incorrectes.
Voici une méthode correcte mais terriblement inefficace (puisque distances sont calculées):
set.seed(0)
x0 <- MASS::mvrnorm(33,1:10,diag(c(seq(1,1/2,l=10)),10))
dM = as.dist(apply(x0, 1, function(i) mahalanobis(x0, i, cov = cov(x0))))
C'est assez facile pour me coder en C, mais je pense que quelque chose de ce basique devrait avoir une solution préexistante. Est-ce qu'il y a un?
Il existe d'autres solutions qui échouent: HDMD::pairwise.mahalanobis()
calcule distances, lorsque seules distances uniques sont requises. compositions::MahalanobisDist()
semble prometteur, mais je ne veux pas que ma fonction provienne d'un package qui en dépend rgl
, ce qui limite sérieusement la capacité des autres à exécuter mon code. À moins que cette implémentation ne soit parfaite, je préfère écrire la mienne. Quelqu'un a-t-il de l'expérience avec cette fonction?
la source
Réponses:
A partir de la solution "succint" d'Ahfoss, j'ai utilisé la décomposition de Cholesky à la place du SVD.
Cela devrait être plus rapide, car la résolution directe d'un système triangulaire est plus rapide que la multiplication matricielle dense avec la covariance inverse ( voir ici ). Voici les références avec les solutions d'Ahfoss et de Whuber dans plusieurs contextes:
Donc Cholesky semble être uniformément plus rapide.
la source
La formule standard pour la distance au carré de Mahalanobis entre deux points de données est
où est un vecteur p × 1 correspondant à l'observation i . En règle générale, la matrice de covariance est estimée à partir des données observées. Sans compter l'inversion de matrice, cette opération nécessite p 2 + p multiplications et p 2 + 2 p additions, chacune répétéexi p×1 i p2+p p2+2p fois.n(n−1)/2
Considérez la dérivation suivante:
où . Notez quexTiΣ-1qi=Σ−12xi . Cela repose sur le fait queΣ-1xTiΣ−12=(Σ−12xi)T=qTi est symétrique, ce qui tient au fait que pour toute matrice diagonalisable symétriqueA=PEΣ- 12 ,A = PEPT
Si nous laissons , et notons que Σ - 1 est symétrique, nous voyons queA = Σ- 1 Σ- 1 doit également être symétrique. SiXest lamatricen×pdes observations etQest lamatricen×ptelle que laithligne deQsoitqΣ- 12 X n × p Q n × p jet h Q , alors Q peut être exprimé succinctement par X Σ - 1qje Q . Ceci et les résultats précédents impliquent queXΣ−12
les seules opérations calculées n ( n - 1 ) / 2 fois sont p multiplications et 2 p additions (par opposition aux p 2 + p multiplications et p 2 + 2 p
la source
pair.diff()
passe et également donner un exemple numérique avec des impressions de chaque étape de votre fonction? Merci.Essayons l'évidence. De
il s'ensuit que nous pouvons calculer le vecteur
en temps et la matriceO(p2)
en temps , très probablement en utilisant des opérations de tableau rapide intégrées (parallélisables), puis formez la solution commeO(pn2+p2n)
où est le produit extérieur par rapport à + : ( a ⊕ b ) i j = a i + b j .⊕ + (a⊕b)ij=ai+bj.
UneΣ=Var(X) h
R
implémentation est succinctement parallèle à la formulation mathématique (et suppose, avec elle, que est en fait inversible avec l'inverse écrit h ici):Notez, pour des raisons de compatibilité avec les autres solutions, que seuls les éléments hors diagonale uniques sont renvoyés, plutôt que la matrice de distance carrée entière (symétrique, zéro sur la diagonale). Les diagrammes de dispersion montrent que ses résultats sont en accord avec ceux de
fastPwMahal
.En C ou C ++, la RAM peut être réutilisée et calculée à la volée, évitant tout besoin de stockage intermédiaire de u ⊕ u .u⊕u u⊕u
Les études de synchronisation avec allant de 33 à 5000 et p allant de 10 à 100 indiquent que cette mise en œuvre est 1,5 à 5 fois plus rapide que dans cette plage. L'amélioration s'améliore à mesure que p et n augmentent. Par conséquent, nous pouvons nous attendre à être supérieurs pour des p plus petits . Le seuil de rentabilité se produit autour de p = 7 pour nn 33 5000 p 10 100 1.5 5 p n p p=7 n≥100 . Le fait que les mêmes avantages de calcul de cette solution simple s'appliquent dans d'autres implémentations peut être une question de la façon dont ils tirent parti des opérations de tableau vectorisé.
fastPwMahal
fastPwMahal
la source
apply
etouter
... à l'exception de l'éclatementRcpp
.R
n'y a rien à y gagner.Si vous souhaitez calculer l' échantillon de distance de Mahalanobis, vous pouvez exploiter quelques astuces algébriques. Ils conduisent tous à calculer des distances euclidiennes par paire, supposons donc que nous pouvons utiliserX n×p p O(np)
dist()
pour cela. Soit la matrice de données n × p , que nous supposons être centrée de sorte que ses colonnes aient une moyenne de 0 et avoir un rang p de sorte que la matrice de covariance de l'échantillon soit non singulière. (Le centrage nécessite des opérations O ( n p ) .) Alors la matrice de covariance échantillon est S = X T X / n .Les distances Mahalanobis par paire de de l'échantillon sont les mêmes que les distances euclidiennes par paire de X L pour toute matrice L satisfaisant L L T = S - 1 , par exemple la racine carrée ou le facteur de Cholesky. Cela découle d'une algèbre linéaire et conduit à un algorithme nécessitant le calcul de S , S - 1 et une décomposition de Cholesky. La complexité la plus défavorable est O ( n p 2 + p 3 ) .X
Plus profondément, ces distances sont liées aux distances entre les principales composantes de l'échantillon de . Soit X = U D V T désignent la SVD de X . Ensuite , S = V D 2 V T / n et S - 1 / 2 = V D - 1 V T n 1 / 2 . Donc , X S - 1 / deux = U V T n 1X X=UDVT X
Voici une implémentation R de la deuxième méthode que je ne peux pas tester sur l'iPad que j'utilise pour écrire cette réponse.
la source
Il s'agit d'une solution beaucoup plus succincte. Il est toujours basé sur la dérivation impliquant la matrice de covariance de racine carrée inverse (voir mon autre réponse à cette question), mais utilise uniquement la base R et le package de statistiques. Il semble être légèrement plus rapide (environ 10% plus rapide dans certains benchmarks que j'ai exécutés). Notez qu'il renvoie la distance Mahalanobis, par opposition à la distance Maha au carré.
Cette fonction nécessite une matrice de covariance inverse et ne renvoie pas d'objet distance - mais je soupçonne que cette version allégée de la fonction sera plus généralement utile pour empiler les utilisateurs d'échange.
la source
SQRT
par la décomposition de Choleskychol(invCovMat)
.Si vous n'utilisez que les fonctionnalités de Fortran77 dans l'interface, votre sous-programme est toujours suffisamment portable pour les autres.
la source
Il existe un moyen très simple de le faire en utilisant le package R "biotools". Dans ce cas, vous obtiendrez une matrice de mahalanobis à distance carrée.
la source
Ceci est le code étendu de mon ancienne réponse déplacé ici à partir d'un autre fil .
Je fais depuis longtemps le calcul d'une matrice symétrique carrée de distances de Mahalanobis par paire dans SPSS via un approche de matrice de chapeau en utilisant la résolution d'un système d'équations linéaires (car c'est plus rapide que l'inversion de la matrice de covariance).
Je ne suis pas un utilisateur R donc j'ai juste essayé de reproduire cette recette @ahfoss ici dans SPSS avec "ma" recette, sur une donnée de 1000 cas par 400 variables, et j'ai trouvé mon chemin beaucoup plus rapidement.
Ainsi, centrez les colonnes de la matrice de données, calculez la matrice chapeau, multipliez par (n-1) et effectuez l'opération à l'opposé du double-centrage. Vous obtenez la matrice des distances de Mahalanobis au carré.
Dans nos paramètres, la matrice "double-centrée" est spécifiquement le chapeau matrice (multipliée par n-1), et non les produits scalaires euclidiens, et la matrice de distance carrée résultante est donc la matrice de distance Mahalanobis carrée, pas la matrice de distance euclidienne carrée.
En notation matricielle: Soit la diagonale de H ( n - 1 ) , vecteur de colonne. Propager la colonne dans la matrice carrée: ; ensuiteH H(n−1) D2mahal=H+H′−2H(n−1) .
H= {H,H,...}
Le code dans SPSS et sonde de vitesse est ci-dessous.
Ce premier code correspond à la fonction @ahfoss
fastPwMahal
de la réponse citée . C'est l'équivalent mathématiquement. Mais je calcule la matrice symétrique complète des distances (via les opérations matricielles) tandis que @ahfoss a calculé un triangle de la matrice symétrique (élément par élément).Ce qui suit est ma modification pour le rendre plus rapide:
Enfin, "l'approche matricielle chapeau". Pour la vitesse, je calcule la matrice du chapeau (les données doivent être centrées en premier)X(X′X)−1X′ (X′X)−1X′
solve(X'X,X')
la source
La formule que vous avez publiée ne calcule pas ce que vous pensez que vous calculez (une statistique U).
Dans le code que j'ai publié, j'utilise
cov(x1)
comme matrice de mise à l'échelle (c'est la variance des différences par paire des données). Vous utilisezcov(x0)
(il s'agit de la matrice de covariance de vos données d'origine). Je pense que c'est une erreur de votre part. L'intérêt de l'utilisation des différences par paires est qu'il vous libère de l'hypothèse que la distribution multivariée de vos données est symétrique autour d'un centre de symétrie (ou d'avoir à estimer ce centre de symétrie d'ailleurs, car ilcrossprod(x1)
est proportionnel àcov(x1)
). De toute évidence, en utilisant,cov(x0)
vous perdez cela.Ceci est bien expliqué dans le document auquel j'ai lié ma réponse originale.
la source
Matteo Fasiolo
et (je suppose)whuber
dans ce fil. Le vôtre est différent. Je serais intéressé à comprendre ce que vous calculez, mais il est clairement différent de la distance de Mahalanobis telle qu'elle est généralement définie.cov(x0)