#>source("c:/Work07/ShortCourse07/Programs/NormalMixtureGibbs.txt",print.eval=TRUE)#
#Gibbs Sampler for a normal mixture#
#--------------Initial values-------------------------------
mu1<--1			#True Values
mu2<-1			#True Values
mu1start<-mu2		#Start of Gibbs Chain
mu2start<-mu2		#Start of Gibbs Chain
sigma<-2
c<-10
p<-.25
nsim<-2000
mG1<-array(mu1start,c(nsim,1))
mG2<-array(mu2start,c(nsim,1))
#----------------------Data------------------------------
n<-25
xdata<-mu1+(runif(n)>p)*(mu2-mu1)+rnorm(n,mean=0,sd=sigma)
Z<-rbinom(n,1,p)
#------------------Gibbs Sampler--------------------------
for(i in 2:nsim)
{
mG1[i]<-rnorm(1,mean=sum(xdata*Z)/((1/c)+sum(Z)),sd=sigma/sqrt((1/c)+sum(Z)))
mG2[i]<-rnorm(1,mean=sum(xdata*(1-Z))/((1/c)+sum(1-Z)),sd=sigma/sqrt((1/c)+sum(1-Z)))
for(j in 1:n)
	{
	 p1<- p*dnorm(xdata[j],mG1[i],sigma)
	 p2<- (1-p)*dnorm(xdata[j],mG2[i],sigma)
	 Z[j]<-rbinom(1,1,p1/(p1+p2))
	}
}
#-----------------------------Plot-----------------------------
par(mfrow=c(1,3))
plot(function(x)p*dnorm(x,mu1,sigma)+(1-p)*dnorm(x,mu2,sigma),mu1-3*sigma,mu2+3*sigma,ylab="",lwd=2)
plot(mG1,mG2,type="l",xlab="mu1",ylab="mu2")
yrange<-c(mu1-2,mu2+2)
plot(cumsum(mG1)/c(1:nsim),type="l",ylim=yrange,ylab="",lwd=2)
par(new=T)
plot(cumsum(mG2)/c(1:nsim),type="l",ylim=yrange,ylab="",lwd=2)