How to code a simple latent mixture model in JAGS

  • Nathaniel Phillips

    I am struggling with how to program a simple latent mixture model in JAGS. I have a vector of 10 data points and I want to model the data as either coming from a normal distribution or a beta distribution. I want to simultaneously estimate the posterior densities of the parameters of each model, while estimating the probability that the data came from each model:

    Ideally, I would like to model the data as:

    Data[i] ~ dt(mu, sd, df)*equals(Model.Index, 0) + dbeta(alpha, beta) * equals(Model.Index, 1)

    Where Model.Index is a bernoulli indicator variable. But this syntax is not JAGS friendly.

    I have included sample code (from R) below. If someone could help me fix the code such that it would allow me to estimate both model parameters and the probability of each model that would be great!

    Start Sample R Code


    Data <- rnorm(10, .5, .1)

    Model.String <- "

    for (i in 1:10) {

    # This is what I want to do: calculate the likelihood of the data depending on whether Model = 0 (the student t model), or Model = 1 (the beta model)

    # Data[i] ~ dt(mu, sd, df)*equals(Model.Index, 0) + dbeta(alpha, beta) * equals(Model.Index, 1)

    # But JAGS doesn't like this

    # These work individually, but How do I get JAGS to estimate the probability of either model?
    Data[i] ~ dt(mu, sd, df)T(0,1)
    # Data[i] ~ dbeta(alpha, beta)


    Model.Index ~ dbern(.5)
    mu ~ dnorm(0, .001)
    sd ~ dgamma(.001, .001)
    df ~ dgamma(.001, .001)
    alpha ~ dunif(0, 100)
    beta ~ dunif(0, 100)


    writeLines(Model.String, con = "model.txt")
    data <- list("Data") <- c("mu", "sd", "df", "Model.Index", "alpha", "beta")

    samples <- jags(data, =, model.file = "model.txt", n.chains = 1, n.iter = 100, n.burnin = 1, n.thin = 1, DIC = T)

    End Sample R code
  • Russell Almond

    Russell Almond - 2014-03-26

    I'm a little bit confused by the difference between your code and your description. From your code, it doesn't look like you are really doing a mixture model (where Model.Index would be Model.Index[i]) but rather trying to fit two different models and then compute the posterior probablity of the two parameterizations. If that is your intent, you might be better off first computing the posteriors for each of the two proposed models and then the modelling probabilities. Your two submodels are conjugate, so you don't really need JAGS. David Draper has written a lot on this subject.

    I suppose you could also monitor the Deviance and work with that, but as the parameters of the two models are really independent given the model.index, the deviances from the two models at given cycles will have nothing to do with each other. In fact, the joint MCMC is probably less efficient at estimating the posterior probabilities of the two models than estimating the deviance separately for each one and then doing the model averaging step in R.

    If you really mean a mixture model, the way I've been doing it is to introduce unobservable variables which is what the observation would be under each model and a switching parameter.

    I would use something like this:

    model {
    ## Hyperparameters
    pi ~ dbeta(.5,.5)
    mu ~ dnorm(0, .001)
    sd ~ dgamma(.001, .001)
    df ~ dgamma(.001, .001)
    alpha ~ dunif(0, 100)
    beta ~ dunif(0, 100)
    ## Observation specific
    for (i in 1:I) {
      Model.index[i] ~dbern(pi)
      Latent[i,1] ~ dt(mu, sd, df)T(0,1)
      Latent[i,2] ~ dbeta(alpha,beta)
      Data[i] <- Latent[i,Model.index[i]+1] ## Index is 1 based, Bernouli is 0 base.

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.

No, thanks