¿Existe una muestra de Monte Carlo / MCMC implementada que pueda tratar con máximos locales aislados de distribución posterior?

10

Actualmente estoy usando un enfoque bayesiano para estimar los parámetros de un modelo que consta de varias EDO. Como tengo 15 parámetros para estimar, mi espacio de muestreo es de 15 dimensiones y mi búsqueda de distribución posterior parece tener muchos máximos locales que están muy aislados por grandes regiones de muy baja probabilidad.

Esto lleva a problemas de mezcla de mis cadenas Monte Carlo, ya que es muy poco probable que una cadena "salte" de un máximo local y accidentalmente golpee a uno de los otros máximos.

Parece que hay mucha investigación en esta área, ya que es fácil encontrar documentos que aborden este problema (ver más abajo), pero encontrar una implementación real es difícil. Solo encontré paquetes relacionados con la dinámica molecular, pero no con la inferencia bayesiana. ¿Existen implementaciones de muestreadores MC (MC) que puedan manejar máximos locales aislados por ahí?

Me veo obligado a trabajar con Matlab ya que eso es lo que está escrito en mi modelo ODE, por lo que las propuestas con respecto a Matlab son bienvenidas ;-). Sin embargo, si hay una "aplicación asesina" en algún otro idioma, tal vez pueda convencer a mi PI para que cambie ;-).

Actualmente estoy trabajando con una muestra de Monte Carlo de rechazo tardío / adaptativo escrita por Haario, Laine et al. , y esa es también la única muestra que pude encontrar hasta ahora que es más sofisticada que el algoritmo estándar de Metropolis-Hastings


Los enfoques notables parecen ser:

EDITAR Actualizado el 2017-Mar-07 con lo que he aprendido mientras tanto

Múltiples cadenas similares con diferentes puntos de partida.

Adaptación entre cadenas. Utilice la matriz de covarianza empírica de las muestras agrupadas generadas por múltiples cadenas independientes para actualizar las matrices de covarianza de las distribuciones de propuestas de la cadena. (1)

Cadenas múltiples con diferente temple

Templado: algún tipo de "temperatura" parece modificar el paisaje posterior, haciendo que la mezcla de cadenas sea más probable. (Todavía no me he sumergido mucho en esto) (1) El propósito del temple es aplanar el paisaje de probabilidad (de alta dimensión) formado por la distribución de probabilidad posterior. Se logra comúnmente tomando la probabilidad posterior a la potencia de , donde el paisaje posterior se aplana para (3, p.298). Esto significa que, en lugar de calcular la probabilidad posterior de un estado , dados los datos se calcula la probabilidad posterior moderada1/TT>1p(θD)θD

p(θD)1/T(p(Dθ)p(θ))1/T

Cuanto más alto se elige , los picos más planos y más amplios en el paisaje de probabilidad se vuelven. Por lo tanto, los valores más altos de conducen a una mayor probabilidad de que el muestreador cambie de un máximo local a otro. Sin embargo, no es la distribución posterior buscada si . Por lo tanto, la cadena de muestras de esa distribución debe usarse para permitir el muestreo de después.TTp(θD)1/TT1p(θD)

Las muestras de la distribución posterior original sin templar, las muestras dadas de una versión templada de esa distribución se pueden obtener por varios métodos:

  • Metropolis acoplado MCMC Run múltiples cadenas simultáneamente, cada uno con un valor diferente pero constante para . Cambia los estados de dos cadenas probabilísticamente. Utilice solo las muestras de la cadena con para las estimaciones posteriores; las otras cadenas solo se aseguran de que se muestren todos los picos. Árbitro. (4) tiene un algoritmo paralelo y cita un artículo de conferencia y un libro de texto para la idea (5,6)TT=1

  • MCMC de mundo pequeño. La muestra cambia entre dos propuestas. La mayoría de las veces se usa una distribución de propuesta con una pequeña varianza, rara vez se usa una propuesta con una gran varianza. La elección entre estas dos propuestas es estocástica. Las propuestas con gran varianza también pueden extraerse de otra cadena que solo realiza saltos muy grandes, muestreando la mayor cantidad posible del espacio muestral de manera aproximada. (2,7)

Hamiltoniano Monte Carlo (HMC)

No sé mucho sobre eso, pero la muestra No-U-Turn (NUTS) de JAGS parece usarlo. Ver ref. (8) Alex Rogozhnikov ha creado un tutorial visual sobre el tema.


Referencias

(1) Craiu et al., 2009: Aprenda de su vecino: MCMC adaptativo regional y de cadena paralela. J Am Stat Assoc 104: 488, págs. 1454-1466. http://www.jstor.org/stable/40592353

(2) Guam et al., 2012: Small World MCMC with temple: Ergocity and espectral gap. https://arxiv.org/abs/1211.4675 ( solo en arXiv )

(3): Brooks y col. (2011) Manual de Markov Chain Monte Carlo. Prensa CRC.

(4): Altekar et al. (2004): Metrópolis paralela acopló la cadena de Markov Monte Carlo para la inferencia filogenética bayesiana. Bioinformática 20 (3) 2004, pp. 407–415, http://dx.doi.org/10.1093/bioinformatics/btg427

(5): Geyer CJ (1991) Markov chain Monte Carlo máxima verosimilitud. En: Keramidas (ed.), Ciencias de la computación y estadística: Actas del 23º Simposio sobre la interfaz . Interface Foundation, Fairfax Station, págs. 156–163.

(6): Gilks ​​WR y Roberts GO (1996). Estrategias para mejorar MCMC. En: Gilks ​​WR, Richardson S y Spiegelhalter (eds) Markov encadena Monte Carlo en la práctica . Chapman & Hall, pág. 89-114.

(7): Guan Y, et al. Markov Chain Monte Carlo en pequeños mundos. Estadísticas e informática (2006) 16 (2), págs. 193-202. http://dx.doi.org/10.1007/s11222-006-6966-6

(8): Hoffmann M y Gelman A (2014): The No-U-Turn Sampler: Adaptively Lengths Path Paths in Hamiltonian Monte Carlo. Journal of Machine Learning Research , 15, págs. 1351-1381. https://arxiv.org/abs/1111.4246

akraf
fuente

Respuestas:

1

Ninguna de las estrategias anteriores es particularmente adecuada para múltiples óptimas.

Una mejor opción son MCMC de evolución diferencial y MCMC derivados como DREAM. Estos algoritmos funcionan con varias cadenas MCMC que se mezclan para generar propuestas. Si tiene al menos una cadena en cada óptimo, podrán saltar de manera eficiente entre los óptimos. Una implementación en R está disponible aquí https://cran.r-project.org/web/packages/BayesianTools/index.html

Florian Hartig
fuente