Convergencia MCMC

Question

Intento ajustar un modelo bayesiano. Este modelo tiene muchos parámetros (unos 150).

Ejecuto MCMC (10000 iters) con un periodo de adelgazamiento para evitar el problema de la correlación y con la esperanza de que me ahorre algo de tiempo de cálculo. Así obtengo la muestra que me permite calcular las estimaciones de Monte Carlo. Sabía que una inferencia adecuada implica comprobar la convergencia MCMC. Llevo a cabo esta tarea con una evaluación gráfica (sé que hay enfoques más formales, este es un primer intento) con trace plot. Me di cuenta de que hay pocos parámetros que no muestran ninguna convergencia. ¿Qué significa esto?

Cuando intento ver cuál era el tamaño efectivo de la muestra, obtengo para estos parámetros valores muy bajos (como 6) y esto es señal de muestras correlacionadas.

No sé cómo seguir.

ACTUALIZACIÓN: Este es mi modelo en BUGS. Es una regresión dirichlet, con la variable de respuesta y compuesta por 4 categorías (por lo que K=4). Lo estoy aplicando a un conjunto de datos relacionados con los delitos violentos en el estado de EE.UU.. Así que mi idea era utilizar un parámetro beta para cada observación i (los estados en EE.UU.) Ejemplo de datos:

     state  violent_crime    murder  rape robbery aggravated_assault   median_income
1   Alabama       20727    276  2005    4701              13745     42917 
2   alaska        4684     41   771     629               3243      70898

Modelo:

model
{
    for (i in 1:n) {
        beta1[i,1:K]~dmnorm(zero[],t[1:K,1:K])

        Y[i,1:K] ~ddirich(alpha[i,1:K])
        for ( k in 1:4)
        {
        xb[i,k]<-exp(beta0[k]+beta1[i,k]*x[i])
        }
        for( k in 1:3)
        {
        mu[i,k]<-xb[i,k]/sum(xb[i,1:4])
        }
        mu[i,K]<-1/sum(xb[i,1:4])
        for( k in 1:K)
        {
        alpha[i,k]<-mu[i,k]*phi
        }

        }
for(i in 1:K)
     { 
 for (j in 1:K) 
       {  B[i,j] <- 0.00001*equals(i,j)
      }}
#prior
beta0[1:K]~dmnorm(zero[],B[,])

    for( k in 1:K)
    {
     zero[k] <- 0
    }

#hyper-prior
t[1:K,1:K]~dwish(R[1:K,1:K],4)
phi~dgamma(0.0001,0.0001)

}

Answer 1

Significa que lo más probable es que su cadena no haya convergido. Con esto quiero decir que debe tener cuidado con el todo cadena, no sólo preocuparse por las dimensiones con bajo número efectivo de muestras.

Soluciones

Burn-in (descartar la parte inicial de la cadena) - véase también esta pregunta

A veces, un número efectivo bajo de muestras se debe simplemente a que la cadena comenzó en una región de baja probabilidad y no encontró la cuenca de convergencia (la región de alta probabilidad o el conjunto típico) hasta más tarde. Yo hago no recomiendan que ahora solo juegues con el burn-in, ya que con una sola cadena es difícil saber lo que pasa.

En su lugar, es posible que desee ejecutar varias cadenas en paralelo, en lugar de una sola cadena larga (véase el libro de MacKay, Capítulo 29 ). Con múltiples cadenas suele ser más fácil detectar la falta de convergencia, aunque aquí hay diferentes escuelas de pensamiento sobre el número de cadenas (yo suelo hacer 3 o 4).

En cuanto al rodaje, también hay varias opiniones. Algunas personas, como Geyer dicen que no tiene sentido, siempre y cuando estés seguro de que empiezas en una región de alta probabilidad:

Cualquier punto que no te importe tener en una muestra es un buen punto de partida.

Sin embargo, esto es más fácil de decir que de hacer en 150 dimensiones. Los paquetes estadísticos más avanzados llegan a quemar hasta el 50% de la cadena (véase Stan ), pero una de las razones de un rodaje tan largo es que también se utiliza para adaptar algunos de los parámetros del muestreador.

En tu caso, te recomiendo que inicialices tu muestreador a partir de una buena estimación (por ejemplo, en algún lugar alrededor de la moda es mejor que nada, aunque la moda en sí probablemente no esté en el conjunto típico), y que hagas algún burn-in ya que en altas dimensiones es difícil saber dónde está la probabilidad masa reside (no es lo mismo que la densidad de probabilidad). En mi opinión, es mejor quemar más que menos (en el peor de los casos, simplemente se están desechando algunas muestras efectivas, mientras que si no se quema es posible que se conserven muestras que no son representativas de la distribución objetivo).

Ver más

Es bastante sencillo: comience donde terminó su cadena MCMC (o a partir de alguna otra buena suposición, como se ha mencionado anteriormente), descarte la cadena anterior y tome más muestras en total.

Afina mejor tu sampler

No sé qué método MCMC está utilizando, pero la mayoría de los muestreadores tienen varios hiperparámetros sintonizables (por ejemplo, longitud(es) de salto para Metropolis-Hastings con propuesta gaussiana, matriz de masa para Hamiltonian Monte Carlo, etc.). El rendimiento de la mayoría de los muestreadores MCMC depende en gran medida de la elección de estos parámetros, por lo que podría echar un vistazo a su distribución muestreada y averiguar si puede mejorarlos.

Cambiar el muestreador

A veces, el muestreador que está utilizando puede no ser la mejor opción. Por ejemplo, si su problema tiene que ver con variables continuas con una distribución objetivo suave y diferenciable, probablemente querrá utilizar alguna variante de Hamiltonian Monte Carlo (el estado del arte es NUTS tal como se aplica en Stan ).

Cambiar la parametrización

Como sugiere @Björn en los comentarios, quizá quieras probar una parametrización diferente de tu modelo. Las reparametrizaciones sencillas consisten en una transformación lineal o no lineal de las variables; la primera, por ejemplo, para reducir la correlación entre variables, y la segunda, por ejemplo, para eliminar las colas largas y gruesas de la curva posterior, de las que puede ser difícil extraer muestras. Las reparametrizaciones más complejas incluyen la adición de variables auxiliares que podrían ayudar a la mezcla (por ejemplo, tendiendo puentes entre diferentes regiones de la variable posterior), pero esto depende en gran medida del modelo.