muestreo de gibbs

Algoritmo de Montecarlo

En estadística , el muestreo de Gibbs o un muestreador de Gibbs es un algoritmo Monte Carlo de cadena de Markov (MCMC) para muestrear a partir de una distribución de probabilidad multivariada específica cuando el muestreo directo de la distribución conjunta es difícil, pero el muestreo 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); aproximar la distribución marginal de una de las variables, o algún subconjunto de las variables (por ejemplo, los parámetros desconocidos o las 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 tanto, no es necesario muestrearlas.

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

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

Introducción

El muestreo de Gibbs lleva el 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 ampliadas (ver más abajo), se puede considerar un marco general para el muestreo de un gran conjunto de variables, muestreando cada variable (o en algunos casos, cada grupo de variables) por turno, y puede incorporar Metropolis– Algoritmo de Hastings (o métodos como el muestreo por cortes ) 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 tomar muestras directamente, pero la distribución condicional de cada variable se conoce y es fácil (o al menos más fácil) tomar muestras. El algoritmo de muestreo de Gibbs genera una instancia a partir de la distribución de cada variable, condicionada 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 generalmente se especifican como una colección de distribuciones condicionales.

Implementación

El muestreo de Gibbs, en su encarnación básica, es un caso especial del algoritmo Metropolis-Hastings . El objetivo del muestreo de Gibbs es que, dada una distribución multivariada , es más sencillo muestrear a partir de una distribución condicional que marginar integrando sobre una distribución conjunta . Supongamos que queremos obtener muestras de una distribución conjunta . Denota la enésima muestra 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. Llame a esta siguiente muestra . Como es un vector, tomamos muestras de cada componente del vector, , a partir de la distribución de ese componente condicionada a todos los demás componentes muestreados hasta el momento. Pero hay un problema: condicionamos los componentes de ' hasta , y luego condicionamos los componentes de ' , comenzando desde hasta . Para lograr esto, tomamos muestras de 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 enésimo en la muestra enésima, no en la muestra enésima. ${\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 veces. ${\displaystyle k}$

Propiedades

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

• Las muestras se aproximan a la distribución conjunta de todas las variables.
• La distribución marginal de cualquier subconjunto de variables se puede aproximar simplemente considerando las muestras de 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 se pueden determinar aleatoriamente o mediante algún otro algoritmo como el de maximización de expectativas .
• 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 preparación ) y luego considerar solo una de cada cuatro muestras al promediar los valores para calcular una expectativa. Por ejemplo, se pueden ignorar las primeras 1000 muestras y luego promediar cada 100 muestras, desechando el resto. La razón de 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 algún tiempo alcanzar esa distribución estacionaria; (2) las muestras sucesivas no son independientes entre sí, sino que forman una cadena de Markov con cierta correlación. A veces, se pueden usar algoritmos para determinar la cantidad de autocorrelación entre muestras y el valor (el período entre muestras que realmente se usan) 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 " paseo aleatorio " 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 ; vea abajo.

Relación de 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 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 de acuerdo con 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 restablecen 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 puede ignorarse, ya que la mayoría de los métodos de muestreo no la requieren.

Inferencia

El muestreo de Gibbs se utiliza comúnmente para inferencia estadística (por ejemplo, determinar el mejor valor de un parámetro, como determinar el número de personas que probablemente comprarán en una tienda en particular en un día determinado, el candidato por el que probablemente votará un votante, 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 tomar muestras 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 con mayor frecuencia; 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 "agrupaciones" para obtener una estimación significativa de la moda). Sin embargo, lo más común es que el valor esperado se elige el valor ( media o promedio) de los valores muestreados; Este es un estimador 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) es capaz de devolver solo un 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 efectivamente representa el extra. masa de probabilidad en esa dirección. (Si una distribución es multimodal, es posible que el valor esperado no devuelva un punto significativo y cualquiera de las modas suele ser una mejor opción).

Aunque algunas de las variables normalmente corresponden a parámetros de interés, otras son variables poco interesantes ("molestas") introducidas en el modelo para expresar adecuadamente las relaciones entre variables. Aunque los valores muestreados representan la distribución conjunta de todas las variables, las variables molestas pueden simplemente ignorarse al calcular los valores o modas esperados; esto equivale a marginar a las variables molestas. Cuando se desea un valor para múltiples variables, el valor esperado simplemente se calcula para 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.

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

En ocasiones, el muestreo de Gibbs también puede manejar modelos lineales generalizados (es decir, variaciones de la regresión lineal ). Por ejemplo, la regresión probit para determinar la probabilidad de una elección binaria determinada (sí/no), con antecedentes normalmente distribuidos 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 (normalmente 7 a 9) de distribuciones normales. Sin embargo, lo más habitual es que se utilice 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 probar de la siguiente manera. Defina if para todos y denote 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 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.}$

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

