Paarweise Mahalanobis-Entfernungen

18

Ich muss den Mahalanobis-Abstand in R zwischen jedem Beobachtungspaar in einer n×p Matrix von Kovariaten berechnen. Ich benötige eine effiziente Lösung, dh es werden nur Abstände berechnet und vorzugsweise in C / RCpp / Fortran usw. implementiert. Ich gehe davon aus, dass , die Populationskovarianzmatrix, unbekannt ist, und verwende die Stichprobe Kovarianzmatrix an seiner Stelle.Σn(n-1)/2Σ

Diese Frage interessiert mich besonders, da es anscheinend keine "Konsensus" -Methode zur paarweisen Berechnung von Mahalanobis-Abständen in R gibt, dh sie ist weder in der distFunktion noch in der cluster::daisyFunktion implementiert . Die mahalanobisFunktion berechnet keine paarweisen Abstände ohne zusätzliche Arbeit vom Programmierer.

Dies wurde hier bereits paarweise nach Mahalanobis Entfernung in R gefragt , aber die Lösungen dort scheinen falsch zu sein.

Hier ist eine korrekte, aber furchtbar ineffiziente Methode (da Entfernungen berechnet werden):n×n

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

Das ist einfach genug, um mich in C zu programmieren, aber ich bin der Meinung, dass dieses Basic eine bereits vorhandene Lösung haben sollte. Ist dort eines?

Es gibt andere Lösungen, die zu kurz kommen: HDMD::pairwise.mahalanobis()Berechnet Abstände, wenn nur eindeutige Abstände erforderlich sind. scheint vielversprechend zu sein, aber ich möchte nicht, dass meine Funktion von einem Paket kommt , das davon abhängt , was die Fähigkeit anderer stark einschränkt, meinen Code auszuführen. Sofern diese Implementierung nicht perfekt ist, schreibe ich lieber meine eigene. Hat jemand Erfahrung mit dieser Funktion?n×nn(n-1)/2compositions::MahalanobisDist()rgl

ahfoss
quelle
Herzlich willkommen. Können Sie die beiden Matrizen der Distanz in Ihrer Frage drucken? Und was ist für Sie "ineffizient"?
TTNPHNS
1
Verwenden Sie nur die Beispiel-Kovarianzmatrix? Wenn ja, dann ist dies äquivalent zu 1) Zentrieren von X; 2) Berechnen der SVD des zentrierten X, sagen wir UDV '; 3) Berechnung paarweiser Abstände zwischen den Reihen von U.
vqv
Vielen Dank, dass Sie diese Frage gestellt haben. Ich denke, dass Ihre Formel nicht korrekt ist. Siehe meine Antwort unten.
User603
@vqv Ja, Beispiel-Kovarianzmatrix. Der ursprüngliche Beitrag wurde entsprechend bearbeitet.
Ahfoss
Siehe auch sehr ähnliche Frage stats.stackexchange.com/q/33518/3277 .
ttnphns

Antworten:

21

Ausgehend von Ahfoss '"Succint" -Lösung habe ich die Cholesky-Zerlegung anstelle der SVD verwendet.

cholMaha <- function(X) {
 dec <- chol( cov(X) )
 tmp <- forwardsolve(t(dec), t(X) )
 dist(t(tmp))
}

