Eyes: Normal Mixture Model

Bowmaker et al (1985) 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).

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


   yi ~ Normal(λTi, τ)
   
   Ti ~ Categorical(P).


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 (1994) points out that when using this model, there is a danger that at some iteration, all the data will go into one component of themixture, and this state will be difficult to escape from --- this matches our experience. obert suggests a re-parameterisation, a simplified version of which is to assume

   λ2 = λ1+ θ, θ > 0.
   
λ1, θ, τ, P, are given independent ``noninformative" priors, including a uniform prior for P on (0,1). The appropriate graph and the BUGS code are given below.
[eyes2]


   model
   {
      for( i in 1 : N ) {
         y[i] ~ dnorm(mu[i], tau)
         mu[i] <- lambda[T[i]]
         T[i] ~ dcat(P[])
      }   
      P[1:2] ~ ddirich(alpha[])
      theta ~ dunif(0.0, 1000)
      lambda[2] <- lambda[1] + theta
      lambda[1] ~ dnorm(0.0, 1.0E-6)
      tau ~ dgamma(0.001, 0.001) sigma <- 1 / sqrt(tau)
   }

Data
list(y = c(529.0, 530.0, 532.0, 533.1, 533.4, 533.6, 533.7, 534.1, 534.8, 535.3,
         535.4, 535.9, 536.1, 536.3, 536.4, 536.6, 537.0, 537.4, 537.5, 538.3,
         538.5, 538.6, 539.4, 539.6, 540.4, 540.8, 542.0, 542.8, 543.0, 543.5,
         543.8, 543.9, 545.3, 546.2, 548.8, 548.7, 548.9, 549.0, 549.4, 549.9,
         550.6, 551.2, 551.4, 551.5, 551.6, 552.8, 552.9,553.2), N = 48, alpha = c(1, 1),
   T = c(1, NA, NA, NA, NA, NA, NA, NA, NA, NA,
            NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,         
            NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
            NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
            NA, NA, NA, NA, NA, NA, NA, 2))
Inits for chain 1
list(lambda = c(535, NA), theta = 5, tau = 0.1)
   
Inits for chain 2
list(lambda = c(100, NA), theta = 50, tau = 1)

Results
[eyes3]