Algoritmo de muestreo anidado

El algoritmo de muestreo anidado es un enfoque computacional para los problemas de estadística bayesiana de comparación de modelos y generación de muestras a partir de distribuciones posteriores. Fue desarrollado en 2004 por el físico John Skilling. ^{[ 1 ]}

El teorema de Bayes se puede utilizar para la selección de modelos , cuando se tiene un par de modelos en competencia.${\displaystyle M_{1}}$y${\displaystyle M_{2}}$para datos${\displaystyle D}$Una de estas afirmaciones puede ser cierta (aunque se desconoce cuál), pero ambas no pueden serlo simultáneamente. La selección de modelos bayesianos proporciona un método para evaluar el factor de Bayes , que indica el mérito relativo de cada modelo.

La probabilidad posterior de${\displaystyle M_{1}}$se puede calcular como:

${\displaystyle {\begin{aligned}P(M_{1}\mid D)&={\frac {P(D\mid M_{1})P(M_{1})}{P(D)}}\\&={\frac {P(D\mid M_{1})P(M_{1})}{P(D\mid M_{1})P(M_{1})+P(D\mid M_{2})P(M_{2})}}\\&={\frac {1}{1+{\frac {P(D\mid M_{2})}{P(D\mid M_{1})}}{\frac {P(M_{2})}{P(M_{1})}}}}\end{aligned}}}$

Las probabilidades previas${\displaystyle M_{1}}$y${\displaystyle M_{2}}$ya se conocen, ya que son elegidos por el investigador con antelación. Sin embargo, el factor de Bayes restante${\displaystyle P(D\mid M_{2})/P(D\mid M_{1})}$no es tan fácil de evaluar, ya que en general requiere marginalizar parámetros molestos. En general,${\displaystyle M_{1}}$tiene un conjunto de parámetros que se pueden agrupar y llamar${\displaystyle \theta }$, y${\displaystyle M_{2}}$tiene su propio vector de parámetros que pueden ser de diferente dimensionalidad, pero aún así se denomina${\displaystyle \theta }$. La marginación para${\displaystyle M_{1}}$es

${\displaystyle P(D\mid M_{1})=\int d\theta \,P(D\mid \theta ,M_{1})P(\theta \mid M_{1})}$

y asimismo para${\displaystyle M_{2}}$Esta integral suele ser intratable analíticamente, y en estos casos es necesario emplear un algoritmo numérico para encontrar una aproximación. El algoritmo de muestreo anidado fue desarrollado por John Skilling específicamente para aproximar estas integrales de marginalización, y tiene el beneficio adicional de generar muestras de la distribución posterior.${\displaystyle P(\theta \mid D,M_{1})}$. ^{[ 2 ]} Es una alternativa a los métodos de la literatura bayesiana ^{[ 3 ]} como el muestreo de puente y el muestreo de importancia defensiva.

Aquí se presenta una versión simple del algoritmo de muestreo anidado, seguida de una descripción de cómo calcula la densidad de probabilidad marginal.${\displaystyle Z=P(D\mid M)}$dónde${\displaystyle M}$es${\displaystyle M_{1}}$o${\displaystyle M_{2}}$:

Comience con${\displaystyle N}$agujas${\displaystyle \theta _{1},\ldots ,\theta _{N}}$muestreado de anterior. para${\displaystyle i=1}$a${\displaystyle j}$hacer % El número de iteraciones j se elige por estimación. ${\displaystyle L_{i}:=\min(}$valores de probabilidad actuales de los puntos${\displaystyle )}$; ${\displaystyle X_{i}:=\exp(-i/N);}$${\displaystyle w_{i}:=X_{i-1}-X_{i}}$${\displaystyle Z:=Z+L_{i}\cdot w_{i};}$ Guarda el punto con menor probabilidad como punto de muestra con peso${\displaystyle w_{i}}$. Actualizar el punto con menor probabilidad mediante muestreo de la distribución previa restringida a probabilidades superiores${\displaystyle L_{i}}$, por ejemplo con el método de Monte Carlo de cadena de Markov . fin de retorno${\displaystyle Z}$;