Es sollte schneller sein, da das Vorwärtslösen eines Dreiecksystems schneller ist als die dichte Matrixmultiplikation mit der inversen Kovarianz ( siehe hier ). Hier sind die Benchmarks für die Lösungen von ahfoss und whuber in verschiedenen Umgebungen:

 require(microbenchmark)
 set.seed(26565)
 N <- 100
 d <- 10

 X <- matrix(rnorm(N*d), N, d)

 A <- cholMaha( X = X ) 
 A1 <- fastPwMahal(x1 = X, invCovMat = solve(cov(X))) 
 sum(abs(A - A1)) 
 # [1] 5.973666e-12  Ressuring!

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X))
Unit: microseconds
expr          min       lq   median       uq      max neval
cholMaha    502.368 508.3750 512.3210 516.8960  542.806   100
fastPwMahal 634.439 640.7235 645.8575 651.3745 1469.112   100
mahal       839.772 850.4580 857.4405 871.0260 1856.032   100

 N <- 10
 d <- 5
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: microseconds
expr          min       lq    median       uq      max neval
cholMaha    112.235 116.9845 119.114 122.3970  169.924   100
fastPwMahal 195.415 201.5620 205.124 208.3365 1273.486   100
mahal       163.149 169.3650 172.927 175.9650  311.422   100

 N <- 500
 d <- 15
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: milliseconds
expr          min       lq     median       uq      max neval
cholMaha    14.58551 14.62484 14.74804 14.92414 41.70873   100
fastPwMahal 14.79692 14.91129 14.96545 15.19139 15.84825   100
mahal       12.65825 14.11171 39.43599 40.26598 41.77186   100

 N <- 500
 d <- 5
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: milliseconds
expr           min        lq      median        uq       max neval
cholMaha     5.007198  5.030110  5.115941  5.257862  6.031427   100
fastPwMahal  5.082696  5.143914  5.245919  5.457050  6.232565   100
mahal        10.312487 12.215657 37.094138 37.986501 40.153222   100

Cholesky scheint also gleichmäßig schneller zu sein.

Matteo Fasiolo
quelle
3
+1 Gut gemacht! Ich schätze die Erklärung, warum diese Lösung schneller ist.
Whuber
Wie gibt maha () Ihnen die paarweise Distanzmatrix statt nur die Distanz zu einem Punkt?
sheß
1
Sie haben recht, nicht wahr? Meine Bearbeitung ist also nicht ganz relevant. Ich werde es löschen, aber vielleicht werde ich eines Tages eine paarweise Version von maha () zum Paket hinzufügen. Vielen Dank für den Hinweis.
Matteo Fasiolo
1
Das wäre schön! Sich auf etwas freuen.
sheß
9

Die Standardformel für den quadratischen Mahalanobis-Abstand zwischen zwei Datenpunkten lautet

D12=(x1-x2)TΣ-1(x1-x2)

wobei ein p × 1- Vektor ist, der der Beobachtung i entspricht . Typischerweise wird die Kovarianzmatrix aus den beobachteten Daten geschätzt. Ohne die Matrixinversion zu zählen, erfordert diese Operation p 2 + p Multiplikationen und p 2 + 2 p Additionen, die jeweils n ( n - 1 ) / 2 Mal wiederholt werden .xichp×1ichp2+pp2+2pn(n-1)/2

Betrachten Sie die folgende Ableitung:

D12=(x1-x2)TΣ-1(x1-x2)=(x1-x2)TΣ-12Σ-12(x1-x2)=(x1TΣ-12-x2TΣ-12)(Σ-12x1-Σ-12x2)=(q1T-q2T)(q1-q2)

wo . Man beachte, dassxTiΣ-1 istqich=Σ-12xich. Dies beruht auf der Tatsache, dassΣ-1xichTΣ-12=(Σ-12xich)T=qichT ist symmetrisch, was aufgrund der Tatsache gilt, dass für jede symmetrische diagonalisierbare MatrixA=PEPT,Σ-12EIN=PEPT

EIN12T=(PE12PT)T=PTTE12TPT=PE12PT=EIN12

Wenn wir und beachten, dass Σ - 1 symmetrisch ist, sehen wir, dass Σ - 1EIN=Σ-1Σ-1 muss auch symmetrisch sein. WennXist dien×pMatrix von Beobachtungen undQist dien×pMatrixso dass dieithReihe vonQistQi, dannQkann kurzbündig ausgedrückt werdenXΣ-1Σ-12Xn×pQ.n×pichthQ.qichQ. . Dies und die vorherigen Ergebnisse implizieren diesXΣ-12

Die einzigen Operationen, die n ( n - 1 ) / 2 malberechnet werden,sind p Multiplikationen und 2 p Additionen (im Gegensatz zu den p 2 + p Multiplikationen und p 2 + 2 p

Dk=ich=1p(Q.kich-Q.ich)2.
n(n-1)/2p2pp2+pp2+2pZusätze in dem obigen Verfahren), in einem Algorithmus führen , die die Rechenkomplexität Ordnung anstelle das ursprüngliche O ( p 2 n 2 ) .O(pn2+p2n)O(p2n2)
require(ICSNP) # for pair.diff(), C implementation

