Eyes: normal mixture models

bowmaker:etal:85 analyse data on the peak sensitivity wavelengths for individual microspectrophotometric records on a small set of monkey's eyes. Data for one monkey (S14 in the paper) are given below (500 has been subtracted from each of the 48 measurements).

Part of the analysis involves fitting a mixture of two normal distributions with common variance to this distribution, so that each observation is assumed drawn from one of two groups. Let be the true group of the ith observation, where group j has a normal distribution with mean and precision . We assume an unknown fraction P of observations are in group 2, 1-P in group 1. The model is thus

We note that this formulation easily generalises to additional components to the mixture, although for identifiability an order constraint must be put onto the group means.

robert:94 points out that when using this model, there is a danger that at some iteration, all the data will go into one component of the mixture, and this state will be difficult to escape from -- this matches our experience. Robert suggests a re-parameterisation, a simplified version of which is to assume

are given independent ``noninformative" priors, including a uniform prior for P on (0,1). The appropriate graph is shown below, and the BUGS code is given over the page.

 
Figure 3:   Graphical model for eyes example

Eyes: model specification in BUGS

model eyes;
  
const
     N=48;
var
    y[N],        # observations
    T[N],        # true groups (labelled 1,2)
    lambda[2],   # means of two groups
    theta,       # scaled positive shift between groups
    tau,         # sampling precision
    sigma,       # sampling standard deviation
    P[2],        # proportion in first group
    alpha[2];    # prior parameters for proportions

data y in "eyes.dat";
inits in "eyes.in";
{
   for (i in 1:N){
       y[i]  ~ dnorm(lambda[T[i]],tau);       
       T[i]  ~ dcat(P[])
   }
   sigma     <- 1/sqrt(tau);
   tau        ~ dgamma(0.01,0.01);
   lambda[1]  ~ dnorm(0,1.0E-6); 
   lambda[2] <- lambda[1]+theta;
   theta      ~ dnorm(0,1.0E-6) I(0,);
   P[]        ~ ddirch(alpha[]);    # prior for mixing proportion
   alpha[1]  <- 1;                  # uniform prior   
   alpha[2]  <- 1;
}

Analysis

A BUGS run of 1000 iterations was extremely fast (4 seconds) and gave the following results which are compared with the maximum likelihood estimates of bowmaker:etal:85

We note the appropriately wider intervals provided by the full Bayesian analysis. We also point out that even with this re-parameterization we have experienced problems with some mixture models, in that a component may contain no observations at some iteration. One solution is to force at least one pre-specified observation to be in each component.