Denotemos observaciones generadas a partir de la distribución muestral y sea un soporte previo en el espacio de parámetros . Entonces uno de los objetivos centrales de la estadística bayesiana 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, asumimos 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 el . El ingrediente esencial del muestreador de Gibbs es la ené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 i-ésimo paso dentro de un ciclo ^[4]

El siguiente algoritmo detalla una muestra genérica de Gibbs:

${\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 funciona mediante el esquema iterativo de Monte Carlo dentro de un ciclo. El número de muestras extraídas por el algoritmo anterior formula cadenas de Markov con la distribución invariante como densidad objetivo . ${\displaystyle S}$ ${\displaystyle \{\theta ^{(s)}\}_{s=1}^{S}}$ ${\displaystyle \pi (\theta |y)}$

Ahora, para cada , defina las siguientes cantidades teóricas de 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 ,}$

a saber, 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 una cantidad aleatoria una vez que la conocemos , a posteriori. Desaparece si y sólo si y son marginalmente independientes, a posteriori. La información mutua se puede interpretar como la cantidad que se transmite desde el -ésimo paso al -ésimo paso 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 extensiones

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.

Muestra de Gibbs bloqueada

• Un muestreador de Gibbs bloqueado agrupa dos o más variables y toma muestras de su distribución conjunta condicionadas a todas las demás variables, en lugar de tomar muestras de cada una individualmente. Por ejemplo, en un modelo de Markov oculto , 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 .

Muestra de Gibbs colapsada

• Un muestreador de Gibbs colapsado integra ( margina ) una o más variables cuando toma muestras de alguna otra variable. Por ejemplo, imagine que un modelo consta de tres variables A , B y C. Un muestreador de Gibbs simple tomaría muestras 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 de 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 p ( A  |  C ) y p ( C  |  A ) y no muestreando B en absoluto. La distribución sobre una variable A que surge al colapsar una variable principal B se llama distribución compuesta ; el muestreo de esta distribución es generalmente manejable cuando B es el conjugado anterior de 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 de Liu (1994). ^[5]

Implementación de un muestreador de Gibbs colapsado

Colapso de las distribuciones de Dirichlet

En los modelos bayesianos jerárquicos con variables categóricas , como la asignación latente de Dirichlet y varios otros modelos utilizados en el procesamiento del lenguaje natural , es bastante común colapsar las distribuciones de Dirichlet que normalmente se utilizan como distribuciones previas sobre las variables categóricas. El resultado de este colapso introduce dependencias entre todas las variables categóricas que dependen de un Dirichlet anterior determinado, 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 colapsado. Las reglas son las siguientes:

1. Colapsar un nodo anterior de Dirichlet afecta solo a los nodos padre e hijo del anterior. Dado que los padres suelen ser una constante, normalmente sólo debemos preocuparnos de los niños.
2. Colapsar un anterior de Dirichlet introduce dependencias entre todos los hijos categóricos que dependen de ese anterior, pero no dependencias adicionales entre ningún otro hijo categórico. (Es importante tener esto en cuenta, por ejemplo, cuando hay varios anteriores de Dirichlet relacionados por el mismo hiperprior. Cada anterior de Dirichlet se puede contraer de forma independiente y afecta solo a sus hijos directos).
3. Después del colapso, la distribución condicional de uno de los hijos dependientes sobre los demás asume una forma muy simple: la probabilidad de ver un valor dado es proporcional a la suma del hiperprior correspondiente a este valor, y el recuento de todos los demás nodos dependientes supone el mismo valor. Los nodos que no dependen del mismo anterior no deben contarse. La misma regla se aplica en otros métodos de inferencia iterativos, como el Bayes variacional o la maximización de expectativas ; sin embargo, si el método implica mantener recuentos parciales, entonces los recuentos parciales del valor en cuestión se deben sumar en todos los demás nodos dependientes. A veces, este recuento parcial resumido 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 es proporcional a la probabilidad!) de todos los niños dados sus padres. Consulte el artículo sobre la distribución multinomial de Dirichlet para obtener una discusión detallada.
5. En el caso en que la pertenencia al grupo de los nodos que dependen de un Dirichlet anterior determinado 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 tema ), todavía se calculan los mismos recuentos esperados. , pero debe hacerse con cuidado para incluir 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 temático.

Colapsando otros anteriores conjugados

En general, cualquier anterior conjugado se puede contraer, si sus únicos hijos tienen distribuciones conjugadas. Las matemáticas relevantes se analizan 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 aún producirá una distribución t de Student, siempre que ambas sean conjugadas, es decir, media gaussiana, varianza 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 parecerá en algunos aspectos a la distribución compuesta, aunque tendrá un producto de varios factores, uno para cada nodo hijo.

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 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 da en el artículo sobre distribuciones compuestas .

