Muestreo de Gibbs

Algoritmo de Monte Carlo

En estadística , el muestreo de Gibbs o un muestreador de Gibbs es un algoritmo de Monte Carlo de cadena de Markov (MCMC) para muestrear a partir de una distribución de probabilidad multivariante especificada cuando el muestreo directo a partir de la distribución conjunta es difícil, pero el muestreo a partir de la distribución condicional es más práctico. Esta secuencia se puede utilizar para aproximar la distribución conjunta (por ejemplo, para generar un histograma de la distribución); para aproximar la distribución marginal de una de las variables, o algún subconjunto de las variables (por ejemplo, los parámetros desconocidos o variables latentes ); o para calcular una integral (como el valor esperado de una de las variables). Normalmente, algunas de las variables corresponden a observaciones cuyos valores se conocen y, por lo tanto, no es necesario muestrearlas.

El muestreo de Gibbs se utiliza habitualmente como medio de inferencia estadística , especialmente inferencia bayesiana . Es un algoritmo aleatorio (es decir, un algoritmo que hace uso de números aleatorios ) y es una alternativa a los algoritmos deterministas para la inferencia estadística, como el algoritmo de expectativa-maximización (EM).

Al igual que con otros algoritmos MCMC, el muestreo de Gibbs genera una cadena de Markov de muestras, cada una de las cuales está correlacionada con muestras cercanas. Como resultado, se debe tener cuidado si se desean muestras independientes. Por lo general, las muestras del comienzo de la cadena (el período de rodaje ) pueden no representar con precisión la distribución deseada y, por lo general, se descartan.

Introducción

El muestreo de Gibbs recibe su nombre del físico Josiah Willard Gibbs , en referencia a una analogía entre el algoritmo de muestreo y la física estadística . El algoritmo fue descrito por los hermanos Stuart y Donald Geman en 1984, unas ocho décadas después de la muerte de Gibbs, ^[1] y se popularizó en la comunidad estadística para calcular la distribución de probabilidad marginal, especialmente la distribución posterior. ^[2]

En su versión básica, el muestreo de Gibbs es un caso especial del algoritmo Metropolis-Hastings . Sin embargo, en sus versiones extendidas (ver más abajo), puede considerarse un marco general para el muestreo de un conjunto grande de variables mediante el muestreo de cada variable (o en algunos casos, de cada grupo de variables) por turno, y puede incorporar el algoritmo Metropolis-Hastings (o métodos como el muestreo por sectores ) para implementar uno o más de los pasos de muestreo.

El muestreo de Gibbs es aplicable cuando la distribución conjunta no se conoce explícitamente o es difícil muestrearla directamente, pero la distribución condicional de cada variable es conocida y es fácil (o al menos, más fácil) muestrearla. El algoritmo de muestreo de Gibbs genera una instancia de la distribución de cada variable a su vez, condicional a los valores actuales de las otras variables. Se puede demostrar que la secuencia de muestras constituye una cadena de Markov , y la distribución estacionaria de esa cadena de Markov es simplemente la distribución conjunta buscada. ^[3]

El muestreo de Gibbs está particularmente bien adaptado para muestrear la distribución posterior de una red bayesiana , ya que las redes bayesianas normalmente se especifican como una colección de distribuciones condicionales.

Implementación

El muestreo de Gibbs, en su forma básica, es un caso especial del algoritmo Metropolis-Hastings . El objetivo del muestreo de Gibbs es que, dada una distribución multivariante, es más sencillo tomar muestras de una distribución condicional que marginalizar mediante la integración sobre una distribución conjunta . Supongamos que queremos obtener muestras de de una distribución conjunta . Denotemos la muestra n por . Procedemos de la siguiente manera: ${\displaystyle \left.k\right.}$ ${\displaystyle \mathbf {X} =(x_{1},\dots ,x_{n})}$ ${\displaystyle p(x_{1},\dots ,x_{n})}$ ${\displaystyle i}$ ${\displaystyle \mathbf {X} ^{(i)}=\left(x_{1}^{(i)},\dots ,x_{n}^{(i)}\right)}$

1. Comenzamos con algún valor inicial . ${\displaystyle \mathbf {X} ^{(0)}}$
2. Queremos la siguiente muestra. Llamemos a esta siguiente muestra . Como es un vector, muestreamos cada componente del vector, , de la distribución de ese componente condicionada a todos los demás componentes muestreados hasta ahora. Pero hay un truco: condicionamos los componentes de hasta , y luego condicionamos los componentes de , comenzando desde hasta . Para lograr esto, muestreamos los componentes en orden, comenzando desde el primer componente. Más formalmente, para muestrear , lo actualizamos de acuerdo con la distribución especificada por . Usamos el valor que tenía el componente n en la n muestra, no la n muestra. ${\displaystyle \mathbf {X} ^{(i+1)}}$ ${\displaystyle \mathbf {X} ^{(i+1)}=\left(x_{1}^{(i+1)},x_{2}^{(i+1)},\dots ,x_{n}^{(i+1)}\right)}$ ${\displaystyle x_{j}^{(i+1)}}$ ${\displaystyle \mathbf {X} ^{(i+1)}}$ ${\displaystyle x_{j-1}^{(i+1)}}$ ${\displaystyle \mathbf {X} ^{(i)}}$ ${\displaystyle x_{j+1}^{(i)}}$ ${\displaystyle x_{n}^{(i)}}$ ${\displaystyle x_{j}^{(i+1)}}$ ${\displaystyle p\left(x_{j}^{(i+1)}|x_{1}^{(i+1)},\dots ,x_{j-1}^{(i+1)},x_{j+1}^{(i)},\dots ,x_{n}^{(i)}\right)}$ ${\displaystyle (j+1)}$ ${\displaystyle i}$ ${\displaystyle (i+1)}$
3. Repita los pasos anteriores varias veces. ${\displaystyle k}$