En cada iteración,${\displaystyle X_{i}}$es una estimación de la cantidad de masa previa cubierta por el hipervolumen en el espacio de parámetros de todos los puntos con probabilidad mayor que${\displaystyle \theta _{i}}$El factor de ponderación${\displaystyle w_{i}}$es una estimación de la cantidad de masa previa que se encuentra entre dos hipersuperficies anidadas.${\displaystyle \{\theta \mid P(D\mid \theta ,M)=P(D\mid \theta _{i-1},M)\}}$y${\displaystyle \{\theta \mid P(D\mid \theta ,M)=P(D\mid \theta _{i},M)\}}$. El paso de actualización${\displaystyle Z:=Z+L_{i}w_{i}}$calcula la suma sobre${\displaystyle i}$de${\displaystyle L_{i}w_{i}}$para aproximar numéricamente la integral

${\displaystyle {\begin{aligned}P(D\mid M)&=\int P(D\mid \theta ,M)P(\theta \mid M)\,d\theta \\&=\int P(D\mid \theta ,M)\,dP(\theta \mid M)\end{aligned}}}$

En el límite${\displaystyle j\to \infty }$, este estimador tiene un sesgo positivo de orden${\displaystyle 1/N}$^{[ 4 ]} que se puede eliminar usando${\displaystyle (1-1/N)}$en lugar de la${\displaystyle \exp(-1/N)}$en el algoritmo anterior.

La idea es subdividir el rango de${\displaystyle f(\theta )=P(D\mid \theta ,M)}$y estimar, para cada intervalo${\displaystyle [f(\theta _{i-1}),f(\theta _{i})]}$, qué tan probable es a priori que un elemento elegido al azar${\displaystyle \theta }$se correspondería con este intervalo. Esto puede considerarse como una forma bayesiana de implementar numéricamente la integración de Lebesgue . ^{[ 5 ] [ 6 ]}

El punto con menor probabilidad se puede actualizar con algunos pasos de cadena de Markov Monte Carlo de acuerdo con el prior, aceptando solo pasos que mantengan la probabilidad por encima${\displaystyle L_{i}}$El procedimiento original descrito por Skilling (presentado anteriormente en pseudocódigo ) no especifica qué algoritmo concreto debe utilizarse para elegir nuevos puntos con mayor probabilidad, pero se han desarrollado varios algoritmos. ^{[ 7 ]}

Los propios ejemplos de código de Skilling (como uno en Sivia y Skilling (2006), ^{[ 8 ]} disponible en el sitio web de Skilling ) eligen un punto existente aleatorio y seleccionan un punto cercano elegido por una distancia aleatoria del punto existente; si la probabilidad es mejor, entonces el punto es aceptado, de lo contrario es rechazado y el proceso se repite. Posteriormente, se han desarrollado varios algoritmos MCMC adaptados para el muestreo anidado, incluido el muestreo por rebanadas, ^{[ 5 ]} popularizado por PolyChord, y el Monte Carlo hamiltoniano restringido . ^{[ 9 ]}

Una línea alternativa de algoritmos se basa en el muestreo por rechazo . Mukherjee et al. (2006) ^{[ 10 ]} encontraron tasas de aceptación más altas al seleccionar puntos aleatoriamente dentro de un elipsoide dibujado alrededor de los puntos existentes; esta idea fue refinada por el algoritmo MultiNest ^{[ 11 ]} , que maneja mejor las distribuciones posteriores multimodales utilizando múltiples elipsoides construidos a partir del agrupamiento de los puntos vivos. Los métodos de rechazo pueden ser eficientes hasta 20-30 dimensiones. ^{[ 7 ]}

Existen ejemplos de implementación que demuestran el algoritmo de muestreo anidado, disponibles para su descarga pública y escritos en varios lenguajes de programación .

