Error when using bridge_sampler on a JAGS-Model with a latent mixture structure

[dizyd]

Hello,

first of all thank you for your work and the very helpful package you created. Now to my problem (this might be a pure user-related problem):

I get the following error when I run the bridge_sampler function based on samples from JAGS-Model with a latent mixture structure and I am not sure how to fix it:

Error in eigen(if (doDykstra) R else Y, symmetric = TRUE) : 
  infinite or missing values in 'x'
Additional warning:
effective sample size cannot be calculated, has been replaced by number of samples. 

The example model and data I used are

# Data
data <- list(y = rnorm(20,sample(c(0,10),20,T),3), N = 20)

# Model
model <- "model {

  for (i in 1:N) {

    y[i] ~ dnorm(mu[(z[i]+1)],tau)

    z[i] ~ dbern(phi)
  }

  ## Prior
  phi    ~ dbeta(1,1)
  mu[1]  ~ dnorm(0, .01)
  mu[2]  ~ dnorm(0, .01)
  tau    ~ dgamma(.01, .01)

}"

And the log_posterior() function is:

# log_posterior
log_posterior <- function(samples, data) {

    phi  <- samples[["phi"]]
  tau    <- samples[["tau"]]
  mu1  <- samples[["mu[1]"]]
  mu2  <- samples[["mu[2]"]]
  z       <- samples[,paste0("z[", seq_along(data$y), "]")]

  sum(dnorm(data$y, ifelse(z == 1,mu1,mu2), 1/sqrt(tau), log = TRUE)) +
    sum(dbinom(x = z, size = 1, prob = phi, log = TRUE)) +
    dgamma(tau,0.1,.01,log = TRUE) +
    dbeta(phi,1,1,log=TRUE) +
    dnorm(mu1,0,.01,log=TRUE) +
    dnorm(mu2,0,.01,log=TRUE)

}

The problem has to do with the z variable, since the functions work when I replaceifelse(z == 1,mu1,mu2) in log_posterior() with e.g. ifelse(sample(c(0,1),0.5) == 1,mu1,mu2).

I haven't found anything in the manual or other issues regarding this error. I already tried running more samples and chains as recommended in another Issue. Therefore, my question, if you have any idea what I am doing wrong?

The full example script is:

# Load packages
library(tidyverse)
library(bridgesampling)
library(R2jags)
library(runjags)

# Data
data <- list(y = rnorm(20,sample(c(0,10),20,T),3), N = 20)

# Model
model <- "model {

  for (i in 1:N) {

    y[i] ~ dnorm(mu[(z[i]+1)],tau)

    z[i] ~ dbern(phi)
  }

  ## Prior
  phi    ~ dbeta(1,1)
  mu[1]  ~ dnorm(0, .01)
  mu[2]  ~ dnorm(0, .01)
  tau    ~ dgamma(.01, .01)

}"

# With z    ------------------------------------------------------------------

# Samples
getSamples_z <- function(data,model){

  s <- jags(data,
            parameters.to.save = c("mu","tau","z","phi"),
            model.file = textConnection(model),
            n.chains = 3, n.iter = 52000,
            n.burnin = 2000, n.thin = 1)

  return(s)

}  

samples_z    <- getSamples_z(data,model)

# log_posterior
log_posterior_z <- function(samples, data) {

  phi    <- samples[["phi"]]
  tau    <- samples[["tau"]]
  mu1    <- samples[["mu[1]"]]
  mu2    <- samples[["mu[2]"]]
  z      <- samples[,paste0("z[", seq_along(data$y), "]")]

  sum(dnorm(data$y, ifelse(z == 1,mu1,mu2), 1/sqrt(tau), log = TRUE)) +
    sum(dbinom(x = z, size = 1, prob = phi, log = TRUE)) +
    dgamma(tau,0.1,.01,log = TRUE) +
    dbeta(phi,1,1,log=TRUE) +
    dnorm(mu1,0,.01,log=TRUE) +
    dnorm(mu2,0,.01,log=TRUE)

}

# Bridge Sampling

cn    <- colnames(samples_z$BUGSoutput$sims.matrix)
cn    <- cn[cn != "deviance"]
lb_z  <- c(-Inf,-Inf,0, 0, rep(0,data$N))
ub_z  <- c(Inf,Inf,1,Inf,rep(1,data$N))
names(lb_z) <- names(ub_z) <- cn

# compute log marginal likelihood via bridge sampling for H0
Mbridge <- bridge_sampler(samples       = samples_z,
                          data          = data,
                          log_posterior = log_posterior_z,
                          lb            = lb_z,
                          ub            = ub_z,
                          silent        = TRUE)

[quentingronau]

Hi @dizyd,

The reason why this does not work is that your model contains discrete variables (i.e., the z's). To use bridge sampling in this case, you need to marginalize out the discrete parameters in the log_posterior function yourself and then remove the z's from the samples object before you pass the samples to the bridge_sampler function. This is mentioned in our JSS paper (p. 22). One explanation of how to marginalize out the mixture indicator variable can be found here.

Cheers, Quentin

[dizyd]

Hello @quentingronau,

thank you very much for your fast response and the references!

Cheers, David