Probabilidad marginal de la salida de Gibbs

Estoy reproduciendo desde cero los resultados en la Sección 4.2.1 de

Probabilidad marginal de la salida de Gibbs

Siddhartha Chib

Revista de la Asociación Americana de Estadística, vol. 90, núm. 432. (diciembre de 1995), págs. 1313-1321.

Es una mezcla de modelos normales con el número conocido de componentes. $k\geq 1$

$$f(x\mid w,\mu,\sigma^2) =\prod_{i=1}^n\sum_{j=1}^k \mathrm{N}(x_i\mid\mu_j,\sigma_j^2) \, . \qquad (*)$$

La muestra de Gibbs para este modelo se implementa utilizando la técnica de aumento de datos de Tanner y Wong. Un conjunto de variables de asignación suponiendo que los valores 1, \ dots, k se introduce, y se especifica que \ Pr (Z_I = j \ mediados w) = w_j y f (x_i \ mediados z , \ mu, \ sigma ^ 2) = \ mathrm {N} (x_i \ mid \ mu_ {z_i}, \ sigma ^ 2_ {z_i}) . De ello se deduce que la integración sobre el z_i 's da la probabilidad original (*) .$z=(z_1,\dots,z_n)$$1,\dots,k$$\Pr(z_i=j\mid w)=w_j$z i ( ∗ $f(x_i\mid z,\mu,\sigma^2)=\mathrm{N}(x_i\mid\mu_{z_i},\sigma^2_{z_i})$$z_i$$(*)$

El conjunto de datos está formado por velocidades de $82$ galaxias de la constelación Corona Borealis.

set.seed(1701)

x <- c(  9.172,  9.350,  9.483,  9.558,  9.775, 10.227, 10.406, 16.084, 16.170, 18.419, 18.552, 18.600, 18.927,
        19.052, 19.070, 19.330, 19.343, 19.349, 19.440, 19.473, 19.529, 19.541, 19.547, 19.663, 19.846, 19.856,
        19.863, 19.914, 19.918, 19.973, 19.989, 20.166, 20.175, 20.179, 20.196, 20.215, 20.221, 20.415, 20.629,
        20.795, 20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814, 21.921, 21.960, 22.185, 22.209,
        22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888, 22.914, 23.206, 23.241, 23.263, 23.484,
        23.538, 23.542, 23.666, 23.706, 23.711, 24.129, 24.285, 24.289, 24.366, 24.717, 24.990, 25.633, 26.960,
        26.995, 32.065, 32.789, 34.279 )

nn <- length(x)