fastPwMahal = function(data) {

    # Calculate inverse square root matrix
    invCov = solve(cov(data))
    svds = svd(invCov)
    invCovSqr = svds$u %*% diag(sqrt(svds$d)) %*% t(svds$u)

    Q = data %*% invCovSqr

    # Calculate distances
    # pair.diff() calculates the n(n-1)/2 element-by-element
    # pairwise differences between each row of the input matrix
    sqrDiffs = pair.diff(Q)^2
    distVec = rowSums(sqrDiffs)

    # Create dist object without creating a n x n matrix
    attr(distVec, "Size") = nrow(data)
    attr(distVec, "Diag") = F
    attr(distVec, "Upper") = F
    class(distVec) = "dist"
    return(distVec)
}
ahfoss
quelle
Interessant. Tut mir leid, ich weiß nicht R. Können Sie erläutern, was pair.diff()funktioniert, und auch ein numerisches Beispiel mit Ausdrucken aller Schritte Ihrer Funktion angeben ? Vielen Dank.
TTNPHNS
Ich habe die Antwort bearbeitet, um die Ableitung einzuschließen, die diese Berechnungen rechtfertigt, aber ich habe auch eine zweite Antwort veröffentlicht, die Code enthält, der viel prägnanter ist.
Ahfoss
7

Versuchen wir das Offensichtliche. Von

Dichj=(xich-xj)Σ-1(xich-xj)=xichΣ-1xich+xjΣ-1xj-2xichΣ-1xj

es folgt, dass wir den Vektor berechnen können

uich=xichΣ-1xich

in Zeit und der MatrixO(p2)

V=XΣ-1X

in , höchstwahrscheinlich unter Verwendung integrierter schneller (parallelisierbarer) Array-Operationen, und bilden Sie dann die Lösung alsO(pn2+p2n)

D=uu-2V

wobei das äußere Produkt in Bezug auf + ist : ( a b ) i j = a i + b j .+(einb)ichj=einich+bj.

Eine RImplementierung entspricht genau der mathematischen Formulierung (und setzt damit voraus, dass tatsächlich invertierbar ist , wenn hier h invers geschrieben wird ):Σ=Var(X)h

mahal <- function(x, h=solve(var(x))) {
  u <- apply(x, 1, function(y) y %*% h %*% y)
  d <- outer(u, u, `+`) - 2 * x %*% h %*% t(x)
  d[lower.tri(d)]
}

Beachten Sie aus Gründen der Kompatibilität mit den anderen Lösungen, dass nur die eindeutigen nicht diagonalen Elemente zurückgegeben werden und nicht die gesamte quadratische Distanzmatrix (symmetrisch, null auf der Diagonale). Scatterplots zeigen, dass die Ergebnisse mit denen von übereinstimmen fastPwMahal.

In C oder C ++ kann RAM wiederverwendet werden und berechnet im laufenden Betrieb , erübrigt jegliche Notwendigkeit zur Zwischenlagerung von u u .uuuu

Timing-Studien mit Bereich von 33 bis 5000 und p im Bereich von 10 bis 100 zeigen, dass diese Implementierung 1,5- bis 5- mal schneller ist als in diesem Bereich. Die Verbesserung wird besser, wenn p und n zunehmen. Folglich können wir erwarten , dass wir für kleinere p überlegen sind . Die Gewinnschwelle liegt bei p = 7 für n 100n335000p101001.55fastPwMahalpnfastPwMahalpp=7n100. Ob sich die gleichen Rechenvorteile dieser einfachen Lösung auch auf andere Implementierungen beziehen, hängt davon ab, wie gut sie vektorisierte Array-Operationen nutzen.

whuber
quelle
Sieht gut aus. Ich nehme an, es könnte noch schneller gemacht werden, wenn man nur die unteren Diagonalen berechnet, obwohl ich mir keine Möglichkeit ausdenken kann, dies in R zu tun, ohne die schnelle Leistung von applyund outer... zu verlieren, außer zu brechen Rcpp.
Ahfoss
Apply / Outer haben keinen Geschwindigkeitsvorteil gegenüber Plain-Vanilla-Loops.
User603
@ user603 Das verstehe ich grundsätzlich - aber mach das Timing. Darüber hinaus besteht der Hauptzweck der Verwendung dieser Konstrukte darin, eine semantische Hilfe für die Parallelisierung des Algorithmus bereitzustellen: Es ist wichtig, wie unterschiedlich sie ausgedrückt werden . (Es kann sich lohnen, an die ursprüngliche Frage zu erinnern, in der nach C / Fortran / usw gesucht wird.) Ahfoss, ich habe darüber nachgedacht, die Berechnung auch auf das untere Dreieck zu beschränken, und bin Rmir einig, dass darin anscheinend nichts zu gewinnen ist.
whuber
5

