Al aproximar una posterior utilizando MCMC, ¿por qué no guardamos las probabilidades posteriores pero utilizamos las frecuencias de los valores de los parámetros después?

Question

Actualmente estoy estimando los parámetros de un modelo definido por varias ecuaciones diferenciales ordinarias (EDO). Lo intento con un enfoque bayesiano aproximando la distribución posterior de los parámetros dados algunos datos utilizando Markov Chain Monte Carlo (MCMC).

Un muestreador MCMC genera una cadena de valores de parámetros en la que utiliza la probabilidad posterior (no normalizada) de un determinado valor de parámetro para decidir (estocásticamente) si añadirá ese valor a la cadena o volverá a añadir el valor anterior. Pero, parece ser la práctica que las probabilidades posteriores reales no necesitan ser guardadas, más bien se genera un histograma n-dimensional de los valores de los parámetros resultantes y las estadísticas de resumen como las regiones de mayor densidad (HDRs) de una ditribución posterior de los parámetros se calcula a partir de este histograma. Al menos eso es lo que creo que aprendí de Libro tutorial de Kruschkes sobre la inferencia bayesiana .

Mi pregunta: ¿No sería más sencillo guardar las probabilidades posteriores de los valores de los parámetros muestreados junto con estos y aproximar la distribución posterior a partir de estos valores y no de las frecuencias de los valores de los parámetros en la cadena MCMC? El problema de la fase de burn-in no surgiría ya que el muestreador seguiría muestreando inicialmente regiones de baja probabilidad con más frecuencia de la que "merecerían" por sus probabilidades posteriores pero ya no tendría el problema de dar valores de probabilidad indebidamente altos a estas.

Answer 1

Esta es una pregunta interesante, con diferentes cuestiones:

1. Los algoritmos MCMC no siempre reciclan el cálculo de la densidad posterior en todos los valores propuestos, pero algunas técnicas de reducción de la varianza como la Rao-Blackwellización sí lo hacen. Por ejemplo, en un Artículo de Biometrika de 1996 con George Casella, proponemos utilizar todos los valores simulados, $\theta_i$ $(i=1,\ldots,T)$ aceptada o no, mediante la introducción de pesos $\omega_i$ que convierten la media $$\sum_{i=1}^T \omega_ih(\theta_i)\big/\sum_{i=1}^T \omega_i$$ en un estimador casi insesgado. (El casi debido a la normalización por la suma de los pesos).
2. El MCMC se utiliza a menudo en problemas de gran dimensión (de parámetros). Proponer una aproximación a la totalidad de la posterioridad basada en los valores de la densidad observada en algunos valores de los parámetros es todo un reto, incluida la cuestión de la constante normalizadora mencionada en la respuesta y los comentarios de Tim. Se puede imaginar una aproximación que sea una mezcla de estimación kernel no paramétrica (como en, por ejemplo krigging ) y la regresión, pero los expertos con los que hablé de esta solución [hace unos años] eran bastante escépticos. El problema es que el estimador resultante sigue siendo no paramétrico y, por lo tanto, "disfruta" de velocidades de convergencia no paramétricas que son más lentas que las velocidades de convergencia de Monte Carlo, tanto peores cuanto mayor sea la dimensión.
3. Otro uso potencial de la disponibilidad de los valores posteriores $\pi(\theta|\mathcal{D})$ es ponderar cada valor simulado por su posterior asociado, como en $$\frac{1}{T}\sum_{t=1}^T h(\theta_t) \pi(\theta_t|\mathcal{D})$$ Desgraciadamente, esto crea un sesgo, ya que los valores simulados ya se han simulado a partir de la posterioridad: $$\mathbb{E}[h(\theta_t) \pi(\theta_t|\mathcal{D})]=\int h(\theta) h(\theta_t) \pi(\theta_t|\mathcal{D})^2\text{d}\theta$$ Incluso sin un problema de normalización, esas simulaciones deberían dirigirse a $\pi(\theta|\mathcal{D})^{1/2}$ y utilizar un peso proporcional a $\pi(\theta|\mathcal{D})^{1/2}$ pero no conozco resultados que defiendan este cambio de objetivo. Como se menciona en los comentarios, esto está relacionado con templado en el sentido de que todas las simulaciones producidas dentro de un ciclo de templado simulado pueden ser recicladas para fines de Monte Carlo (integración) de esta manera. Sin embargo, un problema numérico es manejar varias funciones de importancia de la forma $\pi(\theta)^{1/T}$ con las constantes de normalización que faltan.

Answer 2

Como has notado correctamente, las probabilidades con las que estamos tratando son sin normalizar . Básicamente, utilizamos MCMC para calcular el factor de normalización en el teorema de Bayes. No podemos utilizar las probabilidades porque no están normalizados. El procedimiento que sugieres: guardar las probabilidades no normalizadas y luego dividirlas por su suma es incorrecto.

Permítanme mostrárselo con un ejemplo. Imagina que utilizas Monte Carlo para extraer diez valores de la distribución Bernoulli parametrizada por $p=0.9$ son los siguientes:

1 0 1 1 1 1 1 1 1 1

también tienes las probabilidades correspondientes:

0.9 0.1 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

En este caso las probabilidades están normalizadas, pero dividirlas por su suma (que por axiomas de la probabilidad es igual a la unidad) no debería cambiar nada. Desgraciadamente, utilizando tu procedimiento hace cambiar los resultados a:

> f/sum(f)
 [1] 0.10975610 0.01219512 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610

¿Por qué? La respuesta es simple, en su muestra cada "probabilidad" guardada f aparece con probabilidad f ¡así que estás ponderando las probabilidades por sí mismas!