Asumimos que , 's y ' s son independientes a priori con μ j σ 2 j ( w 1 , … , w k ) ∼ D i r ( a 1 , … , a k$w$$\mu_j$$\sigma^2_j$

$$(w_1,\dots,w_k) \sim \mathrm{Dir}(a_1,\dots,a_k) \, , \quad \mu_j \sim \mathrm{N}(\mu_0,\sigma_0^2) \, , \quad \sigma^2_j\sim\mathrm{IG}\!\left(\frac{\nu_0}{2},\frac{\delta_0}{2}\right) \, .$$

k <- 3

mu0 <- 20
va0 <- 100

nu0 <- 6
de0 <- 40

a <- rep(1, k)

Usando el teorema de Bayes, los condicionales completos son en el que con

$$\begin{align*} w \mid \mu,\sigma^2,z,x &\sim \mathrm{Dir}(a_1+n_1,\dots,a_k+n_k) \\ \mu_j \mid w, \sigma^2,z,x &\sim \mathrm{N}\!\left( \frac{n_j m_j\sigma_0^2+\mu_0\sigma_j^2}{n_j\sigma^2_0+\sigma^2_j}, \frac{\sigma^2_0\sigma^2_j}{n_j\sigma^2_0+\sigma^2_j}\right) \\ \sigma_j^2 \mid w,\mu,z,x &\sim \mathrm{IG}\!\left( \frac{\nu_0+n_j}{2},\frac{\delta_0+\delta_j}{2}\right) \\ \Pr(z_i=j\mid w,\mu,\sigma^2,x) &\propto w_j \times \frac{1}{\sigma_j}e^{-(x_i-\mu_j)^2/2\sigma_j^2} \end{align*}$$ $$n_j = |L_j| \, , \qquad m_j = \begin{cases}\frac{1}{n_j}\sum_{i\in L_j} x_i &\;\mathrm{if}\; n_j>0 \\ 0 &\;\mathrm{otherwise.} \end{cases}\, , \qquad \delta_j = \sum_{i\in L_j} (x_i-\mu_j)^2 \, ,$$ $L_j=\{i\in\{1,\dots,n\}:z_i=j\}$ .

El objetivo es calcular una estimación de la probabilidad marginal del modelo. El método de Chib comienza con una primera ejecución de la muestra de Gibbs utilizando los condicionales completos.

burn_in <- 1000
run     <- 15000

cat("First Gibbs run (full):\n")

N <- burn_in + run

w  <- matrix(1, nrow = N, ncol = k)
mu <- matrix(0, nrow = N, ncol = k)
va <- matrix(1, nrow = N, ncol = k)
z  <- matrix(1, nrow = N, ncol = nn)

n <- integer(k)
m <- numeric(k)
de <- numeric(k)

rdirichlet <- function(a) { y <- rgamma(length(a), a, 1); y / sum(y) }

pb <- txtProgressBar(min = 2, max = N, style = 3)
z[1,] <- sample.int(k, size = nn, replace = TRUE)
for (t in 2:N) {
    n <- tabulate(z[t-1,], nbins = k)
    w[t,] <- rdirichlet(a + n)
    m <- sapply(1:k, function(j) sum(x[z[t-1,]==j]))
    m[n > 0] <- m[n > 0] / n[n > 0]
    mu[t,] <- rnorm(k, mean = (n*m*va0+mu0*va[t-1,])/(n*va0+va[t-1,]), sd = sqrt(va0*va[t-1,]/(n*va0+va[t-1,])))
    de <- sapply(1:k, function(j) sum((x[z[t-1,]==j] - mu[t,j])^2))
    va[t,] <- 1 / rgamma(k, shape = (nu0+n)/2, rate = (de0+de)/2)
    z[t,] <- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mu[t,], sd = sqrt(va[t,]), log = TRUE))))
    setTxtProgressBar(pb, t)
}
close(pb)

De esta primera ejecución obtenemos un punto aproximado de máxima probabilidad. Dado que la probabilidad es realmente ilimitada, lo que probablemente da este procedimiento es un MAP local aproximado.$(w^*,\mu^*,\sigma^{2*})$

w  <- w[(burn_in+1):N,]
mu <- mu[(burn_in+1):N,]
va <- va[(burn_in+1):N,]
z  <- z[(burn_in+1):N,]
N  <- N - burn_in

log_L <- function(x, w, mu, va) sum(log(sapply(1:nn, function(i) sum(exp(log(w) + dnorm(x[i], mean = mu, sd = sqrt(va), log = TRUE))))))

ts <- which.max(sapply(1:N, function(t) log_L(x, w[t,], mu[t,], va[t,])))

ws <- w[ts,]
mus <- mu[ts,]
vas <- va[ts,]

La estimación de registro de Chib de la probabilidad marginal es

$$\begin{align} \log \widehat{f(x)} &= \log L_x(w^*,\mu^*,\sigma^{2*}) + \log \pi(w^*,\mu^*,\sigma^{2*}) \\ &- \log \pi(\mu^*\mid x) - \log \pi(\sigma^{2*}\mid \mu^*,x) - \log \pi(w^*\mid \mu^*,\sigma^{2*},x) \, . \end{align}$$

Ya tenemos los dos primeros términos.

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w))
    + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE))
    + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

chib <- log_L(x, ws, mus, vas) + log_prior(ws, mus, vas)

La estimación Rao-Blackwellized de es y se obtiene fácilmente desde la primera carrera de Gibbs.$\pi(\mu^*\mid x)$