Propiedades

Si se realiza dicho muestreo, se cumplen estos importantes hechos:

• Las muestras aproximan la distribución conjunta de todas las variables.
• La distribución marginal de cualquier subconjunto de variables se puede aproximar simplemente considerando las muestras para ese subconjunto de variables, ignorando el resto.
• El valor esperado de cualquier variable se puede aproximar promediando todas las muestras.

Al realizar el muestreo:

• Los valores iniciales de las variables pueden determinarse aleatoriamente o mediante algún otro algoritmo como el de expectativa-maximización .
• En realidad no es necesario determinar un valor inicial para la primera variable muestreada.
• Es común ignorar una cierta cantidad de muestras al principio (el llamado período de quemado ) y luego considerar solo cada muestra al promediar los valores para calcular una expectativa. Por ejemplo, se pueden ignorar las primeras 1000 muestras y luego promediar cada muestra número 100, descartando el resto. La razón para esto es que (1) la distribución estacionaria de la cadena de Markov es la distribución conjunta deseada sobre las variables, pero puede llevar un tiempo alcanzar esa distribución estacionaria; (2) las muestras sucesivas no son independientes entre sí, sino que forman una cadena de Markov con cierta cantidad de correlación. A veces, se pueden usar algoritmos para determinar la cantidad de autocorrelación entre muestras y el valor de (el período entre muestras que se usan realmente) calculado a partir de esto, pero en la práctica hay una buena cantidad de " magia negra " involucrada. ${\displaystyle n}$ ${\displaystyle n}$
• El proceso de recocido simulado se utiliza a menudo para reducir el comportamiento de " caminata aleatoria " en la primera parte del proceso de muestreo (es decir, la tendencia a moverse lentamente por el espacio muestral, con una gran cantidad de autocorrelación entre muestras, en lugar de moverse rápidamente, como se desea). Otras técnicas que pueden reducir la autocorrelación son el muestreo de Gibbs colapsado , el muestreo de Gibbs bloqueado y la sobrerelajación ordenada ; consulte a continuación.

Relación entre distribución condicional y distribución conjunta

Además, la distribución condicional de una variable dadas todas las demás es proporcional a la distribución conjunta:

${\displaystyle p(x_{j}\mid x_{1},\dots ,x_{j-1},x_{j+1},\dots ,x_{n})={\frac {p(x_{1},\dots ,x_{n})}{p(x_{1},\dots ,x_{j-1},x_{j+1},\dots ,x_{n})}}\propto p(x_{1},\dots ,x_{n})}$

"Proporcional a" en este caso significa que el denominador no es una función de y, por lo tanto, es el mismo para todos los valores de ; forma parte de la constante de normalización para la distribución sobre . En la práctica, para determinar la naturaleza de la distribución condicional de un factor , es más fácil factorizar la distribución conjunta según las distribuciones condicionales individuales definidas por el modelo gráfico sobre las variables, ignorar todos los factores que no son funciones de (todos los cuales, junto con el denominador anterior, constituyen la constante de normalización) y luego restablecer la constante de normalización al final, según sea necesario. En la práctica, esto significa hacer una de tres cosas: ${\displaystyle x_{j}}$ ${\displaystyle x_{j}}$ ${\displaystyle x_{j}}$ ${\displaystyle x_{j}}$ ${\displaystyle x_{j}}$

1. Si la distribución es discreta, se calculan las probabilidades individuales de todos los valores posibles de y luego se suman para encontrar la constante de normalización. ${\displaystyle x_{j}}$
2. Si la distribución es continua y de forma conocida, también se conocerá la constante de normalización.
3. En otros casos, la constante de normalización generalmente se puede ignorar, ya que la mayoría de los métodos de muestreo no la requieren.

Inferencia

El muestreo de Gibbs se utiliza habitualmente para la inferencia estadística (por ejemplo, para determinar el mejor valor de un parámetro, como la cantidad de personas que probablemente compren en una tienda en particular en un día determinado, el candidato por el que un votante probablemente votará, etc.). La idea es que los datos observados se incorporen al proceso de muestreo creando variables separadas para cada dato observado y fijando las variables en cuestión a sus valores observados, en lugar de realizar un muestreo a partir de esas variables. La distribución de las variables restantes es entonces efectivamente una distribución posterior condicionada a los datos observados.

