Odległości Mahalanobisa parami

18

Muszę obliczyć przykładową odległość Mahalanobisa w R pomiędzy każdą parą obserwacji w macierzy współzmiennych . Potrzebuję rozwiązania, które jest wydajne, tj. Obliczane są tylko odległości, a najlepiej realizowane w C / RCpp / Fortran itp. Zakładam, że , macierz kowariancji populacyjnej, jest nieznana i wykorzystuję próbkę macierz kowariancji na swoim miejscu. $n \times p$ $n(n-1)/2$ $\Sigma$

Szczególnie interesuje mnie to pytanie, ponieważ wydaje się, że nie ma metody „konsensusu” do obliczania par Mahalanobisa w parach odległości R, tj. Nie jest ona zaimplementowana distani w funkcji, ani w cluster::daisyfunkcji. Ta mahalanobisfunkcja nie oblicza odległości parami bez dodatkowej pracy programisty.

Zostało to już zadane tutaj Odległość Pairwise Mahalanobis w R , ale rozwiązania tam wydają się nieprawidłowe.

Oto poprawna, ale strasznie nieefektywna (ponieważ obliczane są odległości $n \times 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))))

Łatwo jest to napisać w C, ale wydaje mi się, że coś podstawowego powinno mieć wcześniej istniejące rozwiązanie. Czy jest jeden

Istnieją inne rozwiązania, które nie są wystarczające : HDMD::pairwise.mahalanobis()oblicza $n \times n$ odległości, gdy wymagane są tylko $n(n-1)/2$ unikalne odległości. compositions::MahalanobisDist()wydaje się obiecujące, ale nie chcę, aby moja funkcja pochodziła z pakietu, który zależy rgl, co poważnie ogranicza możliwość uruchamiania mojego kodu przez innych. Jeśli ta implementacja nie jest idealna, wolę napisać własną. Czy ktoś ma doświadczenie w tej funkcji?

r algorithms distance

— ahfoss
źródło

Witamy. Czy potrafisz wydrukować dwie macierze odległości w swoim pytaniu? A co jest dla ciebie „nieefektywne”?

— ttnphns

1

Czy używasz tylko przykładowej macierzy kowariancji? Jeśli tak, to odpowiada to 1) centrowaniu X; 2) obliczanie SVD wyśrodkowanego X, powiedzmy UDV '; 3) obliczanie parami odległości między rzędami U.

— vqv

Dziękujemy za opublikowanie tego jako pytania. Myślę, że twoja formuła jest nieprawidłowa. Zobacz moją odpowiedź poniżej.

— user603

@vqv Tak, przykładowa macierz kowariancji. Oryginalny post jest edytowany, aby to odzwierciedlić.

— ahfoss,

Zobacz także bardzo podobne pytanie stats.stackexchange.com/q/33518/3277 .

— ttnphns

21

Zaczynając od rozwiązania „succint” firmy ahfoss, użyłem rozkładu Cholesky'ego zamiast SVD.

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

Powinno być szybciej, ponieważ rozwiązywanie do przodu układu trójkątnego jest szybsze niż gęste mnożenie macierzy z odwrotną kowariancją ( patrz tutaj ). Oto testy porównawcze rozwiązań ahfoss i whuber w kilku ustawieniach:

 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 wydaje się być jednakowo szybszy.

— Matteo Fasiolo
źródło

3

+1 Dobra robota! Doceniam wyjaśnienie, dlaczego to rozwiązanie jest szybsze.

— whuber

Jak maha () daje ci macierz odległości w parze, a nie tylko odległość do punktu?

— sheß

1

Masz rację, tak nie jest, więc moja edycja nie jest całkowicie aktualna. Usunę go, ale może kiedyś dodam do pakietu wersję maha (). Dzięki za zwrócenie na to uwagi.

— Matteo Fasiolo

1

Byłoby uroczo! Nie mogę się doczekać.

— sheß

9

Standardowa formuła dla kwadratowej odległości Mahalanobisa między dwoma punktami danych to

D_{12} = (x_{1} - x_{2})^{T} Σ^{- 1} (x_{1} - x_{2})