$$\pi(\mu^*\mid x) = \int \prod_{j=1}^k \mathrm{N}\!\left(\mu_j^* \;\Bigg|\; \frac{n_j m_j\sigma_0^2+\mu_0\sigma_j^2}{n_j\sigma^2_0+\sigma^2_j}, \frac{\sigma^2_0\sigma^2_j}{n_j\sigma^2_0+\sigma^2_j}\right)\,p(\sigma^{2},z\mid x)\,d\sigma^2\,dz \, ,$$

pi.mu_va.z.x <- function(mu, va, z) {
    n <- tabulate(z, nbins = k)
    m <- sapply(1:k, function(j) sum(x[z==j]))
    m[n > 0] <- m[n > 0] / n[n > 0]
    exp(sum(dnorm(mu, mean = (n*m*va0+mu0*va)/(n*va0+va), sd = sqrt(va0*va/(n*va0+va)), log = TRUE)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.mu_va.z.x(mus, va[t,], z[t,]))))

La estimación Rao-Blackwellized de es y se calcula a partir de una segunda ejecución reducida de Gibbs en la que los no se actualizan, sino que se hacen igual a en cada paso de iteración.$\pi(\sigma^{2*}\mid \mu^*,x)$

$$\pi(\sigma^{2*}\mid \mu^*,x) = \int \prod_{j=1}^k \mathrm{IG}\!\left( \sigma^{2*}_j \;\Bigg|\; \frac{\nu_0+n_j}{2},\frac{\delta_0+\delta_j}{2}\right) \, p(z\mid\mu^*,x)\,dz \, ,$$ $\mu_j$$\mu^*_j$

cat("Second Gibbs run (reduced):\n")

N <- burn_in + run

w  <- matrix(1, nrow = N, ncol = k)
va <- matrix(1, nrow = N, ncol = k)
z  <- matrix(1, nrow = N, ncol = nn) 

pb <- txtProgressBar(min = 2, max = N, style = 3)
z[1,] <- sample.int(k, size = nn, replace = TRUE)
for (t in 2:N) {
    n <- tabulate(z[t-1,], nbins = k)
    w[t,] <- rdirichlet(a + n)
    de <- sapply(1:k, function(j) sum((x[z[t-1,]==j] - mus[j])^2))
    va[t,] <- 1 / rgamma(k, shape = (nu0+n)/2, rate = (de0+de)/2)
    z[t,] <- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mus, sd = sqrt(va[t,]), log = TRUE))))
    setTxtProgressBar(pb, t)
}
close(pb)

w  <- w[(burn_in+1):N,]
va <- va[(burn_in+1):N,]
z  <- z[(burn_in+1):N,]
N  <- N - burn_in

pi.va_mu.z.x <- function(va, mu, z) {
    n <- tabulate(z, nbins = k)         
    de <- sapply(1:k, function(j) sum((x[z==j] - mu[j])^2))
    exp(sum(((nu0+n)/2)*log((de0+de)/2) - lgamma((nu0+n)/2) - ((nu0+n)/2+1)*log(va) - (de0+de)/(2*va)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.va_mu.z.x(vas, mus, z[t,]))))

Del mismo modo, la estimación Rao-Blackwellized de es y se calcula a partir de una tercera ejecución reducida de Gibbs en la que los 's y los ' s no se actualizan, sino que se igualan a y respectivamente en cada paso de iteración.$\pi(w^*\mid \mu^*,\sigma^{2*},x)$μ j σ 2 j μ ∗ j σ 2 ∗

$$\pi(w^*\mid \mu^*,\sigma^{2*},x) = \int \mathrm{Dir}(w^* \mid a_1+n_1,\dots,a_k+n_k) \, p(z\mid\mu^*,\sigma^{2*},x)\,dz \, ,$$ $\mu_j$$\sigma^2_j$$\mu^*_j$$\sigma^{2*}_j$

cat("Third Gibbs run (reduced):\n")

N <- burn_in + run

