In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo: algoritmo de Gibbs
#
# Objetivo: gerar uma amostra da distribuicao a posteriori de (lambda,phi,m)
# no modelo Poisson com mudanca no ponto m
# Na descricao do modelo, m varia entre 1 e (n-1)
#
#==============================================================================
# Considere a amostra de tamanho n=30 do modelo
# Poisson com mudanca em m=20
# Para i= 1,...,20, y_i ~ Poi(lambda) e lambda=1,5
# Para i=21,...,30, y_i ~ Poi(phi) e phi=4.
#==============================================================================
# Gerando dados artificiais (a amostra aleatoria)
#==============================================================================
set.seed(12345)
n <- 30
pm <- 20
lambday <- 1.5
phiy <- 4
y <- c(rpois(pm,lambday),rpois(n-pm,phiy))
print(y) # Dados
print(length(y))
[1] 2 3 2 3 1 0 1 1 2 5 0 0 2 0 1 1 1 1 0 4 4 3 8 5 5 3 5 4 2 4 [1] 30
In [2]:
#==============================================================================
# Serie (temporal) gerada
#==============================================================================
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")
ts.plot(y,xlab="Elemento da Amostra",ylab="Valor observado",lty=3)
points(y,pch=15,lwd=2,cex=1.3)
In [3]:
#==============================================================================
# Hiperparametros (distribuicao a priori)
#==============================================================================
a <- 0.001
b <- 0.001
c <- 0.001
d <- 0.001
In [4]:
#==============================================================================
# Algoritmo de Gibbs
#==============================================================================
M <- 10000 # Numero de iteracoes
lambda <- rep(NA,M)
phi <- rep(NA,M)
m <- rep(NA,M)
lambda[1] <- 10 # Valor inicial
phi[1] <- 10 # Valor inicial
m[1] <- 10 # Valor inicial
In [5]:
#==============================================================================
# Declaracao de valores auxiliares
#==============================================================================
s1 <- cumsum(y)
s2 <- sum(y)-s1
s1 <- s1[1:(n-1)] # Diminuindo o comprimento
s2 <- s2[1:(n-1)] # Diminuindo o comprimento
ind <- seq(1,n-1,1) # Note que m esta' entre 1 e (n-1)
for (i in 2:M){
# Gerando de lambda
lambda[i] <- rgamma(1,a+s1[m[i-1]],b+m[i-1])
# Gerando de phi
phi[i] <- rgamma(1,c+s2[m[i-1]],d+n-m[i-1])
# Gerando de m
prob <- (lambda[i]^s1)*(phi[i]^s2)*exp(ind*(phi[i]-lambda[i]))
prob <- prob/sum(prob)
m[i] <- sample.int((n-1),size=1,prob=prob)
}
In [6]:
#==============================================================================
# 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(lambda[indc],xlab="Iteração",ylab=expression(lambda))
ts.plot(phi[indc],xlab="Iteração",ylab=expression(phi))
ts.plot(m[indc],xlab="Iteração",ylab=expression(m))
In [7]:
#==============================================================================
# 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(lambda[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(lambda),ylab="Densidade")
hist(phi[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(phi),ylab="Densidade")
hist(m[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(m),ylab="Densidade")
In [8]:
#==============================================================================
# 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(lambda[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",lambda)))
acf(phi[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",phi)))
acf(m[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",m)))
In [9]:
#==============================================================================
# Eliminando alguns valores iniciais (burn-in, warm-up)
#==============================================================================
indc <- 1001:M
In [10]:
#==============================================================================
# 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(lambda[indc],xlab="Iteração",ylab=expression(lambda))
ts.plot(phi[indc],xlab="Iteração",ylab=expression(phi))
ts.plot(m[indc],xlab="Iteração",ylab=expression(m))
In [11]:
#==============================================================================
# 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(lambda[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(lambda),ylab="Densidade")
hist(phi[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(phi),ylab="Densidade")
hist(m[indc],prob=T,col="darkblue",main="",nclass=50,
xlab=expression(m),ylab="Densidade")
In [12]:
#==============================================================================
# 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(lambda[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",lambda)))
acf(phi[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",phi)))
acf(m[indc],xlab="Defasagem",main="",ylim=c(-0.2,1),
ylab=expression(paste("ACF:",m)))
In [13]:
#==============================================================================
# Estimacao pontual
#==============================================================================
print("===== Media =====")
print(c(mean(lambda[indc]),mean(phi[indc]),mean(m[indc])))
print("===== Mediana =====")
print(c(median(lambda[indc]),median(phi[indc]),median(m[indc])))
print("===== Desvio Padrao =====")
print(c(sd(lambda[indc]),sd(phi[indc]),sd(m[indc])))
[1] "===== Media =====" [1] 1,410796 4,268015 19,237000 [1] "===== Mediana =====" [1] 1,389247 4,229010 19,000000 [1] "===== Desvio Padrao =====" [1] 0,2811932 0,6490724 0,8712850
In [14]:
#==============================================================================
# Estimacao intervalar
#==============================================================================
print("===== LAMBDA =====")
print(quantile(lambda[indc],c(0.025,0.975)))
print("===== PHI =====")
print(quantile(phi[indc],c(0.025,0.975)))
print("===== M =====")
print(quantile(m[indc],c(0.025,0.975)))
[1] "===== LAMBDA ====="
2,5% 97,5%
0,9168375 2,0217546
[1] "===== PHI ====="
2,5% 97,5%
3,097779 5,655631
[1] "===== M ====="
2,5% 97,5%
18 22
In [15]:
#==============================================================================
# Matriz de correlacao
#==============================================================================
print(cor(cbind(lambda[indc],phi[indc],m[indc])))
[,1] [,2] [,3] [1,] 1,0000000 -0,0106998 0,2138541 [2,] -0,0106998 1,0000000 0,1418969 [3,] 0,2138541 0,1418969 1,0000000
In [16]:
#==============================================================================
# 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(lambda[indc],phi[indc],m[indc]),pch=15,
labels=c(expression(lambda),expression(phi),expression(m)))
In [17]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================