$D_{12} = (x_1-x_2)^T \Sigma^{-1} (x_1-x_2)$

gdzie jest wektorem odpowiadającym obserwacji . Zazwyczaj macierz kowariancji jest szacowana na podstawie zaobserwowanych danych. Nie licząc inwersję macierzy, ta operacja wymaga mnożenia i dodatki, każdy powtarzane razy. $x_i$ $p \times 1$ $i$ $p^2+p$ $p^2+2p$ $n(n-1)/2$

Rozważ następujące wyprowadzenie:

\begin{array}{rcl} D_{12} & = & (x_{1} - x_{2})^{T} Σ^{- 1} (x_{1} - x_{2}) \\ = & (x_{1} - x_{2})^{T} Σ^{- \frac{1}{2}} Σ^{- \frac{1}{2}} (x_{1} - x_{2}) \\ = & (x_{1}^{T} Σ^{- \frac{1}{2}} - x_{2}^{T} Σ^{- \frac{1}{2}}) (Σ^{- \frac{1}{2}} x_{1} - Σ^{- \frac{1}{2}} x_{2}) \\ = & (q_{1}^{T} - q_{2}^{T}) (q_{1} - q_{2}) \end{array}

$\begin{eqnarray*} D_{12} &=& (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) \\ &=& (x_1-x_2)^T \Sigma^{-\frac{1}{2}} \Sigma^{-\frac{1}{2}} (x_1-x_2) \\ &=& (x_1^T \Sigma^{-\frac{1}{2}} - x_2^T \Sigma^{-\frac{1}{2}}) (\Sigma^{-\frac{1}{2}}x_1 - \Sigma^{-\frac{1}{2}}x_2) \\ &=& (q_1^T - q_2^T)(q_1 - q_2) \end{eqnarray*}$

gdzie . Zauważ, że $q_i = \Sigma^{-\frac{1}{2}}x_i$ . Zależy to od faktu, że $x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$ jest symetryczny, co wynika z faktu, że dla dowolnej symetrycznej macierzy diagonalnej, $\Sigma^{-\frac{1}{2}}$ $A = PEP^T$

\begin{array}{rcl} A^{{\frac{1}{2}}^{T}} & = & (P E^{\frac{1}{2}} P^{T})^{T} \\ = & P^{T^{T}} E^{{\frac{1}{2}}^{T}} P^{T} \\ = & P E^{\frac{1}{2}} P^{T} \\ = & A^{\frac{1}{2)}} \end{array}

$\begin{eqnarray*} A^{\frac{1}{2}^T} &=& (PE^{\frac{1}{2}}P^T)^T \\ &=& P^{T^T} E^{\frac{1}{2}^T} P^T \\ &=& PE^{\frac{1}{2}}P^T \\ &=& A^{\frac{1}{2}} \end{eqnarray*}$

Jeśli pozwolimy i zauważymy , że jest symetryczny, widzimy, że $A=\Sigma^{-1}$ $\Sigma^{-1}$ musi być również symetryczne. Jeślijestmacierząobserwacji, ajestmacierzątaką, żerzęduwynosi, wówczasmożna zwięźle wyrazić jako $\Sigma^{-\frac{1}{2}}$ $X$ $n \times p$ $Q$ $n \times p$ $i^{th}$ $Q$ $q_i$ $Q$ . To i poprzednie wyniki implikują to $X\Sigma^{-\frac{1}{2}}$

jedynymi operacjami, które są obliczane razy, sąmnożenia i dodawania (w przeciwieństwie domnożenia oraz

D_{k ℓ} = \sum_{i = 1}^{p} (Q_{k i} - Q_{ℓ i})^{2} .

$D_{k\ell} = \sum_{i=1}^p (Q_{ki}-Q_{\ell i})^2.$

n (n - 1) / 2

$n(n-1)/2$

p

$p$

2 p

$2p$

p^{2} + p

$p^2+p$

p^{2} + 2 p

$p^2+2p$ uzupełnienia w powyższej metodzie), w wyniku czego powstaje algorytm o złożoności obliczeniowej

zamiast pierwotnego

.

O (p n^{2} + p^{2} n)

