################################################################################################
#
# WINBUGS EXAMPLE (VOLUME II) - 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).    
#
#	29.0   30.0   32.0   33.1   33.4   33.6   33.7   34.1   34.8   35.3
#	35.4   35.9   36.1   36.3   36.4   36.6   37.0   37.4   37.5   38.3
#	38.5   38.6   39.4   39.6   40.4   40.8   42.0   42.8   43.0   43.5
#	43.8   43.9   45.3   46.2   48.8   48.7   48.9   49.0   49.4   49.9
#	50.6   51.2   51.4   51.5   51.6   52.8   52.9   53.2
#
#		
# 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 ith observation, where 
# group j has a normal distribution with mean lambda[j] and precision tau.
#
# We assume an unknown fraction P of observations are in group 2, 1-P in group 1. 
# The model is thus 
#
#	yi  ~ Normal(lambda[T[i]],tau)
#	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 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 
#
#	lambda[2]= lambdad[1] + theta,  theta  >  0.
#	
# lambda[1], theta, tau, P are given independent ``noninformative" priors, including a 
# uniform prior for P on (0,1).  
#
################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
################################################################################################
source("mixture-routines.R")

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 = length(y)

#n = 200
#z = sample(1:2,size=n,replace=TRUE,prob=c(0.7,0.3))
#means=c(537,537+25)
#vars =c(22,5)
#y = rnorm(n,means[z],sqrt(vars[z]))

par(mfrow=c(1,1))
hist(y,nclass=15)

# hyperparameters
alpha     = 0.001
beta      = 0.001
lambda0   = 0
C0        = 100000
theta0    = 0
V0        = 100000

# initial values
lambda    = rep(0,2)
lambda[1] = 535
theta     = 5
lambda[2] = lambda[1]+theta
sigma2    = 10
p         = 0.7
T         = rep(1,n)

# MCMC setup
burnin = 1000
M      = 1000

# Common variance
set.seed(25879)
output1 = MCMC.commonvariance(y,alpha,beta,lambda0,C0,theta0,V0,p,lambda,sigma2,burnin,M)

# Different variances
alpha     = c(0.001,0.001)
beta      = c(0.001,0.001)
sigma2    = c(10,10)
output2 = MCMC.twovariances(y,alpha,beta,lambda0,C0,theta0,V0,p,lambda,sigma2,burnin,M)

par(mfrow=c(2,4))
# printing output from MCMC scheme with common variances
ts.plot(output1$p,xlab="",ylab="",main="P \n common variance",ylim=c(0,1))
lines(1-output1$p,col=2)
L = min(output1$lambda,output1$lambda+output1$theta)
U = max(output1$lambda,output1$lambda+output1$theta)
ts.plot(output1$lambda,xlab="",ylab="",main="Lambda1 and lambda2 \n common variance",ylim=c(L,U))
lines(output1$lambda+output1$theta,col=2)
ts.plot(output1$sigma2,xlab="",ylab="",main="Sigma \n common variance")
plot(y,output1$mT,ylab="Pr(Z=1)",xlab="data")
abline(h=0.5,col=2,lwd=2)
title(paste("cutoff=",mean(max(y[output1$mT<0.5]),min(y[output1$mT>0.5])),sep=""))

# printing output from MCMC scheme with specific variances
ts.plot(output2$p,xlab="",ylab="",main="P \n specific variances",ylim=c(0,1))
lines(1-output2$p,col=2)
L = min(output2$lambda,output2$lambda+output2$theta)
U = max(output2$lambda,output2$lambda+output2$theta)
ts.plot(output2$lambda,xlab="",ylab="",main="Lambda1 and lambda2 \n specific variances",ylim=c(L,U))
lines(output2$lambda+output2$theta,col=2)
L = min(output2$sigma21,output2$sigma22)
U = max(output2$sigma21,output2$sigma22)
ts.plot(output2$sigma21,xlab="",ylab="",main="Sigma1 \n specific variances",ylim=c(L,U))
lines(output2$sigma22,col=2)
plot(y,output2$mT,ylab="Pr(Z=1)",xlab="data")
abline(h=0.5,col=2,lwd=2)
title(paste("cutoff=",mean(max(y[output2$mT<0.5]),min(y[output2$mT>0.5])),sep=""))


par(mfrow=c(2,3))
# printing output from MCMC scheme with common variances
plot(density(output1$p),xlab="",ylab="",main="P \n common variance",xlim=c(0,1))
L = min(output1$lambda,output1$lambda+output1$theta)
U = max(output1$lambda,output1$lambda+output1$theta)
plot(density(output1$lambda),xlab="",ylab="",main="Lambda1 and lambda2 \n common variance",xlim=c(L,U))
lines(density(output1$lambda+output1$theta),col=2)
plot(density(output1$sigma2),xlab="",ylab="",main="Sigma \n common variance")

# printing output from MCMC scheme with specific variances
plot(density(output2$p),xlab="",ylab="",main="P \n specific variances",xlim=c(0,1))
L = min(output2$lambda,output2$lambda+output2$theta)
U = max(output2$lambda,output2$lambda+output2$theta)
plot(density(output2$lambda+output2$theta),xlab="",ylab="",main="Lambda1 and lambda2 \n common variance",xlim=c(L,U),col=2)
lines(density(output2$lambda),col=1)
plot(density(output2$sigma22),xlab="",ylab="",main="Sigma \n common variance",col=2)
lines(density(output2$sigma21),col=1)


N = 100
yy = seq(min(y),max(y),length=N)
pred1 = NULL
pred2 = NULL
for (i in 1:N){
  pred1 = c(pred1,mean(output1$p*dnorm(yy[i],output1$lambda1,output1$sigma2)+
                    (1-output1$p)*dnorm(yy[i],output1$lambda1+output1$theta,output1$sigma2))
  pred2 = c(pred2,mean(output2$p*dnorm(yy[i],output2$lambda1,output2$sigma21)+
                    (1-output2$p)*dnorm(yy[i],output2$lambda1+output2$theta,output2$sigma22))
}