Реализация многомерной гауссовской функции плотности вероятности для> 2 измерений в C ++

Question

Я работаю над реализацией функции плотности вероятности многомерного гауссиана в C ++, и я застрял в том, как лучше всего справиться со случаями, когда размерность> 2.

PDF гауссова может бытьзаписывается как

, где (A) 'или A' представляет транспонирование «матрицы», созданной путем вычитания среднего из всех элементов x.В этом уравнении k - это число измерений, которое у нас есть, а сигма - ковариационная матрица, которая является матрицей akxk.Наконец, | X |означает определитель матрицы X.

В одномерном случае реализация pdf тривиальна.Даже в двумерном (k = 2) случае это тривиально.Однако, когда мы пойдем дальше двух измерений, реализация будет намного сложнее.

В двумерном случае у нас будет

, где rho - этокорреляция между x и y, с корреляцией, равной

В этом случае я мог бы использовать Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> для реализации первого уравнения или просто вычислить все сам, используя второе уравнениебез использования упрощенного интерфейса линейной алгебры Эйгена.

Мои мысли о попытке многомерного случая, вероятно, начнутся с расширения приведенных выше уравнений до многомерного случая

с

Мои вопросы:

1. Было бы целесообразно / рекомендовано использовать boost::multi_array для n-мерного массива,или я должен вместо этого попытаться использовать Eigen?
2. Должны ли я иметь отдельные функции для одномерных / двумерных случаев или просто абстрагировать все это в многовариантный случай с использованием boost :: multi_array (или соответствующей альтернативы)?

Answer 1

Я немного не в своем духе, но некоторые мысли:

Во-первых, с точки зрения программирования, стандартный ответ - "профиль".То есть сначала закодируйте его более понятным способом.Затем профилируйте исполнение, чтобы увидеть, стоит ли оптимизировать.ИМХО, вероятно, яснее использовать матричную библиотеку, чтобы быть ближе к исходной математике.

С математической точки зрения: я немного сомневаюсь в формуле, которую вы предоставляете для многомерного случая.Это не выглядит правильным для меня.Выражение Z должно быть квадратичной формы, а ваше Z - нет.Если только я что-то упустил.

Вот вариант, который вы не упомянули, но может иметь смысл.Особенно, если вы собираетесь оценивать PDF несколько раз для одного дистрибутива.Начните с расчета основных компонентов вашего дистрибутива.То есть собственное основание для Σ.Основные направления компонентов являются ортогональными.В базисе основных компонентов все кросс-ковариации равны 0, поэтому PDF имеет простую форму.Если вы хотите оценить, измените базу на основе входных данных на базис основных компонентов, а затем выполните более простое вычисление PDF для этого.

Мысль заключается в том, что вы можете рассчитать изменение базисной матрицы и основных компонентов один развпереди, а затем нужно выполнить только одно матричное умножение (изменение базиса) для каждой оценки вместо двух умножений матриц, необходимых для оценки (x-μ)' Σ (x-μ) в стандартном базисе.

Answer 2

Я в основном реализовал exp часть уравнения для трехмерного случая в этом вопросе .Сначала я использовал библиотеку компьютерного зрения под названием OpenCV .Но я заметил, что интерфейс C ++ был очень медленным.После этого я попробовал интерфейс C, который был немного быстрее.Наконец, я решил игнорировать гибкость и удобочитаемость, поэтому я реализовал это без каких-либо библиотек, и это было намного быстрее.

Я хочу сказать следующее: когда важна производительность, вы должны рассмотреть возможность реализации особых случаев длянаиболее часто используемое количество измерений с минимальными накладными расходами.В противном случае выберите ремонтопригодность, а не скорость.

Отказ от ответственности: я ничего не знаю о скорости Eigen или boost::multi_array (на что, вероятно, нацелен этот вопрос?).