w  <- matrix(1, nrow = N, ncol = k)
z  <- matrix(1, nrow = N, ncol = nn) 

pb <- txtProgressBar(min = 2, max = N, style = 3)
z[1,] <- sample.int(k, size = nn, replace = TRUE)
for (t in 2:N) {
    n <- tabulate(z[t-1,], nbins = k)
    w[t,] <- rdirichlet(a + n)
    z[t,] <- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mus, sd = sqrt(vas), log = TRUE))))
    setTxtProgressBar(pb, t)
}
close(pb)

w  <- w[(burn_in+1):N,]
z  <- z[(burn_in+1):N,]
N  <- N - burn_in

pi.w_z.x <- function(w, z) {
    n <- tabulate(z, nbins = k)
    exp(lgamma(sum(a+n)) - sum(lgamma(a+n)) + sum((a+n-1)*log(w)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.w_z.x(ws, z[t,]))))

Después de todo esto, obtenemos una estimación que es más grande que la informada por Chib: con el error de Monte Carlo .- 224.138 $-217.9199$$-224.138$$.086$

Para comprobar si de alguna manera estropeé los muestreadores Gibbs, reimplementé todo usando RJAGS. El siguiente código da los mismos resultados.

x <- c( 9.172,  9.350,  9.483,  9.558,  9.775, 10.227, 10.406, 16.084, 16.170, 18.419, 18.552, 18.600, 18.927, 19.052, 19.070, 19.330,
       19.343, 19.349, 19.440, 19.473, 19.529, 19.541, 19.547, 19.663, 19.846, 19.856, 19.863, 19.914, 19.918, 19.973, 19.989, 20.166,
       20.175, 20.179, 20.196, 20.215, 20.221, 20.415, 20.629, 20.795, 20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814,
       21.921, 21.960, 22.185, 22.209, 22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888, 22.914, 23.206, 23.241, 23.263,
       23.484, 23.538, 23.542, 23.666, 23.706, 23.711, 24.129, 24.285, 24.289, 24.366, 24.717, 24.990, 25.633, 26.960, 26.995, 32.065,
       32.789, 34.279 )

library(rjags)

nn <- length(x)

k <- 3

mu0 <- 20
va0 <- 100

nu0 <- 6
de0 <- 40

a <- rep(1, k)

burn_in <- 10^3

N <- 10^4

full <- "
    model {
        for (i in 1:n) {
            x[i] ~ dnorm(mu[z[i]], tau[z[i]])
            z[i] ~ dcat(w[])
        }
        for (i in 1:k) {
            mu[i] ~ dnorm(mu0, 1/va0)
            tau[i] ~ dgamma(nu0/2, de0/2)
            va[i] <- 1/tau[i]
        }
        w ~ ddirich(a)
    }
"
data <- list(x = x, n = nn, k = k, mu0 = mu0, va0 = va0, nu0 = nu0, de0 = de0, a = a)
model <- jags.model(textConnection(full), data = data, n.chains = 1, n.adapt = 100)
update(model, n.iter = burn_in)
samples <- jags.samples(model, c("mu", "va", "w", "z"), n.iter = N)

mu <- matrix(samples$mu, nrow = N, byrow = TRUE)
    va <- matrix(samples$va, nrow = N, byrow = TRUE)
w <- matrix(samples$w, nrow = N, byrow = TRUE)
    z <- matrix(samples$z, nrow = N, byrow = TRUE)

log_L <- function(x, w, mu, va) sum(log(sapply(1:nn, function(i) sum(exp(log(w) + dnorm(x[i], mean = mu, sd = sqrt(va), log = TRUE))))))

ts <- which.max(sapply(1:N, function(t) log_L(x, w[t,], mu[t,], va[t,])))

ws <- w[ts,]
mus <- mu[ts,]
vas <- va[ts,]

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w))
    + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE))
    + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

chib <- log_L(x, ws, mus, vas) + log_prior(ws, mus, vas)

cat("log-likelihood + log-prior =", chib, "\n")

