Inferência Bayesiana - GET/IME/UFF

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

1 Modelo Observacional Bernoulli

Considere o seguinte modelo observacional: \[\begin{eqnarray} Y_i|\theta &\stackrel{iid}\sim& Bern(\theta), \; i=1, \ldots, n. \end{eqnarray}\]

Vamos gerar então um conjunto de dados deste modelo da seguinte forma:

# Yi|theta ~ Bern(theta)
# gerando uma amostra deste modelo
n     = 500 #tamanho da amostra
theta = 0.2 #probabilidade de sucesso da Bernoulli
y     = rbinom(n, 1, theta) #gerando n valores da Bernoulli(theta)

Vamos agora avaliar os valores gerados através de um gráfico de barras:

#analisando a amostra de y gerada
par(mfrow=c(1,2))
barplot(table(y))
barplot(table(y)/sum(table(y)))

Agora suponha que \(\theta\) seja desconhecido. A função de verossimilhança é dada pela seguinte expressão: \[\begin{eqnarray} L(\theta; y_1, \ldots, y_n) &=& \prod_{i=1}^n{p(y_i|\theta)},%\\ %&=& \prod_{i=1}^n{\theta^{y_i}(1-\theta)^{1-y_i}}. \end{eqnarray}\] sendo \(p(y_i|\theta)\) a f.d. da distribuição \(Bern(\theta)\).

Então podemos avaliar como se comporta a função da log-verossimilhança nos dados gerados:

# Y1,...,Yn ~ Bern(theta)
# X = Y1+...+Yn ~ Bin(theta)

par(mfrow=c(1,2))
#analisando a log verossimilhança
thetas        = seq(0,1,l=1000)
logvero       = dbinom(sum(y), n, thetas, log=T)
plot(thetas,logvero, type="l", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="log L(theta; y1, ..., yn)")
abline(v=theta,col="red",lwd=2,lty=3) #plotando o valor verdadeiro

#analisando a log verossimilhança - DANDO UM ZOOM PARA VER MELHOR O COMPORTAMENTO DA FUNÇÃO
plot(thetas[which(thetas>0.1 & thetas<0.4)],logvero[which(thetas>0.1 & thetas<0.4)], type="l", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="log L(theta; y1, ..., yn)")
abline(v=theta,col="red",lwd=2,lty=3) #plotando o valor verdadeiro

1.1 Distribuição a priori Discreta

A priori, considere que o parâmetro \(\theta\) tem a seguinte função de distribuição: \[\begin{eqnarray} p(\theta) &=& 0,7 \times I(\theta=0,1) + 0,2 \times I(\theta=0,2) + 0,1 \times I(\theta=0,4), \end{eqnarray}\] sendo \(I(A)\) uma função indicadora que resulta em 1 quando a condição \(A\) é atendida e 0 caso contrário. Podemos representar graficamente esta crença da seguinte forma:

#analisando a distribuição a priori
thetas        = c(0.1, 0.2, 0.4)
dist_priori   = c(0.7, 0.2, 0.1)
plot(thetas,dist_priori,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="p(theta)")

Uma outra forma de visualizarmos graficamente o comportamento da distribuição a priori é através de um gráfico de barras de uma amostra da distribuição a priori:

#analisando a distribuição a priori através de uma amostra desta distribuição
theta_priori  = sample(thetas, 1000, replace=T, prob = dist_priori)
par(mfrow=c(1,2))
barplot(table(theta_priori))
barplot(table(theta_priori)/sum(table(theta_priori)))

Podemos também refazer o gráfico da função da log-verossimilhança considerando somente os 3 possíveis valores para o parâmetro \(\theta\) conforme imposto pela distribuição a priori:

# Y1,...,Yn ~ Bern(theta)
# X = Y1+...+Yn ~ Bin(theta)

#analisando a log verossimilhança
thetas        = c(0.1, 0.2, 0.4)
logvero       = dbinom(sum(y), n, thetas, log=T)
plot(thetas,logvero,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="log L(theta; y1, ..., yn)")

