Avaliando algumas distribuições a priori e a posteriori (Aula Prática 2)

1 Modelo normal

Considere o seguinte modelo \[\begin{eqnarray} Y_i|\theta &\stackrel{iid}\sim& N(\theta, 1), \quad i=1, \ldots, n\\ \theta &\sim& N(0,1) \end{eqnarray}\]

1 - Mostre que a distribuição a posteriori de \(\theta\) é \[\begin{eqnarray} \theta|Y_1, \ldots, Y_n &\sim& N\left(\frac{\sum_{i=1}^n{Y_i}}{n+1}, \frac{1}{n+1}\right) \end{eqnarray}\]

2 - Simule um conjunto de dados desse modelo e obtenha uma amostra da distribuição a posteriori de \(\theta\).

##########################################################
# Supondo uma amostra de tamanho n da v.a. resposta
# Y_i ~iid N(theta,1), i=1,...n
# theta ~ N(0,1)
# theta|y_1, ..., y_n ~ N( sum(y)/(n+1) , 1/(n+1) )
##########################################################

#Simulando um conjunto de dados
thetav      = 10
n           = 2000
y           = rnorm(n, thetav, sqrt(1))

#Inferindo sobre theta (obtendo uma amostra da distribuição a posteriori)
npost       = 1000
thetapost   = rnorm(npost, sum(y)/(n+1), sqrt(1/(n+1)))

3 - Calcule a média amostral e compare esse valor com a média da distribuição a posteriori. O valor verdadeiro usado para \(\theta\) está próximo da média amostral?

mean(thetapost) #média amostral
## [1] 10.04998
sum(y)/(n+1)    #média teórica
## [1] 10.04997
c(mean(thetapost), sum(y)/(n+1), thetav)
## [1] 10.04998 10.04997 10.00000

4 - Calcule a variância amostral e compare esse valor com a variância da distribuição a posteriori.

var(thetapost) #média amostral
## [1] 0.0004818731
1/(n+1)    #média teórica
## [1] 0.0004997501
c(var(thetapost), 1/(n+1))
## [1] 0.0004818731 0.0004997501

5 - Crie um intervalo usando como limite inferior o quantil amostral de 2,5% e como limite superior o quantil amostral de 97,5%. O valor verdadeiro usado para \(\theta\) pertence ao intervalo criado acima? E se ao invés de usar a amostra, você usar a distribuição a posteriori obtida em (1), o intervalo difere?

#Usando a amostra
q1a = quantile(thetapost, 0.025)
q2a = quantile(thetapost, 0.975)
c(q1a, q2a)
##     2.5%    97.5% 
## 10.00586 10.09083
if( (thetav<=q2a) & (thetav>=q1a) ){print("theta dentro do IC 95%")}else{
  print("theta FORA do IC 95%")}
## [1] "theta FORA do IC 95%"
#Usando a distribuição teórica
q1=qnorm(0.025, sum(y)/(n+1), sqrt(1/(n+1)))
q2=qnorm(0.975, sum(y)/(n+1), sqrt(1/(n+1)))
c(q1, q2)
## [1] 10.00616 10.09379
if( (thetav<=q2) & (thetav>=q1) ){print("theta dentro do IC 95%")}else{
  print("theta FORA do IC 95%")}
## [1] "theta FORA do IC 95%"
rbind(c(q1a, q2a), c(q1, q2))
##          2.5%    97.5%
## [1,] 10.00586 10.09083
## [2,] 10.00616 10.09379

6 - Compare a média a priori com a média a posteriori de \(\theta\).

0
## [1] 0
mean(thetapost)
## [1] 10.04998

7 - Compare a variância a priori com a variância a posteriori de \(\theta\).

1
## [1] 1
var(thetapost)
## [1] 0.0004818731

8 - Faça um histograma da amostra da distribuição a posteriori e represente o valor verdadeiro de \(\theta\), a média amostral, o IC amostral criado anteriormente e a curva da distrbuição a posteriori de \(\theta\).

#visualizando a amostra da distribuição a posteriori
hist(thetapost,freq=F,main="",xlab=expression(theta), 
     ylab="densidade")
abline(v=thetav,lwd=3,lty=1,col="blue")
abline(v=mean(thetapost),lwd=3,lty=1,col="red")
abline(v=quantile(thetapost,0.025),lwd=3,lty=3,col="red")
abline(v=quantile(thetapost,0.975),lwd=3,lty=3,col="red")
curve(dnorm(u, sum(y)/(n+1), sqrt(1/(n+1))),min(thetapost),max(thetapost),add=T, 
      xname = "u", lwd=2, col="purple")