El valor más probable de un parámetro deseado (la moda ) podría entonces seleccionarse simplemente eligiendo el valor de muestra que ocurre más comúnmente; esto es esencialmente equivalente a la estimación máxima a posteriori de un parámetro. (Dado que los parámetros suelen ser continuos, a menudo es necesario "agrupar" los valores muestreados en uno de un número finito de rangos o "bins" para obtener una estimación significativa de la moda). Sin embargo, más comúnmente, se elige el valor esperado ( media o promedio) de los valores muestreados; este es un estimador de Bayes que aprovecha los datos adicionales sobre toda la distribución que están disponibles a partir del muestreo bayesiano, mientras que un algoritmo de maximización como la maximización de expectativas (EM) solo es capaz de devolver un único punto de la distribución. Por ejemplo, para una distribución unimodal, la media (valor esperado) suele ser similar a la moda (valor más común), pero si la distribución está sesgada en una dirección, la media se moverá en esa dirección, lo que explica efectivamente la masa de probabilidad adicional en esa dirección. (Si una distribución es multimodal, el valor esperado puede no devolver un punto significativo y cualquiera de los modos suele ser una mejor opción).

Aunque algunas de las variables corresponden típicamente a parámetros de interés, otras son variables poco interesantes ("molestas") que se introducen en el modelo para expresar adecuadamente las relaciones entre las variables. Aunque los valores muestreados representan la distribución conjunta de todas las variables, las variables molestas pueden simplemente ignorarse al calcular los valores esperados o las modas; esto es equivalente a marginar sobre las variables molestas. Cuando se desea un valor para múltiples variables, el valor esperado simplemente se calcula sobre cada variable por separado. (Sin embargo, al calcular la moda, todas las variables deben considerarse juntas).

El aprendizaje supervisado , el aprendizaje no supervisado y el aprendizaje semisupervisado (también conocido como aprendizaje con valores faltantes) se pueden manejar simplemente fijando los valores de todas las variables cuyos valores se conocen y tomando muestras del resto.

En el caso de los datos observados, habrá una variable para cada observación, en lugar de, por ejemplo, una variable correspondiente a la media o la varianza de la muestra de un conjunto de observaciones. De hecho, generalmente no habrá ninguna variable correspondiente a conceptos como "media de la muestra" o "varianza de la muestra". En cambio, en tal caso habrá variables que representen la media y la varianza verdaderas desconocidas, y la determinación de los valores de la muestra para estas variables resulta automáticamente del funcionamiento del muestreador de Gibbs.

Los modelos lineales generalizados (es decir, variaciones de la regresión lineal ) también pueden manejarse a veces mediante el muestreo de Gibbs. Por ejemplo, la regresión probit para determinar la probabilidad de una elección binaria dada (sí/no), con valores previos distribuidos normalmente colocados sobre los coeficientes de regresión, se puede implementar con el muestreo de Gibbs porque es posible agregar variables adicionales y aprovechar la conjugación . Sin embargo, la regresión logística no se puede manejar de esta manera. Una posibilidad es aproximar la función logística con una mezcla (típicamente 7-9) de distribuciones normales. Sin embargo, más comúnmente, se utiliza Metropolis-Hastings en lugar del muestreo de Gibbs.

Antecedentes matemáticos

Supongamos que se toma una muestra de una distribución que depende de un vector de parámetros de longitud , con distribución previa . Puede ser que sea muy grande y que la integración numérica para encontrar las densidades marginales de la sería computacionalmente costosa. Entonces, un método alternativo para calcular las densidades marginales es crear una cadena de Markov en el espacio repitiendo estos dos pasos: ${\displaystyle \left.X\right.}$ ${\displaystyle \theta \in \Theta \,\!}$ ${\displaystyle \left.d\right.}$ ${\displaystyle g(\theta _{1},\ldots ,\theta _{d})}$ ${\displaystyle \left.d\right.}$ ${\displaystyle \left.\theta _{i}\right.}$ ${\displaystyle \left.\Theta \right.}$

1. Elija un índice aleatorio ${\displaystyle 1\leq j\leq d}$
2. Elija un nuevo valor para según ${\displaystyle \left.\theta _{j}\right.}$ ${\displaystyle g(\theta _{1},\ldots ,\theta _{j-1},\,\cdot \,,\theta _{j+1},\ldots ,\theta _{d})}$

Estos pasos definen una cadena de Markov reversible con la distribución invariante deseada . Esto se puede demostrar de la siguiente manera. Defina si para todos y sea la probabilidad de un salto de a . Entonces, las probabilidades de transición son ${\displaystyle \left.g\right.}$ ${\displaystyle x\sim _{j}y}$ ${\displaystyle \left.x_{i}=y_{i}\right.}$ ${\displaystyle i\neq j}$ ${\displaystyle \left.p_{xy}\right.}$ ${\displaystyle x\in \Theta }$ ${\displaystyle y\in \Theta }$