$O(pn^2 + p^2n)$

O (p^{2} n^{2})

$O(p^2n^2)$

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
źródło

Ciekawy. Przepraszam, nie wiem R. Czy możesz wydać co to pair.diff()znaczy, a także podać numeryczny przykład z wydrukami każdego kroku twojej funkcji? Dzięki.

— ttnphns

Zredagowałem odpowiedź, aby uwzględnić wyprowadzenie uzasadniające te obliczenia, ale zamieściłem również drugą odpowiedź zawierającą kod, który jest znacznie bardziej zwięzły.

— ahfoss,

7

Spróbujmy tego, co oczywiste. Od

D_{i j} = (x_{i} - x_{j})^{'} Σ^{- 1} (x_{i} - x_{j}) = x_{i}^{'} Σ^{- 1} x_{i} + x_{j}^{'} Σ^{- 1} x_{j} - 2 x_{i}^{'} Σ^{- 1} x_{j}

$D_{ij} = (x_i-x_j)^\prime \Sigma^{-1} (x_i-x_j)=x_i^\prime \Sigma^{-1}x_i + x_j^\prime \Sigma^{-1}x_j -2 x_i^\prime \Sigma^{-1}x_j$

wynika z tego, że możemy obliczyć wektor

u_{i} = x_{i}^{'} Σ^{- 1} x_{i}

$u_i = x_i^\prime \Sigma^{-1}x_i$

w czasie i macierzy $O(p^2)$

V = X Σ^{- 1} X^{'}

$V = X \Sigma^{-1} X^\prime$

w czasie , najprawdopodobniej przy użyciu wbudowanych szybkich (równoległych) operacji tablicowych, a następnie utwórz rozwiązanie jako $O(p n^2 + p^2 n)$

D = u \oplus u - 2 V

$D = u \oplus u - 2 V$

gdzie jest iloczynem zewnętrznym w odniesieniu do : $\oplus$ $+$ $(a \oplus b)_{ij} = a_i + b_j.$

RRealizacja zwięźle paralele sformułowanie matematycznego (a zakłada się z nim, że rzeczywiście jest odwracalna z odwrotnym pisemnej tutaj): $\Sigma=\text{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)]
}

Należy zwrócić uwagę, że w celu zapewnienia zgodności z innymi rozwiązaniami zwracane są tylko unikalne elementy o przekątnej, a nie cała kwadratowa macierz odległości (symetryczna, zero na przekątnej). Wykresy rozrzutu pokazują, że jego wyniki są zgodne z wynikami fastPwMahal.

W języku C i C ++, pamięć RAM może być ponownie użyty, a oblicza się na bieżąco, eliminując potrzebę stosowania pośredniego składowania . $u\oplus u$ $u\oplus u$

Badania czasowe z zakresie od do i zakresie od do wskazują, że ta implementacja jest do razy szybsza niż w tym zakresie. Poprawa poprawia się wraz ze wzrostem wartości i . W związku z tym możemy spodziewać się wyższego poziomu dla mniejszych . Próg rentowności występuje w okolicach dla $n$ $33$ $5000$ $p$ $10$ $100$ $1.5$ $5$ fastPwMahal $p$ $n$ fastPwMahal $p$ $p=7$ $n\ge 100$ . To, czy te same zalety obliczeniowe tego prostego rozwiązania dotyczą innych implementacji, może zależeć od tego, jak dobrze wykorzystują one wektoryzowane operacje tablicowe.

— Whuber
źródło

Wygląda dobrze. Zakładam, że można by to uczynić jeszcze szybszym, obliczając tylko dolne przekątne, chociaż nie mogę wymyślić sposobu na zrobienie tego w R bez utraty szybkiej wydajności applyi outer... z wyjątkiem wybuchu Rcpp.

— ahfoss

Apply / outer nie mają przewagi prędkości nad zwykłymi waniliowymi pętlami.

— user603,