• En el sitio web de John Skilling se pueden encontrar ejemplos sencillos en C , R o Python .
• Una versión en Haskell de los códigos simples anteriores está disponible en Hackage .
• En el sitio web de Bojan Nikolic se describe un ejemplo en R diseñado originalmente para ajustar espectros , que también está disponible en GitHub .
• Un NestedSampler forma parte de la caja de herramientas de Python BayesicFitting ^{[ 12 ]} para el ajuste de modelos genéricos y el cálculo de evidencia. Está disponible en GitHub .
• Existe una implementación en C++ , llamada DIAMONDS, en GitHub .
• En GitHub se encuentra disponible un ejemplo paralelo en Python altamente modular para aplicaciones en física estadística y física de la materia condensada .
• pymatnest es un paquete diseñado para explorar el panorama energético de diferentes materiales, calcular variables termodinámicas a temperaturas arbitrarias y localizar transiciones de fase; está disponible en GitHub.
• El paquete de software MultiNest es capaz de realizar muestreo anidado en distribuciones posteriores multimodales. ^{[ 11 ] [ 13 ]} Tiene interfaces para entradas en C++, Fortran y Python, y está disponible en GitHub .
• PolyChord es otro paquete de software de muestreo anidado disponible en GitHub . La eficiencia computacional de PolyChord escala mejor con un aumento en el número de parámetros que MultiNest, lo que significa que PolyChord puede ser más eficiente para problemas de alta dimensión. ^{[ 14 ]} Tiene interfaces para funciones de verosimilitud escritas en Python, Fortran, C o C++. PolyChord se puede usar junto con Cobaya, ^{[ 15 ]} un código basado en Python para el análisis bayesiano de modelos físicos jerárquicos. Cobaya facilita la exploración de posteriores usando varios muestreadores de Monte Carlo, permite la maximización y el reponderado de importancia de las muestras, e incluye interfaces para códigos de teoría cosmológica y verosimilitudes.
• NestedSamplers.jl, un paquete de Julia para implementar algoritmos de muestreo anidado elipsoidal simple y múltiple, está disponible en GitHub .
• Korali es un marco de alto rendimiento para la cuantificación de la incertidumbre, la optimización y el aprendizaje profundo por refuerzo, que también implementa el muestreo anidado.
• El paquete de software UltraNest implementa un algoritmo de muestreo dinámico anidado multielipsoidal generalizado y rápido, compatible con MPI . El usuario también puede elegir algoritmos de muestreo por secciones. Escrito en Python, cuenta con interfaces para Python, C, Fortran, R y Julia, y está disponible en GitHub .

Desde que se propuso el muestreo anidado en 2004, se ha utilizado en muchas áreas científicas, ^{[ 6 ]} en particular en astronomía . Un artículo sugirió usar el muestreo anidado para la selección de modelos cosmológicos y la detección de objetos, ya que "combina de forma única precisión, aplicabilidad general y viabilidad computacional". ^{[ 10 ]} Se ha sugerido un refinamiento del algoritmo para manejar posteriores multimodales como un medio para detectar objetos astronómicos en conjuntos de datos existentes. ^{[ 13 ]} Otras aplicaciones del muestreo anidado se encuentran en el campo de la actualización de elementos finitos, donde el algoritmo se utiliza para elegir un modelo de elementos finitos óptimo , y esto se aplicó a la dinámica estructural . ^{[ 16 ]} Este método de muestreo también se ha utilizado en el campo del modelado de materiales. Se puede utilizar para aprender la función de partición de la mecánica estadística y derivar propiedades termodinámicas . ^{[ 17 ]}

Se han desarrollado diagnósticos específicos para el muestreo anidado con el fin de verificar que una ejecución de muestreo anidado se esté realizando correctamente. Esto incluye una prueba U que comprueba que el rango de la probabilidad del punto de reemplazo se distribuye uniformemente entre los puntos activos, ^{[ 18 ] [ 7 ]} la distancia de salto de la cadena de Markov Monte Carlo, ^{[ 19 ]} la comparación de la consistencia de varias ejecuciones de muestreo anidado independientes, ^{[ 20 ] [ 21 ]} incluyendo ejecuciones repetidas con un mayor número de pasos MCMC. El cálculo también puede comprobarse con técnicas de aplicación genérica, como la calibración basada en simulación. ^{[ 22 ]}

El muestreo anidado dinámico es una generalización del algoritmo de muestreo anidado en el que el número de muestras tomadas en diferentes regiones del espacio de parámetros se ajusta dinámicamente para maximizar la precisión del cálculo. ^{[ 23 ]} Esto puede conducir a mejoras en la precisión y la eficiencia computacional en comparación con el algoritmo de muestreo anidado original, en el que la asignación de muestras no se puede cambiar y, a menudo, se toman muchas muestras en regiones que tienen poco efecto en la precisión del cálculo.

Los paquetes de software de muestreo dinámico anidado disponibles públicamente incluyen:

• dynesty : una implementación en Python de muestreo anidado dinámico que se puede descargar desde GitHub . ^{[ 24 ]}
• dyPolyChord: un paquete de software que se puede utilizar con Python, C++ y Fortran para distribuciones de probabilidad y a priori. ^{[ 25 ]} dyPolyChord está disponible en GitHub .
• UltraNest (ver arriba).

El muestreo dinámico anidado se ha aplicado a una variedad de problemas científicos, incluyendo el análisis de ondas gravitacionales, ^{[ 26 ]} el mapeo de distancias en el espacio ^{[ 27 ]} y la detección de exoplanetas. ^{[ 28 ]}

• Comparación de modelos bayesianos
• Lista de algoritmos