Wenn Sie die Mahalanobis- Beispielentfernung berechnen möchten , gibt es einige algebraische Tricks, die Sie ausnutzen können. Sie alle führen dazu, paarweise euklidische Entfernungen zu berechnen. Nehmen wir also an, wir können sie dafür verwenden dist(). Lassen bezeichnen die n × p Datenmatrix, die wir annehmen , so zentriert werden , dass ihre Spalten Mittelwert 0 und Rang haben p , so dass die Probe Kovarianzmatrix nichtsingulär ist. (Zentrieren erfordert O ( n p ) Operationen.) Dann wird die Probe Kovarianzmatrix S = X T X / n .Xn×ppO(np)

S=XTX/n.

Die paarweise Mahalanobis Probenabstände von ist das gleiche wie die paarweise euklidischen Distanzen von X L für jede Matrix L erfüllen L L T = S - 1 , zum Beispiel der Quadratwurzel oder Cholesky - Faktor. Dies folgt aus einer linearen Algebra und führt zu einem Algorithmus, der die Berechnung von S , S - 1 und eine Cholesky - Zerlegung erfordert . Die Komplexität im ungünstigsten Fall ist O ( n p 2 + p 3 ) .X

XL
LLLT=S-1SS-1O(np2+p3)

Tiefer beziehen sich diese Abstände zwischen den Probenhauptkomponenten auf Abstände . Let X = U D V T den SVD bezeichnet X . Dann S = V D 2 V TXX=UDVTX und S - 1 / 2 = V D - 1 V T n 1 / 2 . So X S - 1 / 2 = U V T n 1

S=VD2VT/n
S-1/2=VD-1VTn1/2.
und die Stichproben-Mahalanobis-Abstände sind nur die paarweisen euklidischen Abstände vonU,skaliert mit einem Faktor von
XS-1/2=UVTn1/2
U , weil der euklidische Abstand rotationsinvariant ist. Dies führt zu einem Algorithmus, der die Berechnung der SVD vonXerfordert, diedie ungünstigste KomplexitätO(np2) aufweist,wennn>p.nXO(np2)n>p

Hier ist eine R-Implementierung der zweiten Methode, die ich auf dem iPad, mit dem ich diese Antwort schreibe, nicht testen kann.

u = svd(scale(x, center = TRUE, scale = FALSE), nv = 0)$u
dist(u)
# these distances need to be scaled by a factor of n
vqv
quelle
2

Dies ist eine sehr viel prägnantere Lösung. Es basiert immer noch auf der Ableitung mit der Kovarianzmatrix der Quadratwurzel (siehe meine andere Antwort auf diese Frage), verwendet jedoch nur die Basis R und das Statistikpaket. Es scheint etwas schneller zu sein (ungefähr 10% schneller in einigen Benchmarks, die ich ausgeführt habe). Beachten Sie, dass es die Mahalanobis-Distanz im Gegensatz zur quadratischen Maha-Distanz zurückgibt.

fastPwMahal = function(x1,invCovMat) {
  SQRT = with(svd(invCovMat), u %*% diag(d^0.5) %*% t(v))
  dist(x1 %*% SQRT)
}

Diese Funktion erfordert eine inverse Kovarianzmatrix und gibt kein Entfernungsobjekt zurück. Ich vermute jedoch, dass diese abgespeckte Version der Funktion für das Stapeln von Exchange-Benutzern allgemeiner nützlich ist.

ahfoss
quelle
3
Dies könnte durch Ersetzen SQRTdurch die Cholesky-Zerlegung verbessert werden chol(invCovMat).
vqv
1

n2

Wenn Sie nur die Fortran77-Funktionen in der Benutzeroberfläche verwenden, ist Ihr Unterprogramm für andere Benutzer noch portabel genug.