@ user603 Rozumiem to w zasadzie - ale wykonuj wyczucie czasu. Co więcej, głównym celem korzystania z tych konstrukcji jest zapewnienie pomocy semantycznej dla równoległości algorytmu: ważna jest różnica w sposobie ich wyrażania . (Być może warto przypomnieć oryginalne pytanie, które dotyczy implementacji C / Fortran / itp.) Ahfoss, pomyślałem o ograniczeniu obliczeń również do dolnego trójkąta i zgadzam się, że Rnie wydaje się, że nic z tego można zyskać.

— whuber

5

Jeśli chcesz obliczyć przykładową odległość Mahalanobisa, istnieje kilka sztuczek algebraicznych, które możesz wykorzystać. Wszystkie prowadzą do obliczenia par euklidesowych odległości, więc załóżmy, że możemy dist()do tego użyć . Niech oznacza macierz danych , która, jak zakładamy, jest wyśrodkowana, tak że jej kolumny mają średnią 0, i ma rangę tak że macierz kowariancji próbki nie jest pojedyncza. (Centrowanie wymaga operacji .) Następnie macierz kowariancji próbki to $X$ $n\times p$ $p$ $O(np)$

S = X^{T} X / n .

$S = X^T X / n.$

Próbki Mahalanobisa w parach są takie same, jak pary euklidesowej odległości $X$ dla każdej macierzy zgodnej , na przykład pierwiastka lub czynnik Choleskiego. Wynika to z pewnej algebry liniowej i prowadzi do algorytmu wymagającego obliczenia , i rozkładu Choleskiego. W najgorszym przypadku złożoność to .

X L

$X L$

L

$L$

L L^{T} = S^{- 1}

$LL^T = S^{-1}$

S

$S$

S^{- 1}

$S^{-1}$

O (n p^{2} + p^{3})

$O(np^2 + p^3)$

$X$ $X=UDV^T$ $X$

S = V D^{2} V^{T} / n

$S=VD^2V^T/n$

S^{- 1 / 2} = V D^{- 1} V^{T} n^{1 / 2} .

$S^{-1/2}=VD^{-1}V^T n^{1/2}.$

X S^{- 1 / 2} = U V^{T} n^{1 / 2}

$X S^{-1/2} = UV^T n^{1/2}$

U

$U$

\sqrt{n}

$\sqrt{n}$

X

$X$

O (n p^{2})

$O(n p^2)$

n > p

$n>p$

Oto implementacja R drugiej metody, której nie mogę przetestować na iPadzie, której używam do napisania tej odpowiedzi.

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
źródło

2

To jest znacznie bardziej zwięzłe rozwiązanie. Wciąż opiera się na wyprowadzeniu obejmującym macierz kowariancji odwrotnego pierwiastka kwadratowego (patrz moja inna odpowiedź na to pytanie), ale używa tylko podstawy R i pakietu statystyk. Wydaje się być nieco szybszy (około 10% szybciej w niektórych testach, które przeprowadziłem). Zauważ, że zwraca odległość Mahalanobisa, w przeciwieństwie do kwadratowej odległości Maha.

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

Ta funkcja wymaga odwrotnej macierzy kowariancji i nie zwraca obiektu odległości - ale podejrzewam, że ta zredukowana wersja funkcji będzie bardziej użyteczna do układania użytkowników na stosie.

— ahfoss
źródło

3

Można to poprawić, zastępując SQRTrozkład Cholesky'ego chol(invCovMat).

— vqv

1

$n^2$

Jeśli używasz tylko funkcji Fortran77 w interfejsie, twój podprogram jest wciąż wystarczająco przenośny dla innych.

— Horst Grünbusch
źródło

1

Jest to bardzo prosty sposób, aby to zrobić za pomocą pakietu R „biotools”. W takim przypadku otrzymasz Matrycę Mahalanobisa o kwadracie odległości.

#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
źródło

Czy możesz mi wyjaśnić, co oznacza kwadratowa macierz odległości? Odpowiednio: interesuje mnie odległość między dwoma punktami / wektorami, więc co mówi macierz?

— Ben,

1

To jest ten rozszerzony kod, który moja stara odpowiedź przeniosła tutaj z innego wątku .

Od dłuższego czasu wykonuję obliczenia kwadratowej macierzy symetrycznej par Mahalanobisa w parach odległości w SPSS metodą macierzy kapelusza, stosując rozwiązanie układu równań liniowych (ponieważ jest ono szybsze niż odwracanie macierzy kowariancji).