9 - Faça um gráfico da distribuição a priori e da posteriori de \(\theta\).

u           = seq(-3,thetav+3,l=2000)
dpriori     = dnorm(u,0,sqrt(1))
dposteriori = dnorm(u, sum(y)/(n+1), sqrt(1/(n+1)))
minimo      = min(dpriori, dposteriori)
maximo      = max(dpriori, dposteriori)
plot(u, dpriori, lwd=2, type="l", xlab=expression(theta), ylab="densidade", bty="n", ylim=c(minimo, maximo))
lines(u, dposteriori, lwd=2, col="red")

10 - Faça um gráfico da função de log-verossimilhança.

#Analisando o log da funcao de verossimilhanca
u = seq(5,15,l=1000)
lvero = rep(0, length(u))
for(j in 1:length(u))
{
  for(i in 1:n)
  {
    lvero[j] = lvero[j] + dnorm( y[i], u[j], sqrt(1) , log=TRUE)
  }
}
plot(u, lvero, main="Função de log-verossimilhança", type="l", lwd=2,
     xlab=expression(theta), ylab="L(theta;y)")
abline(v=thetav, lty=3, lwd=2)

11 - Calcule a probabilidade a posteriori de \(\theta\) assumir um valor menor que 10.

#Calculando de forma teórica
pnorm(10, sum(y)/(n+1), sqrt(1/(n+1)))
## [1] 0.01269661
#Calculando usando a amostra
length( which(thetapost<10) )/length(thetapost)
## [1] 0.018

1.1 Distribuição normal truncada

Se a variável aleatória \(W \sim N(\mu, \sigma^2)I(a< Y < b)\), ou seja, \(W\) tem distribuição normal truncada, então a função de densidade de probabilidade de \(W\) é \[\begin{eqnarray} f(W=w|\mu, \sigma) = \frac{\frac{1}{\sigma}\phi\left(\frac{w-\mu}{\sigma}\right)}{\Phi\left( \frac{b-\mu}{\sigma}\right) - \Phi\left( \frac{a-\mu}{\sigma}\right)}, \end{eqnarray}\] sendo \(\phi()\) a função de densidade da normal padrão e \(\Phi()\) a função de distribuição acumulada da normal padrão.

E o valor esperado de \(W\) é \[\begin{eqnarray} E[W] = \mu + \frac{\phi\left( \frac{a-\mu}{\sigma}\right) - \phi\left( \frac{b-\mu}{\sigma}\right)}{\Phi\left( \frac{b-\mu}{\sigma}\right) - \Phi\left( \frac{a-\mu}{\sigma}\right)} \sigma \end{eqnarray}\]

1.2 Mudando a distribuição a priori

Considere o seguinte modelo \[\begin{eqnarray} Y_i|\theta &\stackrel{iid}\sim& N(\theta, 1), \quad i=1, \ldots, n\\ \theta &\sim& Unif(-15,15) \end{eqnarray}\]

1 - Mostre que a distribuição a posteriori de \(\theta\) é \[\begin{eqnarray} \theta|Y_1, \ldots, Y_n &\sim& N\left(\frac{\sum_{i=1}^n{Y_i}}{n+1}, \frac{1}{n+1}\right) I(-15 < \theta < 15). \end{eqnarray}\]

2 - Obtenha uma amostra da distribuição a posteriori de \(\theta\) e faça um histograma disso. Represente nesse histograma o valor verdadeiro de \(\theta\), a média amostral, um intervalo usando os quantis amostrais de 2,5% e de 97,5% e a curva da distribuição a posteriori.

#Inferindo sobre theta
thetapost2  = rep(0,npost)
mp          = sum(y)/(n+1)
vp          = 1/(n+1)
sp          = sqrt(vp)
li          = -15
ls          = 15 
#for(i in 1:npost)
#{
#   print(i)
#   while( (thetapost2[i]>ls) | (thetapost2[i]<li) ){
#     thetapost2[i]     = rnorm(1, mp, sp); print(thetapost2[i])}
#}
library(truncnorm)
thetapost2 = rtruncnorm(npost, a=li, b=ls, mean=mp, sd=sp)