Horst Grünbusch
quelle
1

Es gibt eine sehr einfache Möglichkeit, dies mit dem R-Paket "biotools" zu tun. In diesem Fall erhalten Sie eine quadratische Mahalanobis-Matrix.

#Manly (2004, p.65-66)

x1 <- c(131.37, 132.37, 134.47, 135.50, 136.17)
x2 <- c(133.60, 132.70, 133.80, 132.30, 130.33)
x3 <- c(99.17, 99.07, 96.03, 94.53, 93.50)
x4 <- c(50.53, 50.23, 50.57, 51.97, 51.37)

#size (n x p) #Means 
x <- cbind(x1, x2, x3, x4) 

#size (p x p) #Variances and Covariances
Cov <- matrix(c(21.112,0.038,0.078,2.01, 0.038,23.486,5.2,2.844, 
        0.078,5.2,24.18,1.134, 2.01,2.844,1.134,10.154), 4, 4)

library(biotools)
Mahalanobis_Distance<-D2.dist(x, Cov)
print(Mahalanobis_Distance)
Jalles10
quelle
Können Sie mir bitte erklären, was eine quadratische Distanzmatrix bedeutet? Jeweils: Ich interessiere mich für den Abstand zwischen zwei Punkten / Vektoren. Was sagt also eine Matrix aus?
Ben
1

Dies ist der Code, den meine alte Antwort von einem anderen Thread hierher verschoben hat .

Ich habe lange Zeit eine quadratische symmetrische Matrix paarweiser Mahalanobis-Abstände in SPSS mithilfe eines Hat-Matrix-Ansatzes berechnet, bei dem ein lineares Gleichungssystem gelöst wird (das schneller ist als das Invertieren der Kovarianzmatrix).

Ich bin kein R-Benutzer, also habe ich gerade versucht, dieses Rezept hier in SPSS zusammen mit "meinem" Rezept auf der Basis von 1000 Fällen mit 400 Variablen zu reproduzieren , und ich habe meinen Weg erheblich schneller gefunden.


H

H(n-1)X(XX)-1XX

Zentrieren Sie also die Spalten der Datenmatrix, berechnen Sie die Hutmatrix, multiplizieren Sie sie mit (n-1) und führen Sie die der doppelten Zentrierung entgegengesetzte Operation aus. Sie erhalten die Matrix der quadratischen Mahalanobis-Abstände.

"Doppelte Zentrierung" ist die geometrisch korrekte Umwandlung von quadratischen Abständen (wie Euklidisch und Mahalanobis) in Skalarprodukte, die aus dem geometrischen Schwerpunkt der Datenwolke definiert werden . Diese Operation basiert implizit auf dem Kosinussatz . Stellen Sie sich vor, Sie haben eine Matrix aus quadratischen euklidischen Abständen zwischen Ihren multivariaten Datenpunkten. Sie ermitteln den Schwerpunkt (multivariates Mittel) der Wolke und ersetzen jeden paarweisen Abstand durch das entsprechende Skalarprodukt (Skalarprodukt), das auf den Abständen basierthh2h1h2cos

In unserer Einstellungen ist die „double-Zentrat“ Matrix spezifisch die Hut - Matrix (multipliziert mit n-1), nicht euklidische Skalarprodukt, und die sich ergebende Matrix quadrierte Abstand wird somit die quadrierten Matrix Mahalanobis - Distanz, nicht euklidische Distanzmatrix quadriert.

HH(n-1)H= {H,H,...}Dmeinheinl2=H+H-2H(n-1)

Der Code in SPSS und Speed ​​Probe ist unten.


Dieser erste Code entspricht der @ahfoss-Funktion fastPwMahalder angegebenen Antwort . Es ist mathematisch äquivalent dazu. Aber ich berechne die gesamte symmetrische Distanzmatrix (über Matrixoperationen), während @ahfoss ein Dreieck der symmetrischen Matrix (Element für Element) berechnet.

