# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo: algoritmo de Metropolis com proposta independente (normal)
#
# Objetivo: gerar uma amostra de uma mistura de normais bivariada.
#
#==============================================================================
library("mnormt")
fx  <- function(x,omega,mu1,mu2,SIGMA1,SIGMA2){
	saida <- (    omega*dmnorm(x,mu1,SIGMA1)
               +(1-omega)*dmnorm(x,mu2,SIGMA2) )
	return(saida)
}

#==============================================================================
# Neste exemplo, temos:
#==============================================================================
omega  <- 0.7
mu1    <- c(4,5)
mu2    <- c(0.7,3.5)
SIGMA1 <- matrix(c(1,0.7,0.7,1),2,2,T)
SIGMA2 <- matrix(c(1,-0.7,-0.7,1),2,2,T)

x1     <- seq(-2,7,length=100)
x2     <- seq(0,8,length=100)
z      <- matrix(NA,100,100)
for(i in 1:100)
	for(j in 1:100)
		z[i,j] <- fx(c(x1[i],x2[j]),omega,mu1,mu2,
                         SIGMA1,SIGMA2)

#==============================================================================
# Curvas de nivel
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
    mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
contour(x1,x2,z,lwd=2,col="blue",xlab=expression(x[1]),
        ylab=expression(x[2]))

#==============================================================================
# Algoritmo de Metropolis: 
# Proposta independente normal com media (3 , 5) e matriz de covariancias 2I_2
#
#==============================================================================
MU  <- c(3,5)        # vetor de medias da proposta
PSI <- diag(c(2,2))  # matriz de covariancias da proposta

set.seed(1234)

M     <- 5000              # Numero de iteracoes
x     <- matrix(NA,M+1,2)
x[1,] <- c(6,8)            # Valor inicial
ind   <- 0  

for(i in 1:M){
	xprop <- rmnorm(1,MU,PSI)
	aux1  <- ( fx(xprop,omega,mu1,mu2,SIGMA1,SIGMA2) / 
		       fx(x[i,],omega,mu1,mu2,SIGMA1,SIGMA2) )
	aux2  <- dmnorm(xprop,MU,PSI)/dmnorm(x[i,],MU,PSI)
	alpha <- min(1,(aux1/aux2))
	if(runif(1)<= alpha){
		x[i+1,] <- xprop
		ind     <- ind+1
	}else{
		x[i+1,] <- x[i,]
	}
}

#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)

[1] 0,379

#==============================================================================
# Dispersao 
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
    mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
contour(x1,x2,z,lwd=2,col="blue",xlab=expression(x[1]),
        ylab=expression(x[2]),xlim=c(-3,8),ylim=c(0,10))
lines(x[,1],x[,2],col=1,lwd=2)

#==============================================================================
# Analises graficas
#==============================================================================
par(mfrow=c(2,2),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
    mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
ts.plot(x[2:M,1],xlab="Iteração",ylab=expression(x[1]))
ts.plot(x[2:M,2],xlab="Iteração",ylab=expression(x[2]))
acf(x[2:M,1],lag.max=50,main="",xlab="Defasagem",xlim=c(0,50),ylim=c(-0.2,1))
acf(x[2:M,2],lag.max=50,main="",xlab="Defasagem",xlim=c(0,50),ylim=c(-0.2,1))

# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================