set.seed(32847)
n<-100
omega<-0.25
muvec<-c(-3,4)
Z<-rbinom(n,1,omega)
Y<-(1-Z)*rnorm(n,muvec[Z+1],1)+Z*rnorm(n,muvec[Z+1],1)
hist(Y)

yvec<-seq(-10,10,by=0.001)
dvec<-(1-omega)*dnorm(yvec,muvec[1],1)+omega*dnorm(yvec,muvec[2],1)
plot(yvec,dvec,type="l")

nburn<-10
nsamp<-2000
nthin<-10
nits<-nburn+nthin*nsamp

#Ordinary Metropolis
#Try different values of sigma.q=0.5,1,2
sigma.q<-0.5
old.y<--5
old.val<-log((1-omega)*dnorm(old.y,muvec[1],1)+omega*dnorm(old.y,muvec[2],1))
plot(yvec,dvec,type="l",xlab="y")
ico<-0;Ysamp<-rep(0,nsamp)
set.seed(3447)
t0<-proc.time()[3]
for(iter in 1:nits){

	new.y<-old.y+rnorm(1)*sigma.q
	new.val<-log((1-omega)*dnorm(new.y,muvec[1],1)+omega*dnorm(new.y,muvec[2],1))
	if(log(runif(1))< new.val-old.val){
		old.y<-new.y
		old.val<-new.val
	}
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		Ysamp[ico]<-old.y
		points(old.y,0,pch=19)
	}

}
t1<-proc.time()[3]
print(t1-t0)

par(mfrow=c(1,2))
hist(Ysamp,breaks=seq(-10,10,by=0.5),ylim=range(0,500),main="Sample from Mixture Distribution")
box()
lines(yvec,dvec*nsamp*0.5,type="l",xlab="y")
plot(Ysamp,ylim=range(-10,10),ylab="Y")

#Does not always sample from both modes

###############################################################################
#Try a population method with M=3 chains
M<-3
Temps<-c(1,5,20)

#Compute the normalizing constants

powered.dens<-function(yv,muv,omv,Tv){
	dv<-exp(log((1-omv)*dnorm(yv,muv[1],1)+omv*dnorm(yv,muv[2],1))/Tv)
	return(dv)
}

tempered.distribution<-function(muv,omv,lv=-100,uv=100,Tv){

	val<-integrate(f=powered.dens,lower=lv,upper=uv,muv=muv,omv=omv,Tv=Tv)
	return(val$val)

}


Norm<-rep(0,M)
yv<-seq(-15,15,by=0.01)
par(mfrow=c(1,1))
plot(yv,powered.dens(yv,muvec,omega,Tv=Temps[1]),type="l",ylab="Density")
stored.dv<-matrix(0,nrow=length(yv),ncol=M)
for(j in 1:M){
	Norm[j]<-tempered.distribution(muvec,omega,Tv=Temps[j])
	dv<-powered.dens(yv,muvec,omega,Tv=Temps[j])
	lines(yv,dv/Norm[j])
	stored.dv[,j]<-dv/Norm[j]
}



nburn<-10
nsamp<-2000
nthin<-10
nits<-nburn+nthin*nsamp

#Population Metropolis
#Try different values of sigma.q=0.5,1,2

sigma.q<-0.5
old.yvec<-new.yvec<-rep(-5,3)
old.vals<-new.vals<-rep(0,3)
for(j in 1:3){
	old.vals[j]<-log(powered.dens(old.yvec[j],muvec,omega,Tv=Temps[j]))-log(Norm[j])
}

ico<-0;Ysamp.pop<-matrix(0,nrow=nsamp,ncol=M)
set.seed(34747)
iacc.exchange<-0
iacc.matrix<-array(0,c(2,M,M))
t0<-proc.time()[3]
for(iter in 1:nits){

	#Mutation Moves
	for(j in 1:M){
		new.yvec[j]<-old.yvec[j]+rnorm(1)*sigma.q
		new.vals[j]<-log(powered.dens(new.yvec[j],muvec,omega,Tv=Temps[j]))-log(Norm[j])
		if(log(runif(1))< new.vals[j]-old.vals[j]){
			old.yvec[j]<-new.yvec[j]
			old.vals[j]<-new.vals[j]
		}
	}

	#Exchange Move
	#Pick two chains at random
	jk<-sort(sample(1:M,size=2))
	iacc.matrix[1,jk[1],jk[2]]<-iacc.matrix[1,jk[1],jk[2]]+1	

	#Swap their current values
	new.yvec[jk[1]]<-old.yvec[jk[2]]
	new.yvec[jk[2]]<-old.yvec[jk[1]]

	#Work out the new density values
	new.vals[jk[1]]<-log(powered.dens(new.yvec[jk[1]],muvec,omega,Tv=Temps[jk[1]]))-log(Norm[jk[1]])	
	new.vals[jk[2]]<-log(powered.dens(new.yvec[jk[2]],muvec,omega,Tv=Temps[jk[2]]))-log(Norm[jk[2]])	

	if(log(runif(1)) < (new.vals[jk[1]]+new.vals[jk[2]])-(old.vals[jk[1]]+old.vals[jk[2]])){
		old.yvec[jk]<-new.yvec[jk]
		old.vals[jk]<-new.vals[jk]
		iacc.exchange<-iacc.exchange+1
		iacc.matrix[2,jk[1],jk[2]]<-iacc.matrix[2,jk[1],jk[2]]+1
	}
	
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		Ysamp.pop[ico,]<-old.yvec
	}

}
t1<-proc.time()[3]
print(t1-t0)

acc.rates<-iacc.matrix[2,,]/iacc.matrix[1,,]

print(paste("Overall Exchange acceptance rate is",round(iacc.exchange/nits,4)))
par(mfrow=c(M,2))
yvec<-seq(-15,15,by=0.01)

for(j in 1:M){
	hist(Ysamp.pop[,j],breaks=seq(-15,15,by=0.5),ylim=range(0,500),xlab="Y",main="Sample from Mixture Distribution")
	box()
	lines(yvec,stored.dv[,j]*nsamp*0.5,type="l",xlab="y")
	plot(Ysamp.pop[,j],ylim=range(-15,15),ylab="Y")
	title(paste("T = ",Temps[j]))
}