In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo: algoritmo de Gibbs
#
# Objetivo: gerar uma amostra da distribuicao normal trivariada
#
#==============================================================================
MU <- c(-5,0,5)
SIGMA <- matrix(c(1,0.9,0.95,0.9,1,0.97,0.95,0.97,1),3,3,T)
print(SIGMA)
indc <- c(2,3)
sig21 <- 1-SIGMA[1,indc]%*%solve(SIGMA[indc,indc])%*%SIGMA[indc,1]
sd1 <- sqrt(sig21)
A1 <- SIGMA[1,indc]%*%solve(SIGMA[indc,indc])
indc <- c(1,3)
sig22 <- 1-SIGMA[2,indc]%*%solve(SIGMA[indc,indc])%*%SIGMA[indc,2]
sd2 <- sqrt(sig22)
A2 <- SIGMA[2,indc]%*%solve(SIGMA[indc,indc])
indc <- c(1,2)
sig23 <- 1-SIGMA[3,indc]%*%solve(SIGMA[indc,indc])%*%SIGMA[indc,3]
sd3 <- sqrt(sig23)
A3 <- SIGMA[3,indc]%*%solve(SIGMA[indc,indc])
[,1] [,2] [,3] [1,] 1,00 0,90 0,95 [2,] 0,90 1,00 0,97 [3,] 0,95 0,97 1,00
In [2]:
#==============================================================================
# Inicio
#==============================================================================
set.seed(1234)
M <- 10000 # Numero de iteracoes
x <- matrix(NA,M,3)
x[1,] <- c(0,0,0) # Valor inicial
for(i in 2:M){
# Gerando da primeira componente: x_1
indc <- c(2,3)
xi <- x[i-1,indc]
media <- MU[1]+A1%*%(xi-MU[indc])
x[i,1] <- rnorm(1,media,sd1)
# Gerando da segunda componente: x_2
indc <- c(1,3)
xi <- c(x[i,1],x[i-1,3])
media <- MU[2]+A2%*%(xi-MU[indc])
x[i,2] <- rnorm(1,media,sd2)
# Gerando da terceira componente: x_3
indc <- c(1,2)
xi <- x[i,1:2]
media <- MU[3]+A3%*%(xi-MU[indc])
x[i,3] <- rnorm(1,media,sd3)
}
In [3]:
#==============================================================================
# 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(x[indc,1],xlab="Iteração",ylab=expression(x[1]))
ts.plot(x[indc,2],xlab="Iteração",ylab=expression(x[2]))
ts.plot(x[indc,3],xlab="Iteração",ylab=expression(x[3]))
In [4]:
#==============================================================================
# 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(x[indc,1],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[1]),ylab="Densidade")
hist(x[indc,2],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[2]),ylab="Densidade")
hist(x[indc,3],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[3]),ylab="Densidade")
In [5]:
#==============================================================================
# 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(x[indc,1],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[1])))
acf(x[indc,2],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[2])))
acf(x[indc,3],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[3])))
In [6]:
#==============================================================================
# Eliminando alguns valores iniciais (burn-in, warm-up)
#==============================================================================
indc <- 1001:M
In [7]:
#==============================================================================
# Analise grafica - Parte 1
#==============================================================================
par(mfrow=c(3,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
ts.plot(x[indc,1],xlab="Iteração",ylab=expression(x[1]))
ts.plot(x[indc,2],xlab="Iteração",ylab=expression(x[2]))
ts.plot(x[indc,3],xlab="Iteração",ylab=expression(x[3]))
In [8]:
#==============================================================================
# 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(x[indc,1],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[1]),ylab="Densidade")
hist(x[indc,2],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[2]),ylab="Densidade")
hist(x[indc,3],prob=T,nclas=50,col="darkblue",main="",
xlab=expression(x[3]),ylab="Densidade")
In [9]:
#==============================================================================
# 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(x[indc,1],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[1])))
acf(x[indc,2],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[2])))
acf(x[indc,3],main="",xlab="Defasagem",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",x[3])))
In [10]:
#==============================================================================
# Medias estimada e verdadeira
#==============================================================================
print(apply(x[indc,],2,mean))
print(MU)
[1] -5,01262899 -0,01447908 4,98543634 [1] -5 0 5
In [11]:
#==============================================================================
# Matrizes de covariancia estimada e verdadeira
#==============================================================================
print(cov(x[indc,]))
print(cor(x[indc,]))
print(SIGMA)
[,1] [,2] [,3]
[1,] 0,9399984 0,8289977 0,8830540
[2,] 0,8289977 0,9201394 0,8918550
[3,] 0,8830540 0,8918550 0,9242828
[,1] [,2] [,3]
[1,] 1,0000000 0,8913802 0,9473736
[2,] 0,8913802 1,0000000 0,9670858
[3,] 0,9473736 0,9670858 1,0000000
[,1] [,2] [,3]
[1,] 1,00 0,90 0,95
[2,] 0,90 1,00 0,97
[3,] 0,95 0,97 1,00
In [12]:
#==============================================================================
# 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(x[indc,],pch=".",
labels=c(expression(x[1]),expression(x[2]),expression(x[3])))
In [13]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================