hist(thetapost2,freq=F,main="",xlab=expression(theta))
abline(v=thetav,lwd=3,lty=1,col="blue")
abline(v=mean(thetapost2),lwd=3,lty=1,col="red")
abline(v=quantile(thetapost2,0.025),lwd=3,lty=3,col="red")
abline(v=quantile(thetapost2,0.975),lwd=3,lty=3,col="red")
curve(dtruncnorm(x, a=li, b=ls, mean=mp, sd=sp), add=T, col="purple", lwd=2)

#u = sort(thetapost2)
#b = pnorm((ls-mp)/sp)
#a = pnorm((li-mp)/sp)
#d = dnorm(u,mp,sp)/ (b-a)
#lines(u, d,col="purple",lwd=2)

2 Outro modelo

Considere que \[\begin{eqnarray} Y_i|\theta &\stackrel{iid}\sim& N(0, V), \quad i=1, \ldots, n\\ \tau = \frac{1}{V} &\sim& G(\alpha,\beta) \end{eqnarray}\]

Calcule a distribuição a posteriori de \(\tau\) e de \(V\) e obtenha uma amostra a posteriori de \(\tau\) e de \(V\). Represente isso graficamente.

##########################
# Outro modelo
# Y_i ~iid N(0,V), i=1,...n
# tau=V^(-1) ~ G(alfa,beta)
# tau|y ~ Gama
##########################

#Simulando um conjunto de dados
V            = 2
tau          = 1/V
n                = 1000
y              = rnorm(n, 0, sqrt(V))
alfa         = 2#0.25
beta         = 1#0.5
#Inferindo sobre tau
taupost      = rgamma(1000, alfa + n/2, beta + sum(y^2)/2)
hist(taupost,freq=F,main="",xlab=expression(tau))
abline(v=tau,lwd=3,lty=1,col="blue")
abline(v=mean(taupost),lwd=3,lty=1,col="red")
abline(v=quantile(taupost,0.025),lwd=3,lty=3,col="red")
abline(v=quantile(taupost,0.975),lwd=3,lty=3,col="red")
u            = sort(taupost)
lines(u, dgamma(u, alfa+n/2, beta + sum(y^2)/2),col="purple",lwd=2)

library(invgamma)
Vpost    = 1/taupost
hist(Vpost,freq=F,main="",xlab="V")
abline(v=V,lwd=3,lty=1,col="blue")
abline(v=mean(Vpost),lwd=3,lty=1,col="red")
abline(v=quantile(Vpost,0.025),lwd=3,lty=3,col="red")
abline(v=quantile(Vpost,0.975),lwd=3,lty=3,col="red")
curve(dinvgamma(x, alfa+n/2, beta + sum(y^2)/2), add=T, col="purple", lwd=2)

#comparando as médias e as variâncias a priori e a posteriori
c(alfa/beta, mean(taupost))
## [1] 2.000000 0.512957
c(alfa/(beta^2), var(taupost))
## [1] 2.0000000000 0.0005067143

3 Modelo binomial

\[\begin{eqnarray*} Y &\sim& Bin(n, \phi)\\ \phi &\sim& Beta(a,b) \end{eqnarray*}\]

##########################################################
# Modelo binomial com diferentes distribuições a priori. #
##########################################################

#--------------------------------------------#
# Distribuicao a Priori Beta: #
#--------------------------------------------#

#graficos da priori para diferentes valores de alpha e beta

curve(dbeta(x,1.5,4),0,1,lwd=2,lty=1,ylim=c(0,5),ylab=expression(p(theta)),
      xlab=expression(theta),col=1, 
      cex.lab=1.2, cex.axis=1.4, bty="n")
curve(dbeta(x,2,0.7),0,1,add=T,lwd=2,lty=6,col=2)
curve(dbeta(x,7,1.5),0,1,add=T,lwd=2,lty=3,col=3)
curve(dbeta(x,3,3)  ,0,1,add=T,lwd=2,lty=5,col=4)
curve(dbeta(x,1,1)  ,0,1,add=T,lwd=2,lty=4,col=5)
legend(.2,5,c('Be(1.5,4)','Be(2,0.5)','Be(7,1.5)','Be(3,3)','Be(1,1)'),
       lty=c(1,6,3,5,4),lwd=c(2,2,2,2),col=c(1,2,3,4,5))

#---------------------------------------------#
# Comparando Priori Beta com Posteriori Beta: #
# Y ~ Bin(n, phi)
# phi ~ Beta(a,b)
# phi |y ~ Beta(a+y, n+b-y)
#---------------------------------------------#
n   = 12
phi = 0.7
y   = rbinom(1, n, phi)
y
## [1] 9
medias_priori       = NULL
vari_priori         = NULL
medias_posteriori   = NULL
vari_posteriori     = NULL