Dessa forma, temos que a distribuição a posteriori é \[\begin{eqnarray} p\left(\theta\left|x = \sum_{i=1}^n{y_i}\right.\right) &=& c \; 0,1^x 0,9^{n-x} \; 0,7 \; I(\theta=0,1) + \\ && c \; 0,2^x 0,8^{n-x} \; 0,2 \; I(\theta=0,2) + \\ && c \; 0,4^x 0,6^{n-x} \; 0,1 \; I(\theta=0,4). \end{eqnarray}\] E como \(p(\theta=0,1|x) + p(\theta=0,2|x) + p(\theta=0,4|x) = 1\), temos que \[\begin{eqnarray} c &=& \left( 0,1^x 0,9^{n-x} \; 0,7 + 0,2^x 0,8^{n-x} \; 0,2 + 0,4^x 0,6^{n-x} \; 0,1 \right)^{-1}. \end{eqnarray}\]

Logo, a representação gráfica da distribuição a posteriori é

x         = sum(y)
dist_post = thetas^x * (1-thetas)^(n-x) * dist_priori
dist_post = dist_post/sum(dist_post)
dist_post

## [1] 9.091279e-09 1.000000e+00 3.536460e-22

plot(thetas,dist_post,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="p(theta|y1, ..., yn)")

Uma outra forma de visualizarmos graficamente o comportamento da distribuição a posteriori é através de um gráfico de barras de uma amostra desta distribuição:

#analisando a distribuição a posteriori através de uma amostra desta distribuição
theta_posteriori = sample(c(0.1, 0.2, 0.4), 1000, replace=T, prob = dist_post)
par(mfrow=c(1,2))
barplot(table(theta_posteriori))
barplot(table(theta_posteriori)/sum(table(theta_posteriori)))

Comparando a função de verossimilhança, com a distribuição a priori com a distribuição a posteriori:

par(mfrow=c(1,3))
plot(thetas,dist_priori,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="p(theta)")
plot(thetas,logvero,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="log L(theta; y1, ..., yn)")
plot(thetas,dist_post,type="h", cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), bty="n", ylab="p(theta|y1, ..., yn)")

1.2 Distribuição a priori Beta

Vimos em aulas passadas que \[\begin{eqnarray} Y_i|\theta &\stackrel{iid}{\sim}& Bern(\theta)\\ \theta &\sim& Beta(a,b)\\ \theta|x = \sum_{i=1}^n{y_i} &\sim& Beta(a + x, b + n - x) \end{eqnarray}\]

Usando o mesmo dado gerado anteriormente, avaliaremos como se comportará a distribuição a posteriori do parãmetro \(\theta\) quando mudamos a crença inicial sobre isso. Usaremos uma distribuição a priori Beta informativa e para isso considere \(a=4\) e \(b=2\).

#criando uma grade de possiveis valores para theta
thetas        = seq(0,1,l=1000) #criando uma grade de valores para theta

#fixando os valores dos hiperparâmetros
a             = 4 
b             = 2 

#comparando o log da verossimilhança com a distribuição a priori com a distribuição a posteriori
lvero         = dbinom(sum(y), n, thetas, log=T) #calculando o log da verossimilhança
dist_priori   = dbeta(thetas, a, b) #calculando a distribuição a priori
dist_post     = dbeta(thetas, a + sum(y), b + n - sum(y)) #calculando a distribuição a posteriori

u             = c(dist_priori, dist_post)
plot(thetas, dist_priori, type="l", bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), ylab="", 
     ylim=c(min(u), max(u)))
lines(thetas, dist_post, lwd=2, col="red")
lines(thetas, exp(lvero)*10, lwd=2, col="purple") #multipliquei a lvero para melhorar a escala
abline(v=theta,lwd=2, col="blue", lty=3)

O que acontece se agora colocarmos uma distribuição a priori não informativa? Consideraremos 2 distrbuições neste caso:

• A distribuição uniforme, que é igual a distribuição Beta com parâmetros \(a=1\) e \(b=1\), para o parâmetro \(\theta\).

• A distribuição proposta por Jeffreys que, conforme visto em aula, será uma \(Beta(a=1/2, b=1/2)\).

#análise de sensibilidade do modelo quanto à distribuição a priori
a2             = 1 #hiperparâmetro da distribuição a priori
b2             = 1 #hiperparâmetro da distribuição a priori
#Beta(a=1, b=1) = Unif(0,1) => dist. a priori não informativa

a3             = 1/2 #hiperparâmetro da distribuição a priori
b3             = 1/2 #hiperparâmetro da distribuição a priori
#Beta(a=1/2, b=1/2) => dist. a priori de Jeffreys

dist_priori2   = dbeta(thetas, a2, b2) #calculando a distribuição a priori
dist_post2     = dbeta(thetas, a2 + sum(y), b2 + n - sum(y)) #calculando a distribuição a posteriori

dist_priori3   = dbeta(thetas, a3, b3) #calculando a distribuição a priori
dist_post3     = dbeta(thetas, a3 + sum(y), b3 + n - sum(y)) #calculando a distribuição a posteriori

par(mfrow=c(1,2))
u             = c(dist_priori, dist_priori2, dist_priori3)
u[abs(u)==Inf] = NA
plot(thetas, dist_priori, type="l", bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), ylab="p(theta)", 
     ylim=c(min(u, na.rm=T), max(u, na.rm=T)))
lines(thetas, dist_priori2, lwd=2, col="red")
lines(thetas, dist_priori3, lwd=2, col="purple") #multipliquei a vero para melhorar a escala
abline(v=theta,lwd=2, col="blue", lty=3)

u               = c(dist_post, dist_post2, dist_post3)
u[abs(u)==Inf]  = NA
plot(thetas, dist_post, type="l", bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     xlab=expression(theta), ylab="p(theta|y1,...,yn)", 
     ylim=c(min(u, na.rm=T), max(u, na.rm=T)))
lines(thetas, dist_post2, lwd=2, col="red")
lines(thetas, dist_post3, lwd=2, col="purple") #multipliquei a vero para melhorar a escala
abline(v=theta,lwd=2, col="blue", lty=3)

Note que mesmo com crenças iniciais tão diferentes, tivemos que as informações sobre o parâmetro baseadas nas distribuições a posteriori são bem semelhantes.

2 Modelo observacional Normal

Considere o seguinte modelo observacional: \[\begin{eqnarray} Y_i|\theta, \sigma^2 &\stackrel{iid}\sim& N(\theta, \sigma^2), \; i=1, \ldots, n. \end{eqnarray}\]

Vamos gerar então um conjunto de dados deste modelo e analisar a função log-verossimilhança:

################################# Modelo 3 #############################
# Y|theta, sigma2 ~ N(theta, sigma2)
########################################################################

#gerando uma amostra deste modelo
n       = 100 #tamanho da amostra
theta   = 10
sigma2  = 1
y       = rnorm(n, theta, sqrt(sigma2))

#analisando a amostra de y gerada
hist(y, freq=F, bty="n", ylab="densidade", xlab="y", main="", cex.lab=1.4,
     cex.axis=1.4)

#analisando o log da verossimilhança perfilada
thetas        = seq(0,20,l=1000) #criando uma grade de valores para theta
lvero         = NULL
for(j in 1:length(thetas))
{
  lvero[j] = sum( dnorm(y, thetas[j], sqrt(sigma2), log=T) ) #calculando o log da verossimilhança
}

plot(thetas, lvero, type="l", cex.lab=1.4, cex.axis=1.4, lwd=2, bty="n", xlab=expression(theta))
abline(v=theta,lwd=2,col="red",lty=3)