pi.mu_va.z.x <- function(mu, va, z, x) {
    n <- sapply(1:k, function(j) sum(z==j))
    m <- sapply(1:k, function(j) sum(x[z==j]))
    m[n > 0] <- m[n > 0] / n[n > 0]
    exp(sum(dnorm(mu, mean = (n*m*va0+mu0*va)/(n*va0+va), sd = sqrt(va0*va/(n*va0+va)), log = TRUE)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.mu_va.z.x(mus, va[t,], z[t,], x))))

cat("log-likelihood + log-prior - log-pi.mu_ =", chib, "\n")

fixed.mu <- "
    model {
        for (i in 1:n) {
            x[i] ~ dnorm(mus[z[i]], tau[z[i]])
            z[i] ~ dcat(w[])
        }
        for (i in 1:k) {
            tau[i] ~ dgamma(nu0/2, de0/2)
            va[i] <- 1/tau[i]
        }
        w ~ ddirich(a)
    }
"
data <- list(x = x, n = nn, k = k, nu0 = nu0, de0 = de0, a = a, mus = mus)
model <- jags.model(textConnection(fixed.mu), data = data, n.chains = 1, n.adapt = 100)
update(model, n.iter = burn_in)
samples <- jags.samples(model, c("va", "w", "z"), n.iter = N)

va <- matrix(samples$va, nrow = N, byrow = TRUE)
    w <- matrix(samples$w, nrow = N, byrow = TRUE)
z <- matrix(samples$z, nrow = N, byrow = TRUE)

pi.va_mu.z.x <- function(va, mu, z, x) {
    n <- sapply(1:k, function(j) sum(z==j))
    de <- sapply(1:k, function(j) sum((x[z==j] - mu[j])^2))
    exp(sum(((nu0+n)/2)*log((de0+de)/2) - lgamma((nu0+n)/2) - ((nu0+n)/2+1)*log(va) - (de0+de)/(2*va)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.va_mu.z.x(vas, mus, z[t,], x))))

cat("log-likelihood + log-prior - log-pi.mu_ - log-pi.va_ =", chib, "\n")

fixed.mu.and.va <- "
    model {
        for (i in 1:n) {
            x[i] ~ dnorm(mus[z[i]], 1/vas[z[i]])
            z[i] ~ dcat(w[])
        }
        w ~ ddirich(a)
    }
"
data <- list(x = x, n = nn, a = a, mus = mus, vas = vas)
model <- jags.model(textConnection(fixed.mu.and.va), data = data, n.chains = 1, n.adapt = 100)
update(model, n.iter = burn_in)
samples <- jags.samples(model, c("w", "z"), n.iter = N)

w <- matrix(samples$w, nrow = N, byrow = TRUE)
    z <- matrix(samples$z, nrow = N, byrow = TRUE)

pi.w_z.x <- function(w, z, x) {
    n <- sapply(1:k, function(j) sum(z==j))
    exp(lgamma(sum(a)+nn) - sum(lgamma(a+n)) + sum((a+n-1)*log(w)))
}

chib <- chib - log(mean(sapply(1:N, function(t) pi.w_z.x(ws, z[t,], x))))

cat("log-likelihood + log-prior - log-pi.mu_ - log-pi.va_ - log-pi.w_ =", chib, "\n")

Mi pregunta es si en la descripción anterior hay algún malentendido del método de Chib o algún error en su implementación.

Hay un pequeño error de programación en el anterior

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w))
    + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE))
    + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

como debería ser en su lugar

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w)) +
      sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE)) +
      sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

Volver a ejecutar el código de esta manera conduce a

> chib
[1] -228.194

¡cuál no es el valor producido en Chib (1995) para ese caso! Sin embargo, en el reanálisis del problema de Neal (1999), menciona que

Según un árbitro anónimo de JASA, la cifra de -224.138 para el registro de la probabilidad marginal para el modelo de tres componentes con variaciones desiguales que se dio en el documento de Chib es un "error tipográfico" con la cifra correcta siendo -228.608.

Entonces esto resuelve el problema de discrepancia.