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 beta_1 e beta_2
# de um modelo de regressao logistica
#
#==============================================================================
# Funcao de verossimilhanca
#==============================================================================
L <- function(beta,y,x){
aux <- exp(beta[1]+beta[2]*x)
prob <- aux/(1+aux)
saida <- prod((prob^y)*((1-prob)^(1-y)))
return(saida)
}
In [2]:
#==============================================================================
# Distribuicao a priori
#==============================================================================
priori <- function(beta){
dnorm(beta[1],0,10)*dnorm(beta[2],0,10)
}
In [3]:
#==============================================================================
# Valores observados
#==============================================================================
y <- c(0,1,0,0,0,0,0,0,1,1,1,0,0,1,0,0,0,0,0,0,1,0,1)
x <- c(66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58)
In [4]:
#==============================================================================
# Proposta normal
#==============================================================================
sdprop1 <- 0.5
sdprop2 <- 0.5
set.seed(1234)
M <- 10000 # Numero de iteracoes
beta <- matrix(NA,M+1,2)
beta[1,] <- c(0,0) # Valor inicial
ind <- 0
for(i in 1:M){
# Proposta de passeio aleatorio
betaprop <- c(rnorm(1,beta[i,1],sdprop1),
rnorm(1,beta[i,2],sdprop2))
numer <- L(betaprop,y,x)*priori(betaprop)
denom <- L(beta[i,],y,x)*priori(beta[i,])
alpha <- min(1,(numer/denom))
if(runif(1)<= alpha){
beta[i+1,] <- betaprop
ind <- ind+1
}else{
beta[i+1,] <- beta[i,]
}
}
In [5]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)
[1] 0,0189
In [6]:
#==============================================================================
# Trajetoria da cadeia
#==============================================================================
beta1 <- seq( 0, 25, length=100)
beta2 <- seq(-0.4, 0, length=100)
z <- matrix(NA,100,100)
for(i in 1:100)
for(j in 1:100)
z[i,j] <- ( L(c(beta1[i],beta2[j]),y,x) *
priori(c(beta1[i],beta2[j])) )
par(mfrow=c(2,2),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(beta[2:(M+1),1],xlab="Iteração",
ylab=expression(beta[1]))
ts.plot(beta[2:(M+1),2],xlab="Iteração",
ylab=expression(beta[2]))
plot(beta[2:(M+1),],type="l",xlab=expression(beta[1]),
ylab=expression(beta[2]))
contour(beta1,beta2,z,col=4,xlab=expression(beta[1]),
ylab=expression(beta[2]))
In [7]:
#==============================================================================
# Mudando a variancia (desvio padrao) da proposta
#==============================================================================
sdprop1 <- 2
sdprop2 <- 0.05
set.seed(1234)
M <- 100000 # Numero de iteracoes
beta <- matrix(NA,M+1,2)
beta[1,] <- c(0,0) # Valor inicial
ind <- 0
for(i in 1:M){
# Proposta de passeio aleatorio
betaprop <- c(rnorm(1,beta[i,1],sdprop1),
rnorm(1,beta[i,2],sdprop2))
numer <- L(betaprop,y,x)*priori(betaprop)
denom <- L(beta[i,],y,x)*priori(beta[i,])
alpha <- min(1,(numer/denom))
if(runif(1)<= alpha){
beta[i+1,] <- betaprop
ind <- ind+1
}else{
beta[i+1,] <- beta[i,]
}
}
In [8]:
#==============================================================================
# Taxa de aceitacao
#==============================================================================
print(ind/M)
[1] 0,16109
In [9]:
#==============================================================================
# Analise grafica
#==============================================================================
par(mfrow=c(2,3),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(beta[2:(M+1),1],xlab="Iteração",
ylab=expression(beta[1]))
ts.plot(beta[2:(M+1),2],xlab="Iteração",
ylab=expression(beta[2]))
plot(beta[2:(M+1),],type="l",xlab=expression(beta[1]),
ylab=expression(beta[2]))
contour(beta1,beta2,z,add=T,lwd=2,col=4,
xlab=expression(beta[1]),ylab=expression(beta[2]))
acf(beta[2:(M+1),1],lag.max=60,xlab="Defasagem",main="",
ylab=expression(paste("ACF:",beta[1])))
acf(beta[2:(M+1),2],lag.max=60,xlab="Defasagem",main="",
ylab=expression(paste("ACF:",beta[2])))
In [10]:
#==============================================================================
# Trajetoria da cadeia
#==============================================================================
par(mfrow=c(2,2),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(beta[2:(M+1),1],xlab="Iteração",
ylab=expression(beta[1]))
ts.plot(beta[2:(M+1),2],xlab="Iteração",
ylab=expression(beta[2]))
plot(beta[2:(M+1),],pch=".",xlab=expression(beta[1]),
ylab=expression(beta[2]))
contour(beta1,beta2,z,add=T,lwd=2,col=4,
xlab=expression(beta[1]),ylab=expression(beta[2]))
In [11]:
#==============================================================================
# Estimacao do beta_1
#==============================================================================
print("===== Media =====")
print(mean(beta[2:(M+1),1]))
print("===== Desvio Padrao =====")
print(sd(beta[2:(M+1),1]))
print("===== Intervalo de Credibilidade =====")
print(quantile(beta[2:(M+1),1],c(.025,.975)))
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")
hist(beta[1002:(M+1),1],xlab=expression(beta[1]),
ylab="Densidade",prob=T,col="darkblue",main="")
[1] "===== Media ====="
[1] 11,8422
[1] "===== Desvio Padrao ====="
[1] 5,070316
[1] "===== Intervalo de Credibilidade ====="
2,5% 97,5%
2,775938 22,751357
In [12]:
#==============================================================================
# Estimacao do beta_2
# Regressao significativa?
#==============================================================================
print("===== Media =====")
print(mean(beta[2:(M+1),2]))
print("===== Desvio Padrao =====")
print(sd(beta[2:(M+1),2]))
print("===== Intervalo de Credibilidade =====")
print(quantile(beta[2:(M+1),1],c(.025,.975)))
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")
hist(beta[1002:(M+1),2],xlab=expression(beta[2]),
ylab="Densidade",prob=T,col="darkblue",main="")
[1] "===== Media ====="
[1] -0,1863503
[1] "===== Desvio Padrao ====="
[1] 0,07444003
[1] "===== Intervalo de Credibilidade ====="
2,5% 97,5%
2,775938 22,751357
In [13]:
#==============================================================================
# Inferencia para pi_i (prob_i) ?
#==============================================================================
prob <- matrix(NA,M,length(x))
for(i in 1:length(x)){
aux <- beta[2:(M+1),1]+beta[2:(M+1),2]*x[i]
prob[,i] <- exp(aux)/(1+exp(aux))
}
par(mfrow=c(2,2),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(5,5,5),
mar=c(4.5,5,3,1),cex.main=2.0,bty="n")
hist(prob[,14],xlab=expression(pi[14]),
ylab="Densidade",prob=T,col="darkblue",
main=expression(paste(x[14],"=53",sep="")))
hist(prob[,15],xlab=expression(pi[15]),
ylab="Densidade",prob=T,col="darkblue",
main=expression(paste(x[15],"=67",sep="")))
hist(prob[,18],xlab=expression(pi[18]),
ylab="Densidade",prob=T,col="darkblue",
main=expression(paste(x[18],"=81",sep="")))
In [14]:
#==============================================================================
# Qual e' o efeito da covariavel ?
#==============================================================================
# Diminuir a temperatura em 10 graus aumenta a chance de falha em OR.10
# OR: Odds ratio (failure)= (1-pi)/pi
#==============================================================================
OR.10 <- exp(-beta[2:(M+1),2]*10)
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")
hist(OR.10,xlab="OR.10",ylab="Densidade",prob=T,ylim=c(0,0.08),
col="darkblue", main="")
print(mean(OR.10))
[1] 8,733546
In [15]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================