Implementando una función de densidad de probabilidad gaussiana multivariante para> 2 dimensiones en C++

Question

Estoy trabajando en la implementación de una función de densidad de probabilidad de un gaussiano multivariante en C++, y estoy atascado en cómo manejar mejor los casos donde dimensión> 2.

la pdf de un gaussiano se puede escribir como

donde (a) 'o a' representa la transpuesta de la 'matriz' creado por restando la media de todos los elementos de x. En esta ecuación, k es el número de dimensiones que tenemos, y sigma representa la matriz de covarianza, que es una matriz k x k. Finalmente, | X | significa el determinante de la matriz X.

En el caso univariante, implementar el pdf es trivial. Incluso en el caso bivariante (k = 2), es trivial. Sin embargo, cuando vamos más allá de dos dimensiones, la implementación es mucho más difícil.

En el caso de dos variables, tendríamos

donde rho es la correlación entre X e Y, con correlación igual a

En este caso, I podría usar Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> para implementar la primera ecuación, o simplemente calcular todo usando la segunda ecuación, sin beneficiarse de la interfaz simplificada de álgebra lineal de Eigen.

Mis pensamientos para una tentativa en el caso multivariado, probablemente comenzará mediante la extensión de las ecuaciones anteriores al caso multivariante

con

Mis preguntas son:

1. ¿Sería apropiado/aconsejable usar un boost::multi_array para la matriz n-dimensional , o debería intentar aprovechar Eigen?
2. ¿Debería tener funciones separadas para los casos univariados/bivariados, o debería simplemente resumirlo todo en el caso multivariado usando boost :: multi_array (o una alternativa adecuada)?

Answer 1

estoy un poco fuera de mi elemento aquí, pero algunos pensamientos:

En primer lugar, desde un punto de vista de programación, la acción es la respuesta "perfil". Es decir, codifíquelo de una manera más clara primero. Luego, perfila tu ejecución para ver si la optimización vale la pena. En mi humilde opinión, probablemente sea más claro usar una biblioteca matricial para estar más cerca de la matemática original.

Desde una vista matemática: Tengo un poco de dudas sobre la fórmula que proporciona para el caso multivariable. No me parece correcto. La expresión Z debe ser una forma cuadrática, y su Z no lo es. A menos que me esté perdiendo algo.

Aquí hay una opción que no mencionó, pero que podría tener sentido. Especialmente si va a evaluar el PDF varias veces para una sola distribución. Comience por calcular la base del componente principal de su distribución. Es decir, una base propia para Σ. Las direcciones de los componentes principales son ortogonales. En la base de componentes principales, las covarianzas cruzadas son todas 0, por lo que el PDF tiene una forma simple. Cuando quiera evaluar, cambie la base de la entrada en la base del componente principal y luego realice el cálculo de PDF más simple sobre eso.

La idea es que se puede calcular el cambio de la matriz base y los componentes principales una vez por adelantado, y luego solo hacer una multiplicación de matriz única (el cambio de base) por evaluación, en lugar de las dos multiplicaciones de matriz necesarias para evalúe el (x-μ)' Σ (x-μ) en la base estándar.

Answer 2

me han implementado básicamente el exp -parte de la ecuación para el caso tridimensional en this question. Usé una biblioteca de visión por computadora llamada OpenCV al principio. Pero noté que la interfaz de C++ era muy lenta. Después probé la interfaz C, que era un poco más rápida. Finalmente, decidí ignorar la flexibilidad y la legibilidad, así que lo implementé sin bibliotecas y fue mucho más rápido.

Lo que trato de decir es esto: cuando el rendimiento es importante, debe considerar la implementación de casos especiales para el número de dimensiones más utilizado con la mínima sobrecarga posible. De lo contrario, elija mantenibilidad sobre velocidad.

Descargo de responsabilidad: No sé nada sobre la velocidad de Eigen o boost::multi_array (¿cuál es probablemente a lo que esta pregunta apunta realmente?).