Мне нужно рассчитать выборочное расстояние Махаланобиса в R между каждой парой наблюдений в матрице ковариат . Мне нужно решение, которое является эффективным, то есть только Е. Рассчитываются расстояний, и желательно, чтобы они были реализованы в C / RCpp / Fortran и т. Д. Я предполагаю, что , ковариационная матрица населенности, неизвестна, и использую выборочную ковариацию матрица на своем месте.
Я особенно заинтересован в этом вопросе, так как, похоже, не существует "консенсусного" метода для расчета попарных расстояний Махаланобиса в R, т.е. он не реализован dist
ни в функции, ни в cluster::daisy
функции. mahalanobis
Функция не вычисляет попарные расстояния без дополнительной работы от программиста.
Об этом уже спрашивали попарно расстояние Махаланобиса в R , но решения там кажутся неверными.
Вот правильный, но ужасно неэффективный (так как как рассчитывается 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))))
Это достаточно просто для написания кода на C, но я чувствую, что у этого базового решения должно быть уже существующее решение. Есть один?
Существуют и другие решения, которые не дотягивают: HDMD::pairwise.mahalanobis()
вычисляются расстояний, когда требуются только уникальных расстояния. compositions::MahalanobisDist()
кажется многообещающим, но я не хочу, чтобы моя функция исходила из пакета, от rgl
которого зависит , что серьезно ограничивает способность других выполнять мой код. Если эта реализация не идеальна, я бы лучше написал свою. У кого-нибудь есть опыт работы с этой функцией?
источник
Ответы:
Исходя из «сукцинтованного» решения Ахфосса, я использовал декомпозицию Холецкого вместо СВД.
Это должно быть быстрее, потому что решение вперед треугольной системы быстрее, чем умножение плотной матрицы с обратной ковариацией ( см. Здесь ). Вот эталонные тесты решений Ahfoss и Whuber в нескольких ситуациях:
Так что Холецкий, кажется, быстрее всех.
источник
Стандартная формула для квадрата расстояния Махаланобиса между двумя точками данных
где - вектор p × 1, соответствующий наблюдению i . Как правило, ковариационная матрица оценивается по наблюдаемым данным. Не считая обращения матрицы, эта операция выполняется р 2 + р умножений и р 2 + 2 р дополнения, каждый из которых повторных п ( п - 1 ) / 2 раз.xi p×1 i p2+p p2+2p n(n−1)/2
Рассмотрим следующий вывод:
где . Обратите внимание, чтоxTiΣ-1Qя= Σ- 12Икся . Это зависит от того факта, чтоΣ-1ИксTяΣ- 12= ( Σ- 12Икся)T= qTя Σ- 12 симметрична, что имеет место в силу того, что для любой симметричной диагонализуемой матрицы ,A = PЕпT
Если мы допустим и заметим, что Σ - 1 симметрична, мы увидим, чтоA = Σ- 1 Σ- 1 Σ- 12 must also be symmetric. If Икс is the n × p matrix of observations and Q is the n × p matrix such that the ят ч row of Q is Qя , then Q can be succinctly expressed as ИксΣ- 12 . This and the previous results imply that
источник
pair.diff()
делает, а также привести числовой пример с распечатками каждого шага вашей функции? Благодарю.Давайте попробуем очевидное. Из
Отсюда следует, что мы можем вычислить вектор
inO(p2) time and the matrix
inO(pn2+p2n) time, most likely using built-in fast (parallelizable) array operations, and then form the solution as
where⊕ is the outer product with respect to + : (a⊕b)ij=ai+bj.
AnΣ=Var(X) actually is invertible with inverse written h here):
R
implementation succinctly parallels the mathematical formulation (and assumes, with it, thatNote, for compability with the other solutions, that only the unique off-diagonal elements are returned, rather than the entire (symmetric, zero-on-the-diagonal) squared distance matrix. Scatterplots show its results agree with those of
fastPwMahal
.In C or C++, RAM can be re-used andu⊕u computed on the fly, obviating any need for intermediate storage of u⊕u .
Timing studies withn ranging from 33 through 5000 and p ranging from 10 to 100 indicate this implementation is 1.5 to 5 times faster than p and n increase. Consequently, we can expect p . The break-even occurs around p=7 for n≥100 . Whether the same computational advantages of this straightforward solution pertain in other implementations may be a matter of how well they take advantage of vectorized array operations.
fastPwMahal
within that range. The improvement gets better asfastPwMahal
to be superior for smallerисточник
apply
andouter
... except for breaking outRcpp
.R
there seems to be nothing to gain by that.If you wish to compute the sample Mahalanobis distance, then there are some algebraic tricks that you can exploit. They all lead to computing pairwise Euclidean distances, so let's assume we can useX denote the n×p data matrix, which we assume to be centered so that its columns have mean 0 and to have rank p so that the sample covariance matrix is nonsingular. (Centering requires O(np) operations.) Then the sample covariance matrix is
dist()
for that. LetThe pairwise sample Mahalanobis distances ofX is the same as the pairwise Euclidean distances of
More deeply, these distances relate to distances between the sample principal components ofX . Let X=UDVT denote the SVD of X . Then
Here is an R implementation of the second method which I cannot test on the iPad I am using to write this answer.
источник
Это гораздо более краткое решение. Он по-прежнему основан на выводе с использованием обратной ковариационной матрицы с квадратным корнем (см. Мой другой ответ на этот вопрос), но использует только базу R и пакет статистики. Кажется, что это немного быстрее (примерно на 10% быстрее в некоторых тестах, которые я проводил). Обратите внимание, что он возвращает расстояние Махаланобиса, а не квадрат Маха.
Эта функция требует обратной ковариационной матрицы и не возвращает объект расстояния - но я подозреваю, что эта урезанная версия функции будет более полезна для обмена стеками пользователей.
источник
SQRT
разложение Холецкогоchol(invCovMat)
.У меня была похожая проблема, решенная написанием подпрограммы на Fortran95. Как и вы, я не хотел вычислять дубликаты средиN2 расстояния. Скомпилированный Fortran95 почти так же удобен с базовыми матричными вычислениями, как R или Matlab, но намного быстрее с циклами. Подпрограммы для разложения Холецкого и замены треугольника могут быть использованы из LAPACK.
Если вы используете в интерфейсе только функции Fortran77, ваша подпрограмма все еще достаточно переносима для других.
источник
Есть очень простой способ сделать это с помощью R Package "biotools". В этом случае вы получите Квадратную Матрицу Махаланобиса.
источник
Это расширенный код, мой старый ответ перенесен сюда из другого потока .
Я долгое время занимался вычислением квадратно-симметричной матрицы парных расстояний Махаланобиса в SPSS с помощью матричного подхода с использованием решения системы линейных уравнений (поскольку это быстрее, чем инвертирование ковариационной матрицы).
Я не пользователь R, поэтому я просто попытался воспроизвести этот рецепт @ahfoss здесь, в SPSS, вместе с «моим» рецептом на данных 1000 случаев по 400 переменным, и я нашел свой путь значительно быстрее.
Более быстрый способ расчета полной матрицы попарных расстояний Махаланобиса - через матрицуЧАС , Я имею в виду, что если вы используете язык высокого уровня (например, R) с довольно быстрыми встроенными функциями умножения матриц и инверсии, вам не понадобятся циклы вообще, и это будет быстрее, чем выполнение регистровых циклов.
Определение . Двойной центрируется матрица квадратов попарных расстояний Махаланобиса равноH (n-1) где матрица шляпы X ( X'Х )- 1Икс' , рассчитывается по центру столбца данных Икс ,
Таким образом, отцентрируйте столбцы матрицы данных, вычислите головную матрицу, умножьте на (n-1) и выполните операцию, противоположную двойному центрированию. Вы получаете матрицу квадратов расстояний Махаланобиса.
«Двойное центрирование» - это геометрически правильное преобразование квадратов расстояний (таких как Евклидово и Махаланобис) в скалярные произведения, определенные из геометрического центроида облака данных. Эта операция неявно основана на теореме косинуса . Представьте, что у вас есть матрица квадратов евклидовых расстояний между вашими многовариантными точками данных. Вы находите центроид (многовариантное среднее) облака и заменяете каждое попарное расстояние на соответствующее скалярное произведение (точечное произведение), оно основано на расстоянияхчас s к центроиду и углу между этими векторами, как показано в ссылке. час2 s стоят на диагонали этой матрицы скалярных произведений и час1час2соз являются недиагональными записями. Затем, используя непосредственно формулу теоремы косинуса, вы легко конвертируете матрицу «двойного центрирования» обратно в матрицу квадратов расстояний.
В наших настройках матрица "двойного центра" - это, в частности, матрица шляпы (умноженная на n-1), а не евклидовы скалярные произведения, и результирующая квадратная матрица расстояний, таким образом, является квадратом матрицы расстояния Махаланобиса, а не квадратом евклидовой матрицы расстояний.
В матричной записи: пустьЧАС быть диагональю H (n-1) столбец вектор. Распространить столбец в квадратную матрицу D2м а ч а л= H+ H'- 2 H ( n - 1 ) ,
H= {H,H,...}
:; тогдаКод в SPSS и датчик скорости ниже.
Этот первый код соответствует @ahfoss функции
fastPwMahal
из процитировал ответ . Это эквивалентно математически. Но я вычисляю полную симметричную матрицу расстояний (через матричные операции), в то время как @ahfoss вычисляет треугольник симметричной матрицы (элемент за элементом).Вот моя модификация, чтобы сделать это быстрее:
И наконец, «подход с использованием шляпной матрицы». Для скорости я вычисляю матрицу шляпы (данные должны быть сначала центрированы)X ( X'Х )- 1Икс' через обобщенное обратное ( X'Х )- 1Икс' полученный в линейном системном решателе
solve(X'X,X')
.источник
The formula you have posted is not computing what you think you are computing (a U-statistics).
В коде, который я разместил, я использую в
cov(x1)
качестве матрицы масштабирования (это дисперсия парных разностей данных). Вы используетеcov(x0)
(это ковариационная матрица ваших исходных данных). Я думаю, что это ошибка с вашей стороны. Весь смысл использования парных различий состоит в том, что это избавляет вас от предположения, что многомерное распределение ваших данных симметрично относительно центра симметрии (или для того, чтобы оценить этот центр симметрии для этого вопроса, посколькуcrossprod(x1)
он пропорционаленcov(x1)
). Очевидно, что при использованииcov(x0)
вы теряете это.Это хорошо объясняется в статье, на которую я ссылаюсь в своем первоначальном ответе.
источник
Matteo Fasiolo
и (я предполагаю)whuber
в этой теме. Твой другой. Мне было бы интересно понять, что вы рассчитываете, но оно явно отличается от расстояния Махаланобиса, которое обычно определяется.cov(x0)
обычно используется в этом контексте и, кажется, согласуется с Croux et использование ал. В статье не упоминается U-статистика , по крайней мере, явно. Они упоминают