Nie jestem użytkownikiem R, więc właśnie próbowałem odtworzyć @ahfoss ' ten przepis tutaj w SPSS wraz z „moim” przepisem, na danych 1000 przypadków przez 400 zmiennych, i znalazłem swoją drogę znacznie szybciej.

Szybszy sposób na obliczenie pełnej macierzy par Mahalanobisa w parach jest już zakończony $\bf H$

$\mathbf{H}(n-1)$ $\bf X(X'X)^{-1}X'$ $\bf X$

Tak więc wyśrodkuj kolumny macierzy danych, oblicz macierz kapelusza, pomnóż ją przez (n-1) i wykonaj operację przeciwną do podwójnego centrowania. Otrzymujesz macierz kwadratowych odległości Mahalanobisa.

$h$ $h^2$ $h_1h_2\cos$

W naszych ustawieniach „podwójnie centrowana” macierz jest w szczególności macierzą kapelusza (pomnożoną przez n-1), a nie euklidesowymi produktami skalarnymi, a wynikowa kwadratowa macierz odległości jest zatem kwadratową macierzą odległości Mahalanobisa, a nie kwadratową macierzą odległości euklidesową.

$H$ $\mathbf{H}(n-1)$ H= {H,H,...} $\mathbf {D_{mahal}^2} = H+H'-2 \mathbf{H}(n-1)$

Kod w SPSS i czujniku prędkości znajduje się poniżej.

Ten pierwszy kod odpowiada @ahfoss funkcji fastPwMahalw cytowanej odpowiedzi . Jest to równoważne matematycznie. Ale obliczam pełną symetryczną macierz odległości (za pomocą operacji macierzowych), podczas gdy @ahfoss obliczył trójkąt macierzy symetrycznej (element po elemencie).

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

Oto moja modyfikacja, aby przyspieszyć:

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

$\bf X(X'X)^{-1}X'$ $\bf (X'X)^{-1}X'$ solve(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
źródło

0

Opublikowana formuła nie oblicza tego, co według Ciebie obliczasz (statystyki U).

W opublikowanym przeze mnie kodzie używam cov(x1)jako macierzy skalowania (jest to wariancja różnic par danych). Używasz cov(x0)(jest to macierz kowariancji oryginalnych danych). Myślę, że to błąd z twojej strony. Chodzi o to, żeby wykorzystać różnice par, ponieważ odciąża cię to od założenia, że wielowymiarowy rozkład twoich danych jest symetryczny wokół środka symetrii (lub żeby oszacować to centrum symetrii pod tym względem, ponieważ crossprod(x1)jest proporcjonalne cov(x1)). Oczywiście przez użycie cov(x0)tracisz to.

Jest to dobrze wyjaśnione w artykule, do którego odsyłam w mojej oryginalnej odpowiedzi.

— użytkownik603
źródło

1

Myślę, że mówimy tutaj o dwóch różnych rzeczach. Moja metoda oblicza odległość Mahalanobisa, którą zweryfikowałem na podstawie kilku innych formuł. Moja formuła została teraz również niezależnie zweryfikowana przez Matteo Fasiolo(i zakładam) whuberw tym wątku. Twój jest inny. Byłbym zainteresowany zrozumieniem tego, co obliczasz, ale wyraźnie różni się on od dystansu Mahalanobisa, jak zwykle definiuje się.

— ahfoss,

@ahfoss: 1) mahalanobis to odległość X do punktu symetrii w ich metryce. W twoim przypadku X są macierzą * (n-1) / 2 różnic par, ich środkiem symetrii jest wektor 0_p, a ich metryką jest to, co nazwałem cov (X1) w moim kodzie. 2) zadaj sobie pytanie, dlaczego w ogóle używasz statystyki U, a jak wyjaśnia papier, zobaczysz, że użycie cov (x0) pokonuje ten cel.

— user603

X

$X$

X

$X$

O_{p}

$O_p$

cov(x0)

S

$S$

G S

$GS$

τ

$\tau$

L Q D

$LQD$