library(mnormt)
library(rgl)
library(mclust)


x	<-	seq(-5, 5, length.out=100)
y	<-	seq(-5, 5, length.out=100)

m1	<-	c(0,1)
S1  <-	matrix(c(1, -0.7, -0.7, 2), 2, 2)
m2	<-	c(-1,-2)
S2  <-	matrix(c(1, 0, 0, 1), 2, 2)

f1	<-	outer(x, y)
for(i in 1:length(x)){
	for(j in 1:length(y)){
		f1[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m1, varcov=S1)
	}	
}

f2	<-	outer(x, y)
for(i in 1:length(x)){
	for(j in 1:length(y)){
		f2[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m2, varcov=S2)
	}	
}

p	<-	0.5
f	<-	p*f1+(1-p)*f2
contour(f, nlevels=30)
persp3d(x,y,f, col=gray(0.7), xlab="x1", ylab="x2", zlab="density")
c	<-	0.055
contour(f, levels=c)


set.seed(701618)
# generate data
n	<-	1000
xx1	<-	rmnorm(n, mean=m1, varcov=S1)
xx2	<-	rmnorm(n, mean=m2, varcov=S2)
uu	<-	rbinom(n, size=1, prob=p)
xx	<-	uu*xx1+(1-uu)*xx2

# get MLE estimates using EM algorithm
start	<-	mstep(modelName="VVV", data=xx, z=unmap(rbinom(n,1,0.5)))
em.fit	<-	em("VVV", xx, parameters=start$parameters)
phat	<-	em.fit$parameters$pro[1]
m1hat	<-	em.fit$parameters$mean[,1]
m2hat	<-	em.fit$parameters$mean[,2]
S1hat	<-	em.fit$parameters$variance$sigma[,,1]
S2hat	<-	em.fit$parameters$variance$sigma[,,2]

f1hat	<-	outer(x, y)
for(i in 1:length(x)){
	for(j in 1:length(y)){
		f1hat[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m1hat, varcov=S1hat)
	}	
}

f2hat	<-	outer(x, y)
for(i in 1:length(x)){
	for(j in 1:length(y)){
		f2hat[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m2hat, varcov=S2hat)
	}	
}

fhat	<-	phat*f1hat+(1-phat)*f2hat
contour(fhat, nlevels=30)
persp3d(x,y,fhat, col=gray(0.6))

par(pty="m", mar=c(2,2,0.5,0.5))
contour(fhat, levels=c, col="black", lwd=4, drawlabels=FALSE)
contour(f, levels=c, col=gray(0.5), lwd=4, drawlabels=FALSE, add=TRUE)




#	get 95% CI for HDR contour using bootstrap
#   this takes a long time and creates a very large file!

# bootstrap simulations
zzz		<-	file("bootstrap.out", "w")

set.seed(1048404)
B	<-	1000
for(k in 1:B){
	star		<-	sample(1:n, size=n, replace=TRUE)
	xxstar		<-	xx[star,]
	start		<-	mstep(modelName="VVV", data=xxstar, z=unmap(rbinom(n,1,0.5)))
	em.fit		<-	em("VVV", xxstar, parameters=start$parameters)
	pstar		<-	em.fit$parameters$pro[1]
	m1star		<-	em.fit$parameters$mean[,1]
	m2star		<-	em.fit$parameters$mean[,2]
	S1star		<-	em.fit$parameters$variance$sigma[,,1]
	S2star		<-	em.fit$parameters$variance$sigma[,,2]

	f1star		<-	outer(x, y)
	for(i in 1:length(x)){
		for(j in 1:length(y)){
			f1star[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m1star, varcov=S1star)
		}	
	}

	f2star	<-	outer(x, y)
	for(i in 1:length(x)){
		for(j in 1:length(y)){
			f2star[i,j]	<-	dmnorm(c(x[i],y[j]), mean=m2star, varcov=S2star)
		}	
	}

	fstar	<-	pstar*f1star+(1-pstar)*f2star
	cat(as.vector(abs(fstar-fhat)), "\n", file= zzz, sep=" ")
	print(k)	
}

close(zzz)

# read in and analyze bootstrap results
boot		<-	read.table("bootstrap.out")
sup.boot	<-	apply(boot,1, max)
q.boot		<-	quantile(sup.boot, probs=0.95)


par(pty="m", mar=c(2,2,0.5,0.5))
contour(fhat, levels=c, col="black", lwd=4, drawlabels=FALSE)
contour(fhat, levels=c(c-q.boot, c+q.boot), col="black", lwd=4, lty=2, add=TRUE, drawlabels=FALSE)
contour(f, levels=c, col=gray(0.5), lwd=4, drawlabels=FALSE, add=TRUE)