Por ejemplo, dada una red Bayes con un conjunto de nodos con distribución gaussiana condicionalmente independientes e idénticamente distribuidos con distribuciones previas conjugadas colocadas en la media y la varianza, la distribución condicional de un nodo dados los demás después de componer tanto la media como la varianza será una Distribución t de Student . De manera similar, el resultado de combinar la gama anterior de varios nodos con distribución de Poisson hace que la distribución condicional de un nodo dados los demás asuma una distribución binomial negativa .

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

En el caso de que los nodos secundarios de los nodos colapsados ​​tengan hijos, la distribución condicional de uno de estos nodos secundarios dados todos los demás nodos en el gráfico deberá 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 define anteriormente y las distribuciones condicionales de todos los nodos secundarios dados sus padres (pero no 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 secundarios de los nodos colapsados ​​son continuos , esta distribución generalmente no tendrá una forma conocida y puede ser difícil tomar muestras a pesar de que se puede escribir una forma cerrada, por las mismas razones descritas anteriormente para los nodos no colapsados. -Distribuciones compuestas conocidas. Sin embargo, en el caso particular de que los nodos secundarios sean discretos , el muestreo es factible, independientemente de si los hijos de estos nodos secundarios 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 .

Muestra de Gibbs con sobrerelajación ordenada

• Un muestreador de Gibbs con sobrerelajación ordenada toma muestras de un número impar dado de valores candidatos para en cualquier paso dado y los clasifica, junto con el valor único de acuerdo con un orden bien definido. Si es el enésimo ^más pequeño en la lista ordenada, entonces se selecciona como el enésimo ^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 tomar como muestra, se puede utilizar una única iteración de muestreo por cortes o el algoritmo de Metropolis-Hastings para tomar muestras de las variables en cuestión. También es posible incorporar variables que no sean variables aleatorias , pero cuyo valor se calcule de forma determinista a partir de otras variables. De esta manera se pueden incorporar modelos lineales generalizados , por ejemplo, regresión logística (también conocidos como " modelos de máxima entropía "). (BUGS, por ejemplo, permite este tipo de mezcla de modelos).

Modos de fallo

Hay dos formas en 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 ½, 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 llegará 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 se atascarán y el muestreo de Gibbs nunca podrá cambiar. a ellos.

El segundo problema puede ocurrir incluso cuando todos los estados tienen una probabilidad distinta de cero y solo hay una 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 ½ 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. Es muy probable que esto dé una respuesta muy cercana a ½. Pero probablemente habría que tomar más que muestras del muestreo de Gibbs para obtener el mismo resultado. Ninguna computadora podría hacer esto en toda su vida. ${\displaystyle {\frac {1}{2(2^{100}-1)}}}$ ${\displaystyle 2^{100}}$

Este problema ocurre sin importar cuánto tiempo dure el período de precalentamiento. 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 cero y distintos de cero. Pero el muestreo de Gibbs alternará entre devolver solo el vector cero durante períodos largos (aproximadamente seguidos) y luego solo vectores distintos de cero durante períodos largos (aproximadamente seguidos). Por tanto, la convergencia hacia la distribución verdadera es extremadamente lenta y requiere mucho más que pasos; tomar tantos pasos no es computacionalmente factible en un período de tiempo razonable. La lenta convergencia 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 está muestreando, entonces el muestreo en bloque equivale a no realizar ningún muestreo de Gibbs, lo que, según la hipótesis, sería difícil). ${\displaystyle 2^{99}}$ ${\displaystyle 2^{99}}$ ${\displaystyle 2^{99}}$

Software

• El software OpenBUGS ( Inferencia bayesiana utilizando Gibbs Sampling ) realiza un análisis bayesiano de modelos estadísticos complejos utilizando la cadena 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 gratuito para realizar inferencias de Gibbs sobre distribuciones arbitrarias especificadas 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. Alemán, S .; Alemán, D. (1984). "Relajación estocástica, distribuciones de Gibbs y restauración bayesiana de imágenes". Transacciones IEEE sobre análisis de patrones e inteligencia artificial . 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 presiona Boca Ratón.{{cite book}}: CS1 maint: multiple names: authors list (link)
4. Lee, Se Yoon (2021). "Muestreador de Gibbs e inferencia variacional de ascenso de coordenadas: 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 caminatas aleatorias en la cadena de Markov Monte Carlo 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 la estadística bayesiana computacional , John Wiley ISBN 978-0-470-04609-8 
• Casella, G.; George, EI (1992). "Explicando el muestreador de Gibbs". El estadístico estadounidense . 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", Revista de la Asociación Estadounidense de Estadística , 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 ", Sociedad Matemática Estadounidense .
• Roberto, CP; Casella, G. (2004), Métodos estadísticos de Monte Carlo (segunda edición), Springer-Verlag.