#install.packages("mvtnorm")
#library("mvtnorm")

post = function(th1,th2){
  pi1*dmvt(c(th1,th2),delta=m1,sigma=S1,df=df1,log=FALSE)+
  pi2*dmvt(c(th1,th2),delta=m2,sigma=S2,df=df2,log=FALSE)+
  (1-pi1-pi2)*dmvt(c(th1,th2),delta=m3,sigma=S3,df=df3,log=FALSE)
}

m1 = c(-1,1)
S1 = matrix(c(1,-.8,-.8,1),2,2)
m2 = c(2,-2)
S2 = matrix(c(1,-.5,-.5,1),2,2)
m3 = c(1,-2)
S3 = matrix(c(1,.5,.5,1),2,2)
pi1 = 0.3
pi2 = 0.2
df1 = 3
df2 = 3
df3 = 3


pdf(file="3comp-mixture.pdf",width=12,height=10)
# Contour plots
# -------------
N         = 200
theta1    = seq(-10,10,length=N)
theta2    = seq(-10,10,length=N)
h         = theta1[2]-theta1[1]
posterior = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    posterior[i,j] = post(theta1[i],theta2[j])

par(mfrow=c(1,2))
contour(theta1,theta2,posterior,nlevels=20,drawlabels=FALSE,
        xlab=expression(theta[1]),ylab=expression(theta[2]))
plot(theta1,apply(h*posterior,1,sum),type="l",xlab=expression(theta),ylab="",ylim=c(0,0.3))
lines(theta2,apply(h*posterior,2,sum),col=2)
legend("topleft",legend=c(expression(theta[1]),expression(theta[2])),col=1:2,lty=1,bty="n")


# SIR in action - bivariate uniform proposal
set.seed(12135)
M = 10000
N = 10000
theta1.draw = runif(M,-10,10)
theta2.draw = runif(M,-10,10)
w = rep(0,M)
for (i in 1:M)
  w[i] = post(theta1.draw[i],theta2.draw[i])
w = w/sum(w)  
k = sample(1:M,size=N,replace=TRUE,prob=w)
theta1.post = theta1.draw[k]
theta2.post = theta2.draw[k]


par(mfrow=c(2,2))
plot(theta1.draw,theta2.draw,xlim=range(theta1.draw),ylim=range(theta2.draw),cex=0.5,pch=16,
xlab=expression(theta[1]),ylab=expression(theta[2]),main="Draws from proposal")
plot(theta1.post,theta2.post,xlim=range(theta1.draw),ylim=range(theta2.draw),cex=0.5,pch=16,
xlab=expression(theta[1]),ylab=expression(theta[2]),main="Draws from target")
contour(theta1,theta2,posterior,drawlabels=FALSE,add=TRUE,col=2,lwd=2)


hist(theta1.post,prob=TRUE,breaks=25,ylab="Posterior density",xlab="",main=expression(theta[1]))
lines(theta1,apply(h*posterior,1,sum),col=2,lwd=2)
hist(theta2.post,prob=TRUE,breaks=25,ylab="Posterior density",xlab="",main=expression(theta[2]))
lines(theta2,apply(h*posterior,2,sum),col=2,lwd=2)

thetas1 = cbind(theta1.post,theta2.post)


# SIR in action - bivariate normal proposal
proposal  = function(theta1,theta2){
  dnorm(theta1,0,3)*dnorm(theta2,0,3)
}
set.seed(12135)
M = 10000
N = 10000
theta1.draw = rnorm(M,0,3)
theta2.draw = rnorm(M,0,3)
w = rep(0,M)
for (i in 1:M)
  w[i] = post(theta1.draw[i],theta2.draw[i])/proposal(theta1.draw[i],theta2.draw[i])
w = w/sum(w)  
k = sample(1:M,size=N,replace=TRUE,prob=w)
theta1.post = theta1.draw[k]
theta2.post = theta2.draw[k]


par(mfrow=c(2,2))
plot(theta1.draw,theta2.draw,xlim=range(theta1.draw),ylim=range(theta2.draw),cex=0.5,pch=16,
xlab=expression(theta[1]),ylab=expression(theta[2]),main="Draws from proposal")
plot(theta1.post,theta2.post,xlim=range(theta1.draw),ylim=range(theta2.draw),cex=0.5,pch=16,
xlab=expression(theta[1]),ylab=expression(theta[2]),main="Draws from target")
contour(theta1,theta2,posterior,drawlabels=FALSE,add=TRUE,col=2,lwd=2)


hist(theta1.post,prob=TRUE,breaks=30,ylab="Posterior density",xlab="",main=expression(theta[1]))
lines(theta1,apply(h*posterior,1,sum),col=2,lwd=2)
hist(theta2.post,prob=TRUE,breaks=30,ylab="Posterior density",xlab="",main=expression(theta[2]))
lines(theta2,apply(h*posterior,2,sum),col=2,lwd=2)

thetas2 = cbind(theta1.post,theta2.post)


par(mfrow=c(1,2))
plot(density(thetas1[,1]),lwd=2,xlab="",ylab="Posterior density",main=expression(theta[1]))
lines(density(thetas2[,1]),col=2,lwd=2)
lines(theta1,apply(h*posterior,1,sum),col=3,lwd=2)
legend("topleft",legend=c("Exact","Uniform proposal","Gaussian proposal"),col=c(3,1,2),bty="n",lty=1,lwd=3)
plot(density(thetas1[,2]),lwd=2,xlab="",ylab="Posterior density",main=expression(theta[2]))
lines(density(thetas2[,2]),col=2,lwd=2)
lines(theta1,apply(h*posterior,2,sum),col=3,lwd=2)

dev.off()