Cálculo de la probabilidad marginal de muestras de MCMC

Esta es una pregunta recurrente (vea esta publicación , esta publicación y esta publicación ), pero tengo un giro diferente.

Supongamos que tengo un montón de muestras de una muestra genérica de MCMC. Para cada muestra $\theta$ , conozco el valor del log de verosimilitud $\log f(\textbf{x} | \theta)$ y del log anterior $\log f(\theta)$ . Si ayuda, también sé el valor de la probabilidad de registro por punto de datos, $\log f(x_i | \theta)$ (esta información ayuda con ciertos métodos, como WAIC y PSIS-LOO).

Quiero obtener una estimación (cruda) de la probabilidad marginal, solo con las muestras que tengo, y posiblemente algunas otras evaluaciones de función (pero sin volver a ejecutar un MCMC ad hoc ).

En primer lugar, despejemos la tabla. Todos sabemos que el estimador armónico es el peor estimador de la historia . Vamonos. Si está haciendo un muestreo de Gibbs con anteriores y posteriores en forma cerrada, puede usar el método de Chib ; pero no estoy seguro de cómo generalizar fuera de esos casos. También hay métodos que requieren que modifique el procedimiento de muestreo (como a través de posteriores templados ), pero no estoy interesado en eso aquí.

El enfoque en el que estoy pensando consiste en aproximar la distribución subyacente con una forma paramétrica (o no paramétrica) $g(\theta)$ , y luego descubrir la constante de normalización $Z$ como un problema de optimización 1-D (es decir, la $Z$ que minimiza algún error entre $Z g(\theta)$ y $f(\textbf{x}|\theta) f(\theta)$ , evaluada sobre las muestras). En el caso más simple, supongamos que la parte posterior es aproximadamente multivariada normal, puedo ajustar $g(\theta)$como normal multivariante y obtener algo similar a una aproximación de Laplace (es posible que desee utilizar algunas evaluaciones de funciones adicionales para refinar la posición del modo). Sin embargo, podría usar como $g(\theta)$ una familia más flexible, como una mezcla variacional de distribuciones multivariadas $t$.

Aprecio que este método solo funciona si $Z g(\theta)$ es una aproximación razonable a $f(\textbf{x}|\theta) f(\theta)$ , pero ¿hay alguna razón o una advertencia de por qué sería muy imprudente hacerlo? ¿Alguna lectura que recomendarías?

El enfoque totalmente no paramétrico utiliza alguna familia no paramétrica, como un proceso gaussiano (GP), para aproximar $\log f(\textbf{x}|\theta) + \log f(\theta)$ (o alguna otra transformación no lineal de la misma, como la raíz cuadrada), y bayesiana cuadratura para integrarse implícitamente sobre el objetivo subyacente (ver aquí y aquí ). Este parece ser un enfoque alternativo interesante, pero análogo en espíritu (también, tenga en cuenta que los médicos de familia serían difíciles de manejar en mi caso).

Desafortunadamente, la extensión de Chib y Jeliazkov (2001) se vuelve rápidamente costosa o muy variable, razón por la cual no se usa mucho fuera de los casos de muestreo de Gibbs.

Si bien hay muchas formas y enfoques para el problema de la estimación constante de normalización (como lo ilustran las charlas bastante diversas en el taller Estimación constante que realizamos la semana pasada en la Universidad de Warwick, diapositivas disponibles allí ), algunas soluciones explotan directamente la salida de MCMC .$\mathfrak{Z}$

1. Como mencionó, el estimador de la media armónica de Newton y Raftery (1994) es casi siempre pobre por tener una varianza infinita. Sin embargo, hay formas de evitar la maldición de varianza infinita utilizando en su lugar un objetivo de soporte finito en la identidad media armónica eligiendoαcomo el indicador de una región HPD para la posterior. Esto asegura una variación finita al eliminar las colas en la media armónica. (Los detalles se encuentran enun artículo que escribí con Darren Wraithy en uncapítulo sobre la normalización de las constantesescritas con Jean-Michel Marin.) En resumen, el método recicla la salida MCMCθ1,...,θMidentificando elβ( 20% dice) valores más grandes del objetivoπ(θ)f(x|θ)y creando

   $$\int \dfrac{\alpha(\theta)}{\pi(\theta)f(x|\theta)}\text{d}\pi(\theta|x)=\frac{1}{\mathfrak{Z}}$$ $\alpha$$\theta_1,\ldots,\theta_M$$\beta$$\pi(\theta)f(x|\theta)$$\alpha$ como un uniforme sobre la unión de las bolas centradas en esas simulaciones de mayor densidad (HPD) $\theta^0_i$ y con el radio , es decir, la estimación de la constante de normalización Z viene dada por Z - 1 = $\rho$$\mathfrak{Z}$ sides la dimensión deθ(se aplican correcciones para las bolas que se cruzan) y siρes lo suficientemente pequeño como para que las bolas nunca se crucen (lo que significa que, en el mejor de los casos, solo un indicador en las bolas es diferente de cero). La explicación para eldenominadorαM2es que esta es una doble suma detérminosβM2: $$\hat{\mathfrak{Z}}^{-1}=\underbrace{\frac{1}{\beta M^2}\sum_{m=1}^M}_{\text{double sum over}\\\beta M\text{ ball centres }\theta_i^0\\\text{and $M$ simulations } \theta_m} \underbrace{\mathbb{I}_{(0,\rho)}(\min_i||\theta_m-\theta^0_i||)\{\pi(\theta_m)f(x|\theta_m)\}^{-1}\big/\overbrace{\pi^{d/2}\rho^d\Gamma(d/2+1)^{-1}}^{\text{volume of ball with radius $\rho$}}}_{\dfrac{\beta M\alpha(\theta_m)}{\pi(\theta_m)f(x|\theta_m)}}$$ $d$$\theta$$\rho$$\alpha M^2$$\beta M^2$ $$\frac{1}{\beta M}\sum_{i=1}^{\beta M} \underbrace{\frac{1}{M}\sum_{m=1}^M {\cal U}(\theta_i^0,\rho)(\theta_m)}_{\text{same as with $\min$}} \times \frac{1}{\pi(\theta_m)f(x|\theta_m)}$$ with each term in $\theta_m$ integrating to ${\mathfrak{Z}}^{-1}$.