#par(mfrow=c(2,2))

#priori 1
a                     = 1.5
b                     = 4
medias_priori[1]      = a/(a+b) #media
vari_priori[1]        = a*b/((a+b)^2 * (a+b+1) ) #variancia
medias_posteriori[1]  = (a+y)/(a+y+n+b-y) #media
vari_posteriori[1]    = (a+y)*(n+b-y)/((a+y+n+b-y)^2 * (a+y+n+b-y+1) ) #variancia
#curve(dbeta(x,alpha,beta),0,1,lty=1,ylim=c(0,5),ylab=expression(p(theta)),
#      xlab=expression(theta),col=4,lwd=2,
#dist a priori
curve(dbeta(x,a,b),0,1,lty=1,ylim=c(0,5),ylab=expression(p(phi)),
      xlab=expression(phi),col=4,lwd=2,
      main=paste("Média = ",round(medias_priori[1],2),
           " Var = ", round(vari_priori[1],2),sep=""))
#dist a posteriori
curve(dbeta(x, a+y, b+n-y),0,1,lty=6,ylim=c(0,5), col=4,add=T,lwd=2)
#dist a priori uniforme
curve(dbeta(x,1,1),0,1,lty=1,ylim=c(0,5),col=1,lwd=2,add=T)
#dist a posteriori
curve(dbeta(x,(y+1),(n-y+1)),0,1,col=1,add=T,lwd=2,lty=6)
abline(v=phi, lwd=2, col="red")

#legend(locator(n=1),c('Be(1.5,4))','Be(10.5,7)'),lty=c(1,6),lwd=c(1,1),col=c(4,4))
#legend(locator(n=1),c('Be(1,1))','Be(10,4)'),lty=c(1,6),lwd=c(1,1),col=c(1,1))

3.1 Mudando os hiperparâmetros da distribuição a priori

#priori 2
a                     = 2
b                     = 0.5
medias_priori[2]      = a/(a+b) #media
vari_priori[2]        = a*b/((a+b)^2 * (a+b+1) ) #variancia
medias_posteriori[2]  = (a+y)/(a+y+n+b-y) #media
vari_posteriori[2]    = (a+y)*(n+b-y)/((a+y+n+b-y)^2 * (a+y+n+b-y+1) ) #variancia
#curve(dbeta(x,alpha,beta),0,1,lty=1,ylim=c(0,5),ylab=expression(p(theta)),
#      xlab=expression(theta),col=4,lwd=2,
#dist a priori
curve(dbeta(x,a,b),0,1,lty=1,ylim=c(0,5),ylab=expression(p(phi)),
      xlab=expression(phi),col=4,lwd=2,
      main=paste("Média = ",round(medias_priori[1],2),
                 " Var = ", round(vari_priori[1],2),sep=""))
#dist a posteriori
curve(dbeta(x, a+y, b+n-y),0,1,lty=6,ylim=c(0,5), col=4,add=T,lwd=2)
#dist a priori uniforme
curve(dbeta(x,1,1),0,1,lty=1,ylim=c(0,5),col=1,lwd=2,add=T)
#dist a posteriori
curve(dbeta(x,(y+1),(n-y+1)),0,1,col=1,add=T,lwd=2,lty=6)
abline(v=phi, lwd=2, col="red")

3.2 Mudando novamente os hiperparâmetros da distribuição a priori

#priori 3
a                     = 7
b                     = 1.5
medias_priori[3]      = a/(a+b) #media
vari_priori[3]        = a*b/((a+b)^2 * (a+b+1) ) #variancia
medias_posteriori[3]  = (a+y)/(a+y+n+b-y) #media
vari_posteriori[3]    = (a+y)*(n+b-y)/((a+y+n+b-y)^2 * (a+y+n+b-y+1) ) #variancia
#curve(dbeta(x,alpha,beta),0,1,lty=1,ylim=c(0,5),ylab=expression(p(theta)),
#      xlab=expression(theta),col=4,lwd=2,
#dist a priori
curve(dbeta(x,a,b),0,1,lty=1,ylim=c(0,5),ylab=expression(p(phi)),
      xlab=expression(phi),col=4,lwd=2,
      main=paste("Média = ",round(medias_priori[1],2),
                 " Var = ", round(vari_priori[1],2),sep=""))