${\displaystyle p_{xy}={\begin{cases}{\frac {1}{d}}{\frac {g(y)}{\sum _{z\in \Theta :z\sim _{j}x}g(z)}}&x\sim _{j}y\\0&{\text{otherwise}}\end{cases}}}$

Entonces

${\displaystyle g(x)p_{xy}={\frac {1}{d}}{\frac {g(x)g(y)}{\sum _{z\in \Theta :z\sim _{j}x}g(z)}}={\frac {1}{d}}{\frac {g(y)g(x)}{\sum _{z\in \Theta :z\sim _{j}y}g(z)}}=g(y)p_{yx}}$

ya que es una relación de equivalencia . Por lo tanto, se satisfacen las ecuaciones de equilibrio detalladas , lo que implica que la cadena es reversible y tiene una distribución invariante . ${\displaystyle x\sim _{j}y}$ ${\displaystyle \left.g\right.}$

En la práctica, el índice no se elige al azar, y la cadena recorre los índices en orden. En general, esto da como resultado un proceso de Markov no estacionario, pero cada paso individual seguirá siendo reversible y el proceso general seguirá teniendo la distribución estacionaria deseada (siempre que la cadena pueda acceder a todos los estados bajo el orden fijo). ${\displaystyle \left.j\right.}$

El sampler de Gibbs en la inferencia bayesiana y su relación con la teoría de la información

Sea , denotamos las observaciones generadas a partir de la distribución de muestreo y , una distribución a priori apoyada en el espacio de parámetros . Entonces, uno de los objetivos centrales de las estadísticas bayesianas es aproximar la densidad posterior. ${\displaystyle y}$ ${\displaystyle f(y|\theta )}$ ${\displaystyle \pi (\theta )}$ ${\displaystyle \Theta }$

${\displaystyle \pi (\theta |y)={\frac {f(y|\theta )\cdot \pi (\theta )}{m(y)}}}$

donde se supone que la probabilidad marginal es finita para todos . ${\displaystyle m(y)=\int _{\Theta }f(y|\theta )\cdot \pi (\theta )d\theta }$ ${\displaystyle y}$

Para explicar el muestreador de Gibbs, suponemos además que el espacio de parámetros se descompone como ${\displaystyle \Theta }$

${\displaystyle \Theta =\prod _{i=1}^{K}\Theta _{i}=\Theta _{1}\times \cdots \Theta _{i}\times \cdots \times \Theta _{K},\quad \quad (K>1)}$,

donde representa el producto cartesiano . Cada espacio de parámetros de componentes puede ser un conjunto de componentes escalares, subvectores o matrices. ${\displaystyle \times }$ ${\displaystyle \Theta _{i}}$

Defina un conjunto que complemente al . Los ingredientes esenciales del muestreador de Gibbs son la -ésima distribución posterior condicional completa para cada ${\displaystyle \Theta _{-i}}$ ${\displaystyle \Theta _{i}}$ ${\displaystyle i}$ ${\displaystyle i=1,\cdots ,K}$

${\displaystyle \pi (\theta _{i}|\theta _{-i},y)=\pi (\theta _{i}|\theta _{1},\cdots ,\theta _{i-1},\theta _{i+1},\cdots ,\theta _{K},y)}$.

Una descripción gráfica del algoritmo de muestreo de Gibbs ^[4]

Descripción esquemática de la igualdad de información asociada con el muestreador de Gibbs en el paso i-ésimo dentro de un ciclo ^[4]

El siguiente algoritmo detalla un muestreador de Gibbs genérico:

${\displaystyle {\text{Initialize: pick arbitrary starting value}}\,\,\theta ^{(1)}=(\theta _{1}^{(1)},\theta _{2}^{(1)},\cdots ,\theta _{i}^{(1)},\theta _{i+1}^{(1)},\cdots ,\theta _{K}^{(1)})}$

${\displaystyle {\text{Iterate a Cycle:}}\,}$

