In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo: algoritmo de Metropolis com proposta de passeio aleatorio (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)
}
In [2]:
#==============================================================================
# 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)
In [3]:
#==============================================================================
# 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]))
In [4]:
#==============================================================================
# Algoritmo de Metropolis:
# Proposta de passeio aleatorio normal
# Consideraremos a proposta com matriz de covariancias diagonal (0,22 ; 0,71)
#==============================================================================
# Matriz de covariancias da proposta
PSI <- diag(c(0.22,0.71))
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,x[i,],PSI)
aux <- ( fx(xprop,omega,mu1,mu2,SIGMA1,SIGMA2) /
fx(x[i,],omega,mu1,mu2,SIGMA1,SIGMA2) )
alpha <- min(1,aux)
if(runif(1)<= alpha){
x[i+1,] <- xprop
ind <- ind+1
}else{
x[i+1,] <- x[i,]
}
}
In [5]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)
[1] 0,6086
In [6]:
#==============================================================================
# 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)
In [7]:
#==============================================================================
# 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))
In [8]:
#==============================================================================
# Algoritmo de Metropolis:
# Proposta de passeio aleatorio normal
# Consideraremos a proposta com matriz de covariancias diagonal (0,71 ; 0,71)
#==============================================================================
# Matriz de covariancias da proposta
PSI <- diag(c(0.71,0.71))
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,x[i,],PSI)
aux <- ( fx(xprop,omega,mu1,mu2,SIGMA1,SIGMA2) /
fx(x[i,],omega,mu1,mu2,SIGMA1,SIGMA2) )
alpha <- min(1,aux)
if(runif(1)<= alpha){
x[i+1,] <- xprop
ind <- ind+1
}else{
x[i+1,] <- x[i,]
}
}
In [9]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)
[1] 0,5426
In [10]:
#==============================================================================
# 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)
In [11]:
#==============================================================================
# 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))
In [12]:
#==============================================================================
# Algoritmo de Metropolis:
# Proposta de passeio aleatorio normal
# Consideraremos a proposta com matriz de covariancias diagonal (1,41 ; 0,71)
#==============================================================================
# Matriz de covariancias da proposta
PSI <- diag(c(1.41,0.71))
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,x[i,],PSI)
aux <- ( fx(xprop,omega,mu1,mu2,SIGMA1,SIGMA2) /
fx(x[i,],omega,mu1,mu2,SIGMA1,SIGMA2) )
alpha <- min(1,aux)
if(runif(1)<= alpha){
x[i+1,] <- xprop
ind <- ind+1
}else{
x[i+1,] <- x[i,]
}
}
In [13]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)
[1] 0,4984
In [14]:
#==============================================================================
# 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)
In [15]:
#==============================================================================
# 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))
In [16]:
#==============================================================================
# No algoritmo de Metropolis com proposta de passeio aleatorio, uma taxa de
# aceitacao alta pode ser ruim. Uma taxa de aceitacao menor pode ser melhor.
# No primeiro caso, a cadeia esta' se movendo localmente com pequenos passos
# gerando uma amostra bastante (auto)correlacionada. No segundo caso, a taxa
# de aceitacao menor pode ocorrer porque os saltos sao maiores permitindo
# explorar o suporte da distribuicao mais rapidamente.
#
#==============================================================================
In [17]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================