2. Otro enfoque es convertir la constante de normalización en un parámetro. Esto suena como una herejía estadística, pero el artículo de Guttmann y Hyvärinen (2012) me convenció de lo contrario. Sin entrar demasiado en detalles, la idea clara es convertir la probabilidad de registro observada n ∑ i = 1 f ( x i | θ ) - n log ∫ exp f ( x | θ ) d x en una probabilidad de registro conjunta n ∑ i = 1 [ $\mathfrak{Z}$

   $$\sum_{i=1}^n f(x_i|\theta) - n \log \int \exp f(x|\theta) \text{d}x$$ que es la probabilidad logarítmica de un proceso de punto de Poisson con función de intensidad exp { f ( x | θ ) + ν + log n $$\sum_{i=1}^n[f(x_i|\theta)+\nu]-n\int\exp[f(x|\theta)+\nu]\text{d}x$$ $$\exp\{ f(x|\theta) + \nu +\log n\}$$ Este es un modelo alternativo en el sentido de que la probabilidad original no aparece como marginal de lo anterior. Solo los modos coinciden, con el modo condicional en ν que proporciona la constante de normalización. En la práctica, la probabilidad del proceso de Poisson anterior no está disponible y Guttmann y Hyvärinen (2012) ofrecen una aproximación por medio de una regresión logística. Para conectarse aún mejor con su pregunta, la estimación de Geyer es un MLE, por lo tanto, la solución a un problema de maximización.
3. $\pi(\theta|x)$$\pi(\theta|x)$, $g(\theta)$y ejecutar una regresión logística en el índice de la distribución detrás de los datos (1 para $\pi(\theta|x)$ y 0 para $g(\theta)$) Con los regresores siendo los valores de ambas densidades, normalizadas o no. Esto está directamente relacionado con el muestreo en puente de Gelman y Meng (1997), que también recicla muestras de diferentes objetivos. Y versiones posteriores, como el MLE de Meng.
4. Un enfoque diferente que obliga a ejecutar un muestreador MCMC específico es el muestreo anidado de Skilling . Si bien yo [y otros] tengo algunas reservas sobre la eficiencia del método, es bastante popular en astrostatistics y cosmología, con software disponible como multinest .
5. Una última solución [potencial, si no siempre posible] es explotar la representación Savage-Dickey del factor Bayes en el caso de una hipótesis nula incrustada. Si el nulo escribe como$H_0: \theta=\theta_0$ sobre un parámetro de interés y si $\xi$ es la parte restante [molestia] del parámetro del modelo, suponiendo un previo de la forma $\pi_1(\theta)\pi_2(\xi)$, el factor Bayes de $H_0$ contra la alternativa escribe como $$\mathfrak{B}_{01}(x)=\dfrac{\pi^\theta(\theta_0|x)}{\pi_1(\theta_0)}$$ dónde $\pi^\theta(\theta_0|x)$ denota la densidad posterior marginal de $\theta$ en el valor específico $\theta_0$. En caso de que la densidad marginal debajo de la nula$H_0: \theta=\theta_0$ $$m_0(x)=\int_\Xi f(x|\theta_0,\xi)\pi_2(\xi)\text{d}\xi$$ está disponible en forma cerrada, se puede derivar la densidad marginal para el modelo sin restricciones $$m_a(x)=\int_{\Theta\times\Xi} f(x|\theta,\xi)\pi_1(\theta)\pi_2(\xi)\text{d}\theta\text{d}\xi$$ del factor Bayes. (Esta representación de Savage-Dickey se basa en versiones específicas de tres densidades diferentes y, por lo tanto, está llena de peligros, sin mencionar siquiera el desafío computacional de producir el posterior marginal).

[Aquí hay un conjunto de diapositivas que escribí sobre la estimación de constantes de normalización para un taller de NIPS en diciembre pasado].