${\displaystyle \quad \quad {\text{Step 1. draw}}\,\,\theta _{1}^{(s+1)}\sim \pi (\theta _{1}|\theta _{2}^{(s)},\theta _{3}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad {\text{Step 2. draw}}\,\,\theta _{2}^{(s+1)}\sim \pi (\theta _{2}|\theta _{1}^{(s+1)},\theta _{3}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad \quad \vdots }$

${\displaystyle \quad \quad {\text{Step i. draw}}\,\,\theta _{i}^{(s+1)}\sim \pi (\theta _{i}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{i-1}^{(s+1)},\theta _{i+1}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad {\text{Step i+1. draw}}\,\,\theta _{i+1}^{(s+1)}\sim \pi (\theta _{i+1}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{i}^{(s+1)},\theta _{i+2}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad \quad \vdots }$

${\displaystyle \quad \quad {\text{Step K. draw}}\,\,\theta _{K}^{(s+1)}\sim \pi (\theta _{K}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{K-1}^{(s+1)},y)}$

${\displaystyle {\text{end Iterate}}}$

Tenga en cuenta que el muestreador de Gibbs se opera mediante el esquema iterativo de Monte Carlo dentro de un ciclo. La cantidad de muestras extraídas por el algoritmo anterior formula cadenas de Markov con la distribución invariante como la densidad objetivo . ${\displaystyle S}$ ${\displaystyle \{\theta ^{(s)}\}_{s=1}^{S}}$ ${\displaystyle \pi (\theta |y)}$

Ahora, para cada , defina las siguientes cantidades de teoría de la información: ${\displaystyle i=1,\cdots ,K}$

${\displaystyle I(\theta _{i};\theta _{-i})={\text{KL}}(\pi (\theta |y)||\pi (\theta _{i}|y)\cdot \pi (\theta _{-i}|y))=\int _{\Theta }\pi (\theta |y)\log {\bigg (}{\frac {\pi (\theta |y)}{\pi (\theta _{i}|y)\cdot \pi (\theta _{-i}|y)}}{\bigg )}d\theta ,}$

${\displaystyle H(\theta _{-i})=-\int _{\Theta _{-i}}\pi (\theta _{-i}|y)\log \pi (\theta _{-i}|y)d\theta _{-i},}$

${\displaystyle H(\theta _{-i}|\theta _{i})=-\int _{\Theta }\pi (\theta |y)\log \pi (\theta _{-i}|\theta _{i},y)d\theta ,}$

es decir, información mutua posterior, entropía diferencial posterior y entropía diferencial condicional posterior, respectivamente. De manera similar, podemos definir cantidades teóricas de información , y intercambiando y en las cantidades definidas. Entonces, se cumplen las siguientes ecuaciones. ^[4] ${\displaystyle I(\theta _{-i};\theta _{i})}$ ${\displaystyle H(\theta _{i})}$ ${\displaystyle H(\theta _{i}|\theta _{-i})}$ ${\displaystyle i}$ ${\displaystyle -i}$ ${\displaystyle K}$

${\displaystyle I(\theta _{i};\theta _{-i})=H(\theta _{-i})-H(\theta _{-i}|\theta _{i})=H(\theta _{i})-H(\theta _{i}|\theta _{-i})=I(\theta _{-i};\theta _{i}),\quad (i=1,\cdots ,K)}$.

La información mutua cuantifica la reducción de la incertidumbre de la cantidad aleatoria una vez que conocemos , a posteriori. Se desvanece si y solo si y son marginalmente independientes, a posteriori. La información mutua se puede interpretar como la cantidad que se transmite del paso -ésimo al paso -ésimo dentro de un solo ciclo del muestreador de Gibbs. ${\displaystyle I(\theta _{i};\theta _{-i})}$ ${\displaystyle \theta _{i}}$ ${\displaystyle \theta _{-i}}$ ${\displaystyle \theta _{i}}$ ${\displaystyle \theta _{-i}}$ ${\displaystyle I(\theta _{i};\theta _{-i})}$ ${\displaystyle i}$ ${\displaystyle i+1}$

Variaciones y ampliaciones

Existen numerosas variaciones del muestreador básico de Gibbs. El objetivo de estas variaciones es reducir la autocorrelación entre muestras lo suficiente como para superar cualquier costo computacional adicional.

Muestreador de Gibbs bloqueado

• Un muestreador de Gibbs bloqueado agrupa dos o más variables y toma muestras de su distribución conjunta condicionada a todas las demás variables, en lugar de tomar muestras de cada una de ellas individualmente. Por ejemplo, en un modelo oculto de Markov , un muestreador de Gibbs bloqueado podría tomar muestras de todas las variables latentes que componen la cadena de Markov de una sola vez, utilizando el algoritmo de avance-retroceso .

Muestreador de Gibbs colapsado

• Un muestreador de Gibbs colapsado integra ( margina ) una o más variables cuando muestrea para alguna otra variable. Por ejemplo, imagine que un modelo consta de tres variables A , B y C. Un muestreador de Gibbs simple muestrearía de p ( A  |  B , C ), luego p ( B  |  A , C ), luego p ( C  |  A , B ). Un muestreador de Gibbs colapsado podría reemplazar el paso de muestreo para A con una muestra tomada de la distribución marginal p ( A  |  C ), con la variable B integrada en este caso. Alternativamente, la variable B podría colapsarse por completo, muestreando alternativamente de p ( A  |  C ) y p ( C  |  A ) y sin muestrear sobre B en absoluto. La distribución sobre una variable A que surge al colapsar una variable padre B se llama distribución compuesta ; el muestreo de esta distribución es generalmente manejable cuando B es el conjugado anterior para A , particularmente cuando A y B son miembros de la familia exponencial . Para obtener más información, consulte el artículo sobre distribuciones compuestas o Liu (1994). ^[5]

Implementación de un sampler de Gibbs colapsado

Distribuciones de Dirichlet colapsantes

En los modelos bayesianos jerárquicos con variables categóricas , como la asignación latente de Dirichlet y otros modelos diversos utilizados en el procesamiento del lenguaje natural , es bastante común colapsar las distribuciones de Dirichlet que se utilizan normalmente como distribuciones previas sobre las variables categóricas. El resultado de este colapso introduce dependencias entre todas las variables categóricas que dependen de una distribución previa de Dirichlet dada, y la distribución conjunta de estas variables después del colapso es una distribución multinomial de Dirichlet . La distribución condicional de una variable categórica dada en esta distribución, condicionada a las demás, asume una forma extremadamente simple que hace que el muestreo de Gibbs sea aún más fácil que si no se hubiera realizado el colapso. Las reglas son las siguientes:

1. El colapso de un nodo anterior de Dirichlet afecta solo a los nodos padre e hijos del anterior. Dado que el padre suele ser una constante, normalmente solo debemos preocuparnos por los hijos.
2. Al colapsar un prior de Dirichlet se introducen dependencias entre todos los hijos categóricos que dependen de ese prior, pero no se introducen dependencias adicionales entre ningún otro hijo categórico. (Esto es importante tenerlo en cuenta, por ejemplo, cuando hay varios priores de Dirichlet relacionados por el mismo hiperprior. Cada prior de Dirichlet se puede colapsar de forma independiente y afecta solo a sus hijos directos).
3. Después de colapsar, la distribución condicional de un hijo dependiente sobre los otros asume una forma muy simple: la probabilidad de ver un valor dado es proporcional a la suma del hiperprior correspondiente para este valor y el recuento de todos los demás nodos dependientes que asumen el mismo valor. Los nodos que no dependen del mismo prior no deben contarse. La misma regla se aplica en otros métodos de inferencia iterativa, como Bayes variacional o maximización de expectativas ; sin embargo, si el método implica mantener recuentos parciales, entonces los recuentos parciales para el valor en cuestión deben sumarse en todos los demás nodos dependientes. A veces, este recuento parcial sumado se denomina recuento esperado o similar. La probabilidad es proporcional al valor resultante; la probabilidad real debe determinarse normalizando todos los valores posibles que puede tomar la variable categórica (es decir, sumando el resultado calculado para cada valor posible de la variable categórica y dividiendo todos los resultados calculados por esta suma).
4. Si un nodo categórico dado tiene hijos dependientes (por ejemplo, cuando es una variable latente en un modelo mixto ), el valor calculado en el paso anterior (recuento esperado más el anterior, o lo que se calcule) debe multiplicarse por las probabilidades condicionales reales ( ¡no un valor calculado que sea proporcional a la probabilidad!) de todos los hijos dados sus padres. Consulte el artículo sobre la distribución multinomial de Dirichlet para obtener una explicación detallada.
5. En el caso en que la pertenencia a un grupo de nodos que dependen de una distribución Dirichlet previa dada pueda cambiar dinámicamente dependiendo de alguna otra variable (por ejemplo, una variable categórica indexada por otra variable categórica latente, como en un modelo de temas ), se siguen calculando los mismos recuentos esperados, pero deben hacerse con cuidado para que se incluya el conjunto correcto de variables. Consulte el artículo sobre la distribución multinomial de Dirichlet para obtener más información, incluso en el contexto de un modelo de temas.

Colapso de otros priores conjugados

En general, cualquier conjugado anterior puede colapsarse, si sus únicos hijos tienen distribuciones conjugadas a él. Las matemáticas relevantes se discuten en el artículo sobre distribuciones compuestas . Si solo hay un nodo hijo, el resultado a menudo asumirá una distribución conocida. Por ejemplo, colapsar una varianza distribuida en gamma inversa de una red con un solo hijo gaussiano producirá una distribución t de Student . (De hecho, colapsar tanto la media como la varianza de un solo hijo gaussiano producirá una distribución t de Student, siempre que ambas sean conjugadas, es decir, media gaussiana, varianza en gamma inversa).

Si hay varios nodos secundarios, todos se volverán dependientes, como en el caso categórico de Dirichlet . La distribución conjunta resultante tendrá una forma cerrada que se asemeja en algunos aspectos a la distribución compuesta, aunque tendrá un producto de varios factores, uno para cada nodo secundario.

Además, y lo más importante, la distribución condicional resultante de uno de los nodos secundarios dados los demás (y también dados los padres de los nodos colapsados, pero no dados los hijos de los nodos secundarios) tendrá la misma densidad que la distribución predictiva posterior de todos los nodos secundarios restantes. Además, la distribución predictiva posterior tiene la misma densidad que la distribución compuesta básica de un solo nodo, aunque con diferentes parámetros. La fórmula general se proporciona en el artículo sobre distribuciones compuestas .

Por ejemplo, dada una red de Bayes con un conjunto de nodos distribuidos de manera idéntica y condicionalmente independientes con distribución gaussiana y con distribuciones previas conjugadas en la media y la varianza, la distribución condicional de un nodo dados los otros después de calcular la media y la varianza será una distribución t de Student . De manera similar, el resultado de calcular la distribución previa gamma de varios nodos distribuidos de Poisson hace que la distribución condicional de un nodo dados los otros asuma una distribución binomial negativa .

En estos casos en los que la composición produce una distribución bien conocida, a menudo existen procedimientos de muestreo eficientes, y su uso será a menudo (aunque no necesariamente) más eficiente que no colapsar y, en su lugar, muestrear por separado los nodos anteriores y secundarios. Sin embargo, en el caso en que la distribución compuesta no sea bien conocida, puede que no sea fácil muestrear, ya que generalmente no pertenecerá a la familia exponencial y, por lo general, no será logarítmicamente cóncava (lo que facilitaría el muestreo mediante muestreo de rechazo adaptativo , ya que siempre existe una forma cerrada).

En el caso en que los nodos hijos de los nodos colapsados ​​tengan hijos, la distribución condicional de uno de estos nodos hijos dados todos los demás nodos en el gráfico tendrá que tener en cuenta la distribución de estos hijos de segundo nivel. En particular, la distribución condicional resultante será proporcional a un producto de la distribución compuesta como se definió anteriormente, y las distribuciones condicionales de todos los nodos hijos dados sus padres (pero no dados sus propios hijos). Esto se desprende del hecho de que la distribución condicional completa es proporcional a la distribución conjunta. Si los nodos hijos de los nodos colapsados ​​son continuos , esta distribución generalmente no tendrá una forma conocida, y puede ser difícil de muestrear a pesar del hecho de que se puede escribir una forma cerrada, por las mismas razones que se describieron anteriormente para distribuciones compuestas no bien conocidas. Sin embargo, en el caso particular de que los nodos hijos sean discretos , el muestreo es factible, independientemente de si los hijos de estos nodos hijos son continuos o discretos. De hecho, el principio involucrado aquí se describe con bastante detalle en el artículo sobre la distribución multinomial de Dirichlet .

Muestreador de Gibbs con sobrerelajación ordenada

• Un muestreador de Gibbs con sobrerelajación ordenada muestrea un número impar dado de valores candidatos para en cualquier paso dado y los ordena, junto con el valor único para de acuerdo con un orden bien definido. Si es el s ^º más pequeño en la lista ordenada, entonces se selecciona como el s ^º más grande en la lista ordenada. Para obtener más información, consulte Neal (1995). ^[6] ${\displaystyle x_{j}^{(i)}}$ ${\displaystyle x_{j}^{(i-1)}}$ ${\displaystyle x_{j}^{(i-1)}}$ ${\displaystyle x_{j}^{(i)}}$

Otras extensiones

También es posible ampliar el muestreo de Gibbs de varias maneras. Por ejemplo, en el caso de variables cuya distribución condicional no es fácil de muestrear, se puede utilizar una única iteración de muestreo por cortes o el algoritmo Metropolis-Hastings para muestrear las variables en cuestión. También es posible incorporar variables que no sean variables aleatorias , pero cuyo valor se calcule de manera determinista a partir de otras variables. Los modelos lineales generalizados , por ejemplo, la regresión logística (también conocidos como " modelos de máxima entropía "), se pueden incorporar de esta manera. (BUGS, por ejemplo, permite este tipo de mezcla de modelos).

Modos de falla

Hay dos formas en las que el muestreo de Gibbs puede fallar. La primera es cuando hay islas de estados de alta probabilidad, sin caminos entre ellos. Por ejemplo, considere una distribución de probabilidad sobre vectores de 2 bits, donde los vectores (0,0) y (1,1) tienen cada uno una probabilidad ⁠1/2⁠ , pero los otros dos vectores (0,1) y (1,0) tienen probabilidad cero. El muestreo de Gibbs quedará atrapado en uno de los dos vectores de alta probabilidad y nunca alcanzará al otro. De manera más general, para cualquier distribución sobre vectores de alta dimensión y valor real, si dos elementos particulares del vector están perfectamente correlacionados (o perfectamente anticorrelacionados), esos dos elementos quedarán atrapados y el muestreo de Gibbs nunca podrá cambiarlos.

El segundo problema puede ocurrir incluso cuando todos los estados tienen una probabilidad distinta de cero y solo hay una única isla de estados de alta probabilidad. Por ejemplo, considere una distribución de probabilidad sobre vectores de 100 bits, donde el vector de todos ceros ocurre con probabilidad ⁠1/2⁠ , y todos los demás vectores son igualmente probables, y por lo tanto tienen una probabilidad de cada uno. Si desea estimar la probabilidad del vector cero, sería suficiente tomar 100 o 1000 muestras de la distribución verdadera. Eso muy probablemente daría una respuesta muy cercana a ⁠ ${\displaystyle {\frac {1}{2(2^{100}-1)}}}$1/2⁠ . Pero probablemente tendrías que tomar más muestras que las del muestreo de Gibbs para obtener el mismo resultado. Ninguna computadora podría hacer esto en toda su vida. ${\displaystyle 2^{100}}$

Este problema ocurre sin importar cuán largo sea el período de rodaje. Esto se debe a que en la distribución verdadera, el vector cero ocurre la mitad del tiempo, y esas ocurrencias se mezclan aleatoriamente con los vectores distintos de cero. Incluso una muestra pequeña verá vectores tanto cero como distintos de cero. Pero el muestreo de Gibbs alternará entre devolver solo el vector cero durante largos períodos (aproximadamente seguidos) y luego solo vectores distintos de cero durante largos períodos (aproximadamente seguidos). Por lo tanto, la convergencia a la distribución verdadera es extremadamente lenta y requiere mucho más que pasos; dar tantos pasos no es computacionalmente factible en un período de tiempo razonable. La convergencia lenta aquí puede verse como una consecuencia de la maldición de la dimensionalidad . Un problema como este se puede resolver muestreando en bloque todo el vector de 100 bits a la vez. (Esto supone que el vector de 100 bits es parte de un conjunto más grande de variables. Si este vector es lo único que se muestrea, entonces el muestreo en bloque es equivalente a no hacer el muestreo de Gibbs en absoluto, lo que por hipótesis sería difícil). ${\displaystyle 2^{99}}$ ${\displaystyle 2^{99}}$ ${\displaystyle 2^{99}}$

Software

• El software OpenBUGS ( Bayesian inference Using Gibbs Sampling ) realiza un análisis bayesiano de modelos estadísticos complejos utilizando cadenas de Markov Monte Carlo .

• JAGS ( Just another Gibbs sampler ) es un programa GPL para el análisis de modelos jerárquicos bayesianos utilizando Markov Chain Monte Carlo.

• Church es un software libre para realizar inferencias de Gibbs sobre distribuciones arbitrarias que se especifican como programas probabilísticos.

• PyMC es una biblioteca Python de código abierto para el aprendizaje bayesiano de modelos gráficos probabilísticos generales .
• Turing es una biblioteca Julia de código abierto para inferencia bayesiana que utiliza programación probabilística .

Notas

1. Geman, S. ; Geman, D. (1984). "Relajación estocástica, distribuciones de Gibbs y restauración bayesiana de imágenes". IEEE Transactions on Pattern Analysis and Machine Intelligence . 6 (6): 721–741. doi :10.1109/TPAMI.1984.4767596. PMID  22499653.
2. Gelfand, Alan E.; Smith, Adrian FM (1 de junio de 1990). "Enfoques basados ​​en muestreo para calcular densidades marginales". Revista de la Asociación Estadounidense de Estadística . 85 (410): 398–409. doi :10.1080/01621459.1990.10476213. ISSN  0162-1459.
3. Gelman, Andrew y Carlin, John B y Stern, Hal S y Dunson, David B y Vehtari, Aki y Rubin, Donald B (2014). Análisis de datos bayesianos . Vol. 2. FL: CRC press Boca Raton.{{cite book}}: CS1 maint: multiple names: authors list (link)
4. Lee, Se Yoon (2021). "Inferencia variacional de ascenso de coordenadas y muestreador de Gibbs: una revisión de la teoría de conjuntos". Comunicaciones en estadística: teoría y métodos . 51 (6): 1549–1568. arXiv : 2008.01006 . doi :10.1080/03610926.2021.1921214. S2CID  220935477.
5. Liu, Jun S. (septiembre de 1994). "El muestreador de Gibbs colapsado en cálculos bayesianos con aplicaciones a un problema de regulación genética". Revista de la Asociación Estadounidense de Estadística . 89 (427): 958–966. doi :10.2307/2290921. JSTOR  2290921.
6. Neal, Radford M. (1995). Supresión de paseos aleatorios en el método Monte Carlo de cadenas de Markov mediante sobrerelajación ordenada (informe técnico). Universidad de Toronto, Departamento de Estadística. arXiv : bayes-an/9506004 . Código Bibliográfico :1995bayes.an..6004N.

Referencias

• Bishop, Christopher M. (2006), Reconocimiento de patrones y aprendizaje automático , Springer, ISBN 978-0-387-31073-2
• Bolstad, William M. (2010), Comprensión de las estadísticas bayesianas computacionales , John Wiley ISBN 978-0-470-04609-8 
• Casella, G.; George, EI (1992). "Explicación del muestreador de Gibbs". The American Statistician . 46 (3): 167. CiteSeerX  10.1.1.554.3993 . doi :10.2307/2685208. JSTOR  2685208.(Contiene un resumen básico y muchas referencias).
• Gelfand, Alan E.; Smith, Adrian FM (1990), "Enfoques basados ​​en muestreo para calcular densidades marginales", Journal of the American Statistical Association , 85 (410): 398–409, doi :10.2307/2289776, JSTOR  2289776, MR  1141740
• Gelman, A. , Carlin JB, Stern HS, Dunson D., Vehtari A., Rubin DB (2013), Análisis de datos bayesianos , tercera edición. Londres: Chapman & Hall .
• Levin, David A.; Peres, Yuval ; Wilmer, Elizabeth L. (2008), " Cadenas de Markov y tiempos de mezcla ", American Mathematical Society .
• Robert, CP; Casella, G. (2004), Métodos estadísticos de Monte Carlo (segunda edición), Springer-Verlag.