Mezcla de distribuciones Normal y Exponencial con pesos y parámetros desconocidos utilizando JAGS

Question

Estoy tratando de modelar un proceso que genera variables aleatorias normales y variables aleatorias exponenciales, con pesos desconocidos entre las dos distribuciones, y parámetros desconocidos de las distribuciones. Estoy tratando de hacer esto usando rjags en R.

He encontrado este puesto contestado por Matt Denwood que mezcla distribuciones normales y uniformes, pero por alguna razón no fui capaz de adaptarla a mi problema (hay que reconocer que sólo cambié la uniforme por la exponencial, sin ser capaz de entender el "truco de los" mencionado en el post).

Distribución de mis datos:

Salida de JAGS:

Mi modelo:

model{

for(i in 1:N){

    # Log density for the normal part:
    ld_norm[i] <- logdensity.norm(x[i], mu, tau)

    # Log density for the uniform part:
    ld_exp[i] <- logdensity.exp(x[i], 1/lambda)

    # Select one of these two densities:
    density[i] <- exp(ld_norm[i]*norm_chosen[i] + ld_exp[i]*(1-norm_chosen[i]))

    # Generate a likelihood for the MCMC sampler:
    Ones[i] ~ dbern(density[i])

    # The latent class part as usual:
    norm_chosen[i] ~ dbern(prob)

}

# Priors:
lambda ~ dnorm(1, 10^-6)
prob ~ dbeta(1,1)
mu ~ dnorm(0, 10^-6)
tau <- pow(sigma, -2)
sigma ~ dnorm(1, 10^-6)

}

¿Alguna idea de por qué este método no me funciona?

Answer 1

Sus posteriors se parecen sospechosamente a sus priors, lo que suele indicar que su modelo no se está ajustando a los datos. Mi mejor suposición es que no has incluido correctamente el vector "Ones" con tus datos, aunque esto es sólo una suposición, ya que tiene que ver con tu código R y no con el modelo JAGS en sí.

Pero, en cualquier caso, puedo intentar explicar el "truco de Ones", que puede ayudar a enfatizar por qué es necesario incluirlo como dato. Básicamente es una forma de "engañar" a JAGS/BUGS para que ajuste una función de probabilidad personalizada a tus datos. Implica el cálculo de la función de probabilidad logarítmica para cada punto de datos utilizando algún método arbitrario - en su caso la densidad[i] se calcula a partir de una distribución normal o exponencial (dependiendo de norm_chosen[i]), y los "datos reales" (en su caso x) y los parámetros se incluyen en alguna parte de este paso. Pero este valor de probabilidad logarítmica no está realmente conectado a ningún muestreador, sólo se calcula. Para que JAGS muestree realmente valores "buenos" de sus parámetros se requiere un paso más: el valor de probabilidad calculado se "ajusta" a un valor "observado" de Uno (es decir, el número 1) utilizando una distribución Bernoulli. Esto permite que los valores de probabilidad calculados más altos sean preferidos por el muestreador en virtud de tener una probabilidad más alta para un ensayo Bernoulli con un resultado conocido de 1. Por otro lado, si el vector Ones no se observa (lo que es perfectamente legal en la sintaxis de JAGS) entonces el resultado del ensayo Bernoulli no es conocido por JAGS, por lo que las probabilidades calculadas muy pequeñas son tan buenas como las probabilidades calculadas más grandes - es decir, su función de probabilidad personalizada es irrelevante. La principal desventaja de este truco (y del truco similar de Zeros en el que se utiliza una distribución Poisson en lugar del ensayo Bernoulli) es que JAGS se ve obligado a utilizar una estrategia de muestreo muy ineficiente.

Para ilustrarlo, vea los 3 modelos siguientes:

m1 <- 'model{

    for(i in 1:N){      
        Y[i] ~ dnorm(mu, tau)       
    }

    mu ~ dnorm(0, 10^-6)
    tau ~ dgamma(0.01, 0.01)

    #data# N, Y
    #monitor# mu, tau

}'

m2 <- 'model{

    for(i in 1:N){      
        ll[i] <- logdensity.norm(Y[i], mu, tau)
        Ones[i] ~ dbern(exp(ll[i]))
    }

    mu ~ dnorm(0, 10^-6)
    tau ~ dgamma(0.01, 0.01)

    #data# N, Y
    #monitor# mu, tau

}'

m3 <- 'model{

    for(i in 1:N){      
        ll[i] <- logdensity.norm(Y[i], mu, tau)
        Ones[i] ~ dbern(exp(ll[i]))
    }

    mu ~ dnorm(0, 10^-6)
    tau ~ dgamma(0.01, 0.01)

    #data# N, Y, Ones
    #monitor# mu, tau

}'

N <- 100
Y <- rnorm(N, 3, sqrt(1/2))

library('runjags')
run.jags(m1)
run.jags(m2)
Ones <- rep(1, N)
run.jags(m3)

Los modelos 1 y 3 son equivalentes, pero el modelo 1 es más eficiente (tamaño efectivo de la muestra de 20000 comparado con unos 12000 para los dos parámetros mu y tau): esta diferencia de eficiencia será más importante con modelos más complejos. El modelo 2 incluye los mismos "datos" Y pero omite los "Unos" y, por tanto, ofrece resultados que no tienen mucho sentido.