In [1]:
# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo: algoritmo de Metropolis com proposta de passeio aleatorio (normal)
#
# Objetivo: gerar uma amostra da distribuicao a posteriori de theta no 
# "problema dos animais"
#
#==============================================================================
# Funcao de verossimilhanca
#==============================================================================
L  <- function(theta,y){
	((2+theta)^y[1])*((1-theta)^y[2])*(theta^y[3])
}
In [2]:
#==============================================================================
# Valores observados
#==============================================================================
y  <- c(125,38,34)
In [3]:
#==============================================================================
# Proposta normal com desvio padrao 'sdprop'
# Utilize tambem sdprop={0.1, 1.0, 2.0} e verifique a respectiva taxa de 
# aceitacao
#==============================================================================
sdprop <- 0.5 

set.seed(1234)

M        <- 10000       # Numero de iteracoes
theta    <- rep(NA,M+1)
theta[1] <- 0.5         # Valor inicial
xi       <- 0 
ind      <- 0  

for(i in 1:M){
	xiprop    <- rnorm(1,xi,sdprop)
	thetaprop <- exp(xiprop)/(1+exp(xiprop))
	aux       <- ( L(thetaprop,y)*thetaprop*(1-thetaprop) / 
			   (L(theta[i],y)*theta[i]*(1-theta[i]))  )
	alpha     <- min(1,aux)
	if(runif(1)<= alpha){
		xi         <- xiprop
		theta[i+1] <- thetaprop
		ind        <- ind+1
	}else{
		theta[i+1] <- theta[i]
	}
}
In [4]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)

[1] 0,4604
In [5]:
#==============================================================================
# Algumas 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(theta[2:(M+1)],xlab="Iteração",ylab=expression(theta))
hist(theta[2:(M+1)],nclass=20,xlab=expression(theta),
     main="",ylab="Densidade",prob=T)
acf(theta[2:(M+1)],lag.max=50,main="",xlab="Defasagem",xlim=c(0,50),
    ylim=c(-0.2,1))
plot(density(theta,adjust=1.5,from=0.4,to=0.85),main=" ",
     lwd=3,ylab="Densidade",xlab=expression(theta))
No description has been provided for this image
In [6]:
#==============================================================================
# Media estimada
#==============================================================================
print(mean(theta[2:(M+1)]))

[1] 0,6242418
In [7]:
#==============================================================================
# Desvio padrao estimado
#==============================================================================
print(sd(theta[2:(M+1)]))

[1] 0,05056477
In [8]:
#==============================================================================
# Intervalo de credibilidade de 95% estimado
#==============================================================================
print(quantile(theta[2:(M+1)],c(.025,.975)))

     2,5%     97,5% 
0,5216856 0,7199032
In [9]:
# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================