matrix. /*Matrix session in SPSS;
        /*note: * operator means matrix multiplication, &* means usual, elementwise multiplication.
get data. /*Dataset 1000 cases x 400 variables
!cov(data%cov). /*compute usual covariances between variables [this is my own matrix function].
comp icov= inv(cov). /*invert it
call svd(icov,u,s,v). /*svd
comp isqrcov= u*sqrt(s)*t(v). /*COV^(-1/2)
comp Q= data*isqrcov. /*Matrix Q (see ahfoss answer)
!seuclid(Q%m). /*Compute 1000x1000 matrix of squared euclidean distances;
               /*computed here from Q "data" they are the squared Mahalanobis distances.
/*print m. /*Done, print
end matrix.

Time elapsed: 3.25 sec

Das Folgende ist meine Modifikation, um es schneller zu machen:

matrix.
get data.
!cov(data%cov).
/*comp icov= inv(cov). /*Don't invert.
call eigen(cov,v,s2). /*Do sdv or eigen decomposition (eigen is faster),
/*comp isqrcov= v * mdiag(1/sqrt(s2)) * t(v). /*compute 1/sqrt of the eigenvalues, and compose the matrix back, so we have COV^(-1/2).
comp isqrcov= v &* (make(nrow(cov),1,1) * t(1/sqrt(s2))) * t(v). /*Or this way not doing matrix multiplication on a diagonal matrix: a bit faster .
comp Q= data*isqrcov.
!seuclid(Q%m).
/*print m.
end matrix.

Time elapsed: 2.40 sec

X(XX)-1X(XX)-1Xsolve(X'X,X')

matrix.
get data.
!center(data%data). /*Center variables (columns).
comp hat= data*solve(sscp(data),t(data))*(nrow(data)-1). /*hat matrix, and multiply it by n-1 (i.e. by df of covariances).
comp ss= diag(hat)*make(1,ncol(hat),1). /*Now using its diagonal, the leverages (as column propagated into matrix).
comp m= ss+t(ss)-2*hat. /*compute matrix of squared Mahalanobis distances via "cosine rule".
/*print m.
end matrix.

[Notice that if in "comp ss" and "comp m" lines you use "sscp(t(data))",
 that is, DATA*t(DATA), in place of "hat", you get usual sq. 
 euclidean distances]

Time elapsed: 0.95 sec
ttnphns
quelle
0

Die von Ihnen veröffentlichte Formel berechnet nicht das, was Sie zu berechnen glauben (eine U-Statistik).

In dem Code, den ich gepostet habe, verwende ich cov(x1)als Skalierungsmatrix (dies ist die Varianz der paarweisen Unterschiede der Daten). Sie verwenden cov(x0)(dies ist die Kovarianzmatrix Ihrer Originaldaten). Ich denke, das ist ein Fehler von Ihrer Seite. Der springende Punkt bei der Verwendung der paarweisen Unterschiede ist, dass Sie nicht davon ausgehen müssen, dass die multivariate Verteilung Ihrer Daten symmetrisch um ein Symmetriezentrum ist (oder dass Sie dieses Symmetriezentrum für diese Angelegenheit schätzen müssen, da crossprod(x1)es proportional zu istcov(x1) ). Offensichtlich cov(x0)verlieren Sie das, indem Sie verwenden.

Dies wird ausführlich in dem Artikel erklärt, auf den ich in meiner ursprünglichen Antwort verwiesen habe.

user603
quelle
1
Ich denke, wir sprechen hier über zwei verschiedene Dinge. Meine Methode berechnet die Mahalanobis-Distanz, die ich anhand einiger anderer Formeln überprüft habe. Meine Formel wurde nun auch unabhängig von Matteo Fasiolound (nehme ich an) whuberin diesem Thread verifiziert . Dein ist anders. Ich wäre daran interessiert zu verstehen, was Sie berechnen, aber es unterscheidet sich deutlich von der Mahalanobis-Distanz, wie sie normalerweise definiert wird.
Ahfoss
@ahfoss: 1) Mahalanobis ist der Abstand des X zu einem Symmetriepunkt in ihrer Metrik. In Ihrem Fall sind die X eine * (n-1) / 2-Matrix oder paarweise Differenzen, ihr Symmetriezentrum ist der Vektor 0_p und ihre Metrik ist das, was ich in meinem Code cov (X1) nannte. 2) Fragen Sie sich, warum Sie überhaupt eine U-Statistik verwenden, und wie im Artikel erklärt, wird die Verwendung von cov (x0) diesen Zweck zunichte machen.
user603
XXOp
cov(x0)SGSτLQ.D