#analisando o log da verossimilhança perfilada
sigma2s       = seq(0.1,5,l=1000) #criando uma grade de valores para sigma2
lvero2       = NULL
for(j in 1:length(sigma2s))
{
  lvero2[j] = sum( dnorm(y, theta, sqrt(sigma2s[j]), log=T) ) #calculando o log da verossimilhança
}
plot(sigma2s, lvero2, type="l", cex.lab=1.4, cex.axis=1.4, lwd=2, bty="n", xlab=expression(sigma^2))
abline(v=sigma2,lwd=2,col="red",lty=3)

A priori, considere que \[\begin{eqnarray} \theta|\phi=\sigma^{-2} &\sim& N\left(\mu_0, \frac{1}{\phi\; C_0}\right)\\ \phi &\sim& G(a, b) \end{eqnarray}\]

Vimos, em sala de aula, que esta distribuição é conjugada com o modelo observacional normal. Podemos gerar uma amostra desta distribuição para avaliá-la, da seguinte forma:

################################# Modelo 3 #############################
# Y|theta, sigma2 ~ N(theta, sigma2)
# theta|phi=1/sigma2 ~ N(mu0, 1/(phi*C0))
# phi ~ G(a, b)
########################################################################

mu0           = 7 #hiperparâmetro da distribuição a priori
C0            = 1 #hiperparâmetro da distribuição a priori
a             = 5 #hiperparâmetro da distribuição a priori
b             = 2 #hiperparâmetro da distribuição a priori

# theta|phi=1/sigma2 ~ N(mu0, 1/(phi*C0))
# phi ~ G(a, b)
np          = 1000
phi_pr      = rgamma(np, a,b)
thetas_pr   = rnorm(np, mu0, sqrt(1/(phi_pr*C0))) 

par(mfrow=c(1,2))
hist(thetas_pr, freq=F, bty="n", ylab="densidade", xlab=expression(theta), main="", cex.lab=1.4,
     cex.axis=1.4)
abline(v=theta, lwd=2, lty=3, col="red")

hist(1/phi_pr, freq=F, bty="n", ylab="densidade", xlab=expression(sigma^2), main="", cex.lab=1.4,  cex.axis=1.4)
abline(v=sigma2, lwd=2, lty=3, col="red")

Note que a crença inicial não contemplou o valor verdadeiro de \(\theta\).

A posteriori, temos que \[\begin{eqnarray} \theta|\phi=\sigma^{-2}, y_1, \ldots, y_n &\sim& N\left(\frac{n \bar{y}+\mu_0 C_0}{n+C_0}, \frac{1}{\phi(n+ C_0)}\right)\\ \phi|y_1, \ldots, y_n &\sim& G\left(a+n/2, b+0,5\left[\sum_{i=1}^n{(y_i-\bar{y})^2} + n C_0 \frac{(\mu_0-\bar{y})^2}{n+C_0}\right]\right) \end{eqnarray}\] E então podemos obter uma amostra desta distribuição da seguinte forma:

# theta|y1, ..., yn,phi=1/sigma2 ~ N((n*ybarra+mu0*C0)/(n+C0), 1/(phi*(n+C0)))
# phi|y1, ..., yn, ~ G(a+n/2, b+0.5[sum((y-ybarra)^2) + n*C0*(mu0-ybarra)^2/(n+C0)])
phi_post      = rgamma(np, a+n/2,b+0.5*(sum((y-mean(y))^2) + n*C0*(mu0-mean(y))^2/(n+C0)))
thetas_post   = rnorm(np, (n*mean(y)+mu0*C0)/(n+C0), sqrt(1/(phi_post*(n+C0)))) 

par(mfrow=c(1,2))
hist(thetas_post, freq=F, bty="n", ylab="densidade", xlab=expression(theta), main="", cex.lab=1.4, cex.axis=1.4)
abline(v=theta, lwd=2, lty=3, col="red")

hist(1/phi_post, freq=F, bty="n", ylab="densidade", xlab=expression(sigma^2), main="", cex.lab=1.4, cex.axis=1.4)
abline(v=sigma2, lwd=2, lty=3, col="red")