#dist a posteriori
curve(dbeta(x, a+y, b+n-y),0,1,lty=6,ylim=c(0,5), col=4,add=T,lwd=2)
#dist a priori uniforme
curve(dbeta(x,1,1),0,1,lty=1,ylim=c(0,5),col=1,lwd=2,add=T)
#dist a posteriori
curve(dbeta(x,(y+1),(n-y+1)),0,1,col=1,add=T,lwd=2,lty=6)
abline(v=phi, lwd=2, col="red")

3.3 Mudando novamente os hiperparâmetros da distribuição a priori

#priori 4
a                     = 3
b                     = 3
medias_priori[4]      = a/(a+b) #media
vari_priori[4]        = a*b/((a+b)^2 * (a+b+1) ) #variancia
medias_posteriori[4]  = (a+y)/(a+y+n+b-y) #media
vari_posteriori[4]    = (a+y)*(n+b-y)/((a+y+n+b-y)^2 * (a+y+n+b-y+1) ) #variancia
#curve(dbeta(x,alpha,beta),0,1,lty=1,ylim=c(0,5),ylab=expression(p(theta)),
#      xlab=expression(theta),col=4,lwd=2,
#dist a priori
curve(dbeta(x,a,b),0,1,lty=1,ylim=c(0,5),ylab=expression(p(phi)),
      xlab=expression(phi),col=4,lwd=2,
      main=paste("Média = ",round(medias_priori[1],2),
                 " Var = ", round(vari_priori[1],2),sep=""))
#dist a posteriori
curve(dbeta(x, a+y, b+n-y),0,1,lty=6,ylim=c(0,5), col=4,add=T,lwd=2)
#dist a priori uniforme
curve(dbeta(x,1,1),0,1,lty=1,ylim=c(0,5),col=1,lwd=2,add=T)
#dist a posteriori
curve(dbeta(x,(y+1),(n-y+1)),0,1,col=1,add=T,lwd=2,lty=6)
abline(v=phi, lwd=2, col="red")

3.4 Comparando

round(medias_priori,3)
## [1] 0.273 0.800 0.824 0.500
round(medias_posteriori, 3)
## [1] 0.600 0.759 0.780 0.667
round(vari_priori,3)
## [1] 0.031 0.046 0.015 0.036
round(vari_posteriori, 3)
## [1] 0.013 0.012 0.008 0.012
#-----------------------------------------------------------------#
# Uma função para comparar os gráficos do modelo binomial beta    #
# (Considerando agora uma amostra de tamanho m da dist. binomial) # 
# e uma função proporcional a fção de verossimilhança             #
#-----------------------------------------------------------------#

bin.beta = function(n,y,alpha,beta){
  alpha.post = alpha + sum(y)
  beta.post  = beta + sum(n-y)
  return(list(ap = alpha.post, bp = beta.post))
}

grafico.bin.beta = function(n,y,alpha,beta){
  par = bin.beta(n,y,alpha,beta)
  ap = par$ap
  bp = par$bp
  curve(dbeta(x,(sum(y)+1),(sum(n-y)+1)),0,1,col=1,xlab=expression(theta),ylab='',bty='n', lwd=2)
  curve(dbeta(x,alpha,beta),0,1,add=T,lty=1,col=4,lwd=2)
  curve(dbeta(x,ap,bp),0,1,add=T,lty=2,col=2,lwd=2)
  legend(0.2, 2.5,c('veross prop','priori','posteriori'),lty=1:3,bty='n',col=c(1,4,2))
}

par(mfrow=c(2,2))
grafico.bin.beta(12,9,1,1)
grafico.bin.beta(12,9,2,2)
grafico.bin.beta(12,9,1,3)
grafico.bin.beta(12,9,3,1)

3.5 Considerando agora uma amostra de tamanho m da distribuição binomial

#-----------------------------------------#
# Efeito do aumento do tamanho da amostra #
#-----------------------------------------#

alpha = 1
beta = 4

n = 10
p = 0.75

par(mfrow=c(2,2))

m = 1
y = rbinom(m,n,p)
grafico.bin.beta(n,y,alpha,beta)

m = 5
y = rbinom(m,n,p)
grafico.bin.beta(n,y,alpha,beta)

m = 10
y = rbinom(m,n,p)
grafico.bin.beta(n,y,alpha,beta)

m = 100
y = rbinom(m,n,p)
grafico.bin.beta(n,y,alpha,beta)