# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo: algoritmo de Metropolis
#==============================================================================
px  <- function(x,n,a){
	(x^(-n/2))*exp(-0.5*a/x)
}

#==============================================================================
# Avaliacao da probabilidade de aceitacao
#==============================================================================
n    <- 5
a    <- 4

# proposto    /  corrente
print(px(39.82,n,a) /  px(1,n,a))

[1] 0,0007023004

#==============================================================================
# Metropolis com proposta U(0,100)
#==============================================================================
set.seed(1234)

M    <- 5000         # Numero de iteracoes
x    <- rep(NA,M+1)
x[1] <- 1            # Valor inicial
ind  <- 0  

for(i in 1:M){
	xprop <- runif(1,0,100)
	aux   <- px(xprop,n,a)/px(x[i],n,a)
	alpha <- min(1,aux)
	if(runif(1)<= alpha){
		x[i+1] <- xprop
		ind    <- ind+1
	}else{
		x[i+1] <- x[i]
	}
}

#==============================================================================
# Taxa de aceitacao
#
# A cadeia fica muitas iteracoes em um valor de x. 
# Isso se deve a proposta ruim!
#==============================================================================
print(ind/M)

[1] 0,0464

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

#==============================================================================
# Mudando a proposta para uma gama(0.09,0.01)
# Media 9, variancia 900 e desvio padrao 30
# Escolha para se ter cuidado com a cauda da funcao objetivo 
#==============================================================================
set.seed(1234)

M    <- 5000         # Numero de iteracoes
x    <- rep(NA,M+1)
x[1] <- 1            # Valor inicial
ind  <- 0  

for(i in 1:M){
	xprop <- rgamma(1,0.09,0.01)
	numer <- px(xprop,n,a)*dgamma(x[i],0.09,0.01)
	denom <- px(x[i],n,a)*dgamma(xprop,0.09,0.01)
	alpha <- min(1,numer/denom)
	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,1892

#==============================================================================
# Algumas analises graficas
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(5,5,5),
    mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
ts.plot(x[2:M],xlab="Iteração",ylab=expression(x),xlim=c(0,5000),ylim=c(0,200))
hist(x[2:M],prob=T,main="",xlab=expression(x),ylab="Densidade",xlim=c(0,200),
    ylim=c(0,0.1))
acf(x[2:M],lag.max=40,xlab="Defasagem",main="",xlim=c(0,40),ylim=c(-0.2,1))

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