# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo: algoritmo de Gibbs
#
# Objetivo: gerar uma amostra da distribuicao a posteriori de (mu,phi) no 
# modelo normal. O parametro phi e' a precisao (inverso da variancia)
#
#==============================================================================
# Considere a amostra de tamanho n=20 do modelo normal gerado com valor 
# mu=5 e sigma^2=2 (phi=1/2).
#==============================================================================
y <- c(4.05,5.13,6.93,4.85,5.72,5.50,5.21,5.17,9.23,3.32,
       4.33,4.20,5.55,7.96,3.49,4.71,6.12,4.10,2.82,5.85)

#==============================================================================
# Dados
#==============================================================================
n      <- length(y)
ybarra <- mean(y)

#==============================================================================
# Priori
#==============================================================================
a0     <- 0.01
b0     <- 0.01
tau02  <- 100

#==============================================================================
# Algoritmo de Gibbs
#==============================================================================
M      <- 10000      # Numero de iteracoes
mu     <- rep(NA,M)
phi    <- rep(NA,M)
mu[1]  <-  0         # Valor inicial
phi[1] <-  1         # Valor inicial

for (i in 2:M){
	# Gerando de phi
	somaq  <- sum( (y - mu[i-1])^2 )
	phi[i] <- rgamma(1,(n+a0)/2,(somaq+b0)/2)

	# Gerando de mu
	delta2 <- 1/(n*phi[i]+1/tau02)
	xi     <- n*ybarra*phi[i]*delta2
	mu[i]  <- rnorm(1,xi,sqrt(delta2))
}

sigma2 <- 1/phi

#==============================================================================
# Analise grafica - Parte 1
#==============================================================================
indc <- 1:M

par(mfrow=c(3,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(mu[indc],xlab="Iteração",ylab=expression(mu))
ts.plot(phi[indc],xlab="Iteração",ylab=expression(phi))
ts.plot(sigma2[indc],xlab="Iteração",
        ylab=expression(sigma^2))

#==============================================================================
# Analise grafica - Parte 2
#==============================================================================
par(mfrow=c(3,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")
hist(mu[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(mu),ylab="Densidade")
hist(phi[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(phi),ylab="Densidade")
hist(sigma2[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(sigma^2),ylab="Densidade")

#==============================================================================
# Analise grafica - Parte 3
#==============================================================================
par(mfrow=c(3,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")
acf(mu[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",mu)))
acf(phi[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",phi)))
acf(sigma2[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",sigma^2)))

#==============================================================================
# Eliminando alguns valores iniciais (burn-in, warm-up)
#==============================================================================
indc <- 1001:M

#==============================================================================
# Analise grafica - Parte 1
#==============================================================================
par(mfrow=c(3,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(mu[indc],xlab="Iteração",ylab=expression(mu))
ts.plot(phi[indc],xlab="Iteração",ylab=expression(phi))
ts.plot(sigma2[indc],xlab="Iteração",
        ylab=expression(sigma^2))

#==============================================================================
# Analise grafica - Parte 2
#==============================================================================
par(mfrow=c(3,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")
hist(mu[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(mu),ylab="Densidade")
hist(phi[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(phi),ylab="Densidade")
hist(sigma2[indc],prob=T,col="darkblue",main="",nclass=50,
     xlab=expression(sigma^2),ylab="Densidade")

#==============================================================================
# Analise grafica - Parte 3
#==============================================================================
par(mfrow=c(3,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")
acf(mu[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",mu)))
acf(phi[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",phi)))
acf(sigma2[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
    ylab=expression(paste("ACF:",sigma^2)))

#==============================================================================
# Estimacao pontual
#==============================================================================
print("===== Media =====")
print(c(mean(mu[indc]),mean(phi[indc]),mean(sigma2[indc])))
print("===== Mediana =====")
print(c(median(mu[indc]),median(phi[indc]),median(sigma2[indc])))
print("===== Desvio Padrao =====")
print(c(sd(mu[indc]),sd(phi[indc]),sd(sigma2[indc])))

[1] "===== Media ====="
[1] 5,1978275 0,4209187 2,6599232
[1] "===== Mediana ====="
[1] 5,1980445 0,4071197 2,4562802
[1] "===== Desvio Padrao ====="
[1] 0,3600814 0,1381719 0,9784594

#==============================================================================
# Estimacao intervalar
#==============================================================================
print(quantile(mu[indc],c(0.025,0.975)))
print(quantile(sigma2[indc],c(0.025,0.975)))

    2,5%    97,5% 
4,480989 5,895333 
    2,5%    97,5% 
1,366496 5,103680 

#==============================================================================
# Matriz de correlacao
#==============================================================================
print(cor(mu[indc],sigma2[indc]))

[1] -0,03196427

#==============================================================================
# Grafico de dispersao
#==============================================================================
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)
pairs(cbind(mu[indc],sigma2[indc]),pch=15,
     labels=c(expression(mu),expression(phi)))

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