0.1 Componentes da inferência Bayesiana

  • Distribuição a priori : utiliza a probabilidade como um meio de quantificar a incerteza sobre quantidades desconhecidas (variáveis), então temos \(f(\theta)\): distribuição a priori para o parâmetro \(\theta\);
  • Verossimilhança : relaciona todas as variáveis num modelo de probabilidade completo, então temos \(L(\theta | \boldsymbol{y})\): função de verossimilhança de \(\theta\) dado o conjunto de dados, vem diretametne de \(f(\boldsymbol{y}|\theta)\):
  • Distribuição a posteriori : quando observamos algumas variáveis (os dados), podemos usar a fórmula de Bayes para encontrar as distribuições de probabilidade condicionais para as quantidades de interesse não observadas, então temos \(f(\theta | \boldsymbol{y})\): distribuição a posteriori para o parâmetro \(\theta\).

Nosso principal objetivo é utilizar a distribuição a posteriori para a nossa tomada de decisões. Pelo Teorema de Bayes, temos:

    1. Caso discreto: neste caso, assumimos que \(\theta\) é uma variável aleatória discreta. \[ P(\theta=\theta_j|\boldsymbol{y})=\frac{P(\theta=\theta_j)f(\boldsymbol{y}|\theta)}{\sum \limits_j P(\theta=\theta_j) f(\boldsymbol{y}|\theta)}, \] onde \(\theta_j, j=1,2,\ldots\) são os valores que \(\theta\) pode assumir, ou seja, o espaço paramétrico de \(\theta\) é discreto,
    1. Caso contínuo: neste caso, assumimos que \(\theta\) é uma variável aleatória contínua. \[ f(\theta|\boldsymbol{y})=\frac{f(\theta)f(\boldsymbol{y}|\theta)}{\int \limits_\Theta f(\theta)f(\boldsymbol{y}|\theta) d\theta}, \] onde \(\Theta\) é o espaço paramétrico de \(\theta\), o espaço paramétrico de \(\theta\) é contínuo

Observações:

  • O caso contínuo é mais comumente utilizado na estatística Bayesiana,
  • A distribuição a priori também pode ser denotada por \(p(\theta)\) ou \(\pi(\theta)\), assim como a distribuição a posteriori denotada por \(p(\theta | \boldsymbol{y})\) ou \(\pi(\theta | \boldsymbol{y})\),
  • Em geral, não é necessário efetuar o cálculo do denominador \(\int \limits_\Theta f(\theta)L(\theta | \boldsymbol{y}) d\theta\) pois se trata de uma constante que não depende de \(\theta\), então temos \[ f(\theta|\boldsymbol{y})\propto \frac{f(\theta)L(\theta | \boldsymbol{y})}{\int \limits_\Theta f(\theta)L(\theta | \boldsymbol{y}) d\theta}, \] donde o símbolo \(\propto\) significa “é proporcional a”,
  • As distribuições a priori podem ser de vários tipos e características:
  • Quanto à propriedade de integrabilidade: existem prioris próprias ou impróprias,
  • Quanto ao nível de informação: existem prioris não informativas ou informativas,
  • Quanto a depender ou não da amostra (dos dados): existem prioris subjetivas ou objetivas

0.1.1 Exercícios - continuação

    1. Seja \(Y\) uma variável aleatória com distribuição Binomial \(Y \sim \text{ Binomial }(n,p)\).
      1. Qual é a estimativa de máxima verossimilhança de \(p\)?
      1. Se \(p\) tem distribuição a priori Beta com parâmetros conhecidos \(a\) e \(b\) então qual é a distribuição a posteriori para \(p\)?
      1. Segundo o item \(b)\), qual é a média a posteriori para \(p\)?

Dicas para o item \(b)\) \[P(Y=y)=?\] Qual é a distribuição de Y (os dados)? \(f(p)=\)? Qual é a distribuição a priori para o parâmetro \(p\)? \(f(p|\boldsymbol{y})=\) ? Qual é a distribuição de \(p\) condicionada aos dados?

    1. Foram gerados \(10\) valores da distribuição de Poisson com taxa \(\lambda=2\), através do seguinte código no software R:
set.seed(15052017)
lambda=2  #este é o valor verdadeiro de lambda
n=10
x=rpois(n,lambda)
x
##  [1] 3 2 0 1 4 4 3 0 3 3
    1. Obtenha a estimativa de máxima verossimilhança para \(\lambda\);
    1. Considerando a distribuição a priori Gamma com média \(1\) e variância \(5\), obtenha a distribuição a posteriori para \(\lambda\);
    1. Segundo o item \(b)\), qual é a média a posteriori para \(\lambda\)?
    1. Segundo o item \(b)\), qual é a variância a posteriori para \(\lambda\)? Veremos mais tarde que os exercícios 4 e 5 envolvem distribuições a priori conjugadas.

Resolução do Exercício 4

Você vai precisar de

  • Apêndice B do Mood - distribuições ;

    1. Considerando Y com distribuição de Binomial: \[ \begin{array}{lll} P(Y=y)=\left(\begin{array}{ll} n\\y \end{array}\right) \times p^y (1-p)^{n-y} \text{ então }\\ \hat{p}=\frac{y}{n}, \text{ ou seja, a e.m.v. de }p \text{ é igual a contagem do número de sucessos sobre o número de tentativas} \end{array} \]

  • Passo a passo

    1. a posteriori para \(p\) vem do caso 3 de prioris conjugadas (veja seção 2.2):

\[ p \vert \boldsymbol{y} \sim \text{Beta}(a+y,b+n-y), \text{ ou seja}\\ f(p \vert \boldsymbol{y})= \frac{1}{B(a+y,b+n-y)}\times p^{a+y-1} \times (1-p)^{b+n-y-1} \times I_{(0,1)}(p) \] - c) média a posteriori - Olhar no apêndice do Mood a fórmula da média: \[ \text{E}(p \vert \boldsymbol{y})= \frac{a+y}{(a+y)+(b+n-y)}\\ =\frac{a+y}{a+b+n} \] - Aplicação com números: Em 16 ensaios de Bernoulli, foram obtidos 7 sucessos. Considere uma priori Beta de parâmetros \(a=0,5\) e \(b=0,5\).

  • Calcule a e.m.v. de \(p\) e mostre graficamente.
n=16
y=7
logL=function(p){
L=dbinom(y,size=n,prob=p)
log(L)
}
optimize(logL, interval=c(0.01,0.99), maximum=TRUE) #encontra o ponto de máximo
## $maximum
## [1] 0.4374995
## 
## $objective
## [1] -1.620156
emv=y/n
print(paste0("emv= ",emv))
## [1] "emv= 0.4375"
p=seq(0.01,0.99,0.01)
temp <- data.frame(p = p, logL = logL(p))
datatable(temp)
ggplot(data = temp, aes(x = p, y = logL)) + geom_line() + stat_peaks(col = "red")

  • Implemente os gráficos da priori, verossimilhança e posteriori:
  • Para colocar as três funções no mesmo gráfico, é necessário aproximar a curva da verossimilhança para o formato de uma densidade
    • com soma igual a 1, apenas para fins ilustrativos de construção de gráfico;
    • Sendo assim, para cada valor calculado da função de verossimilhança, divide-se por uma constante que é dada pela soma de todos os valores da função - denotado por “vero_normalizado” na rotina abaixo:
a=0.5
b=0.5
p=seq(0.01,0.99,0.01)
pri=dbeta(p,shape1=a,shape2=b)

#Aproxima a verossimilhança para o formato de uma densidade (com soma igual a 1, apenas para fins ilustrativos de construção de gráfico)
intervalos=0.01
vero=dbinom(y,size=n,prob=p)
vero_normalizado=vero/sum(vero*intervalos)
pos=dbeta(p,shape1=a+y,shape2=b+n-y)
plot(p,pri,type='l',xlab=expression(p),ylab=expression(f(p)),ylim=c(0,4))
lines(p,vero_normalizado,type='l',col=2) 
lines(p,pos,type='l',col=3)
legend("topright",c(expression("priori "*f(p)),expression("verossimilhança "*L(p*"|"*bold(y))),expression("posteriori "*f(p*"|"*bold(y)))),col=c(1,2,3),bty = "n",lty=c(1,1,1))

Resolução do Exercício 5

Você vai precisar de

  • Apêndice B do Mood - distribuições ;

    1. Considerando uma amostra de tamanho \(n\) com distribuição de Poisson: \[ \begin{array}{lll} P(X=x)=\frac{e^{-\lambda} \lambda^x}{x!} \text{ então}\\ \hat{\lambda}=\bar{x}, \text{ou seja, a e.m.v. de } \lambda \text{ é igual à média amostral} \end{array} \]

  • Passo a passo

  • Aplicando nos dados: \(\hat{\lambda}=2.3\), pois

x=c(3,2,0,1,4,4,3,0,3,3)
lambda_hat=mean(x) #estimativa de maxima verossimilhança
lambda_hat
## [1] 2.3
  • Visualizando a log-verossimilhança numericamente graficamente:
n=length(x)
logL=function(lambda){
L=exp(-n*lambda)*lambda^(sum(x))/prod(factorial(x))
log(L)
}
optimize(logL, interval=c(0.01,10), maximum=TRUE) #encontra o ponto de máximo
## $maximum
## [1] 2.300012
## 
## $objective
## [1] -18.05938
lambda=seq(0.1,5,0.1) #só assume valores positivos
temp <- data.frame(lambda=lambda, logL = logL(lambda))
datatable(temp)
ggplot(data = temp, aes(x = lambda, y = logL)) + geom_line() + stat_peaks(col = "red")

    1. Seja \(\lambda \sim \text{Gamma}(a,r)\), então podemos encontrar os valores de \(a\) e \(r\) baseados na média e variância a priori: \[ f(\lambda)=\frac{a^r}{\Gamma[r]}\lambda^{r-1}e^{-a\lambda} I_{(0,\infty)}(\lambda) \text{ com } \left\{ \begin{array}{lll} \text{E}(\lambda)=\frac{r}{a}=1 \\ \text{VAR}(\lambda)=\frac{r}{a^2}=5 \end{array} \right. \Rightarrow a=r=\frac{1}{5} \\ \] Sendo assim, a posteriori para \(\lambda\) vem do caso 2 de prioris conjugadas (veja seção 2.2):
paste0("somatório de x= ",sum(x))
## [1] "somatório de x= 23"

\[ \lambda \vert \boldsymbol{y} \sim \text{Gamma}(a+n,r+\sum_{i=1}^n x_i), \text{ ou seja, }\\ \text{Como } a=r=\frac{1}{5}, n=10, \text{ e }\sum_{i=1}^n x_i=23\\ \lambda \vert \boldsymbol{y} \sim \text{Gamma}\left(\frac{51}{5},\frac{116}{5}\right) \] - c) e d) média e variância a posteriori - Olhar no apêndice do Mood as fórmulas da média & variância: \[ \begin{array}{lll} \text{E}(\lambda \vert \boldsymbol{x})&=& \frac{\frac{116}{5}}{\frac{51}{5}}\approx 2.2745\\ \text{VAR}(\lambda \vert \boldsymbol{x})&=&\frac{\frac{116}{5}}{\left(\frac{51}{5}\right)^2}\approx 0.2230 \end{array} \] - Graficamente: Fazendo os graficos da priori, verossimilhança e posteriori

a=1/5
r=1/5

pri_lambda=dgamma(lambda,shape=r, scale=1/a) #na parametrização do R, o parâmetro de  escala (a) é invertido

#Aproxima a verossimilhança para o formato de uma densidade (com soma igual a 1, apenas para fins ilustrativos de construção de gráfico)
intervalos=0.1
L_lambda=exp(-n*lambda)*lambda^(sum(x))/prod(factorial(x))
L_lambda_normalizado=L_lambda/sum(L_lambda*intervalos)
pos_lambda=dgamma(lambda,shape=r+sum(x), scale= 1/(a+n)) #na parametrização do R, o parâmetro de escala beta é invertido
plot(lambda,pri_lambda,type='l',xlab=expression(lambda),ylab=expression(f(lambda)),ylim=c(0,1.5))
lines(lambda,L_lambda_normalizado,type='l',col=2) 
lines(lambda,pos_lambda,type='l',col=3)
legend("topright",c(expression("priori "*f(lambda)),expression("verossimilhança "*L(lambda*"|"*bold(x))),expression("posteriori "*f(lambda*"|"*bold(x)))),col=c(1,2,3),bty = "n",lty=c(1,1,1))

Figura 1: Priori, Verossimilhança e Posteriori

Exemplo 1: Ensaios de Bernoulli com distribuição a priori discreta.

  • Uma determinada droga tem taxa de resposta \(\theta\) podendo assumir os seguintes valores a priori: \(0,2\); \(0,4\), \(0,6\) ou \(0,8\), sendo cada um dos valores com mesma probabilidade de ocorrência. Do resultado de uma amostra unitária, obtivemos sucesso. Como nossa crença pode ser revisada? Podemos representar o problema através de uma tabela.

  • Atribui-se priori uniforme discreta para a proporção \(\theta\);

  • E a distribuição a posteriori para esta proporção será uniforme discreta.

Resolução: Você vai precisar de

  • Apêndice B do Mood - distribuições ;
  • Fórmula de Bayes (da seção 1.3) para o caso discreto (pois \(\theta\)=0,2 ou 0,4 ou 0,6 ou 0,8);
  • sem dispensar o termo do denominador - pois é caso discreto. \[ P(\theta=\theta_j \vert y=1) = \frac{P(\theta=\theta_j)P(y=1 \vert \theta=\theta_j)}{\sum_{i=1}^4 P(\theta=\theta_j)P(y=1 \vert \theta=\theta_j )}\\ \text{ mesmo raciocínio do Teorema de Bayes!} \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Bernoulli), ou seja, temos uma amostra de tamanho \(1\) de uma distribuição de Bernoulli com proporção igual a \(\theta\), e pelo enunciado \(y=1\).

Passo I: atribuir priori para \(\theta\), sendo \(\theta\)=0,2 ou 0,4 ou 0,6 ou 0,8:

require(kableExtra)
## Loading required package: kableExtra
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
theta = c(0.2,0.4,0.6,0.8)
priori_theta = c(0.25,0.25,0.25,0.25)

tabela <- data.frame(
          theta,
          priori_theta)

names(tabela) = c("$\\theta_j$","$P(\\theta=\\theta_j)$")

tabela=rbind(tabela,c("Total",1))

 tabela %>%
  kbl(escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
     kable_styling(bootstrap_options = c("striped","hold_position"),
                full_width = F)
\(\theta_j\) \(P(\theta=\theta_j)\)
0.2 0.25
0.4 0.25
0.6 0.25
0.8 0.25
Total 1

Passo II: atribuir distribuição Bernoulli para os dados(amostra de tamanho unitário): \[ Y ~ \sim \text{ Bernoulli }(\theta) \text{ então }\\ P(Y=1) = \theta \]

theta = c(0.2,0.4,0.6,0.8)
L_theta = theta

tabela <- data.frame(
      theta,
      L_theta)
names(tabela) = c("$\\theta_j$","$P(Y=1 \\vert \\theta = \\theta_j)$")

 tabela %>%
  kbl(escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
     kable_styling(bootstrap_options = c("striped","hold_position"),
                full_width = F)
\(\theta_j\) \(P(Y=1 \vert \theta = \theta_j)\)
0.2 0.2
0.4 0.4
0.6 0.6
0.8 0.8

Passo III:Aplicar a fórmula

theta = c(0.2,0.4,0.6,0.8)
priori_theta = c(0.25,0.25,0.25,0.25)
L_theta = theta
post_theta=priori_theta*L_theta/sum(priori_theta*L_theta)

tabela <- data.frame(
          theta,
          priori_theta,
          L_theta,
          post_theta)

names(tabela) = c("$\\theta_j$","$P(\\theta=\\theta_j)$","$P(Y=1 \\vert \\theta = \\theta_j)$","$P(\\theta=\\theta_j \\vert y=1)$")

tabela=rbind(tabela,c("Total",1,"",1))

 tabela %>%
  kbl(escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
     kable_styling(bootstrap_options = c("striped","hold_position"),
                full_width = F)
\(\theta_j\) \(P(\theta=\theta_j)\) \(P(Y=1 \vert \theta = \theta_j)\) \(P(\theta=\theta_j \vert y=1)\)
0.2 0.25 0.2 0.1
0.4 0.25 0.4 0.2
0.6 0.25 0.6 0.3
0.8 0.25 0.8 0.4
Total 1 1

Ilustração:

plot(theta,priori_theta,xlab=expression(theta),ylab=expression(p(theta)),type='p',col=2,ylim=c(0,1),pch=2)
lines(theta,L_theta,type='p',col=3,pch=3)
lines(theta,post_theta,type='p',col=4,pch=4)
legend("topleft",c("priori","verossimilhança","posteriori"),col=c(2,3,4),pch=c(2,3,4),bty = "n")

Exemplos de aplicação. Material de Inferência Estatística Ricardo Ehlers pag. 7:

  1. Considere \(X_1, X_2, \ldots,X_n\) uma amostra aleatória de tamanho \(n\) da distribuição exponencial com parâmetro \(\lambda\).
    1. Encontre a expressão para a estimativa de máxima verossimilhança.

Resposta: Considerando a parametrização para a exponencial: \[ \begin{array}{lll} f(x)=\lambda e^{-\lambda x} \text{ então}\\ \hat{\lambda}=\frac{1}{\bar{x}}, \text{ou seja, a e.m.v. de } \lambda \text{ é dada pelo inverso da média amostral} \end{array} \] - Passo a passo

    1. Faça uma aplicação mostrando o gráfico da função de verossimilhança de \(\lambda\) baseado em uma amostra:
set.seed(05102020)
lambda=5  #valor utilizado para gerar os dados
n=30
x=rexp(n,lambda)
x_bar=mean(x)
paste0("n= ",n)
## [1] "n= 30"
emv=1/x_bar
paste0("média amostral=",x_bar)
## [1] "média amostral=0.148695362229909"
paste0("emv=",emv)
## [1] "emv=6.72515931232494"
require(tidyverse)
require(ggpmisc)
require(DT)
logL=function(lambda){
L=lambda^n *exp(-lambda*sum(x))
log=log(L)
log
}
optimize(logL, interval=c(0.1, 10), maximum=TRUE) #encontra o ponto de máximo
## $maximum
## [1] 6.725157
## 
## $objective
## [1] 27.17567
lambda=seq(0.1,10,0.1)
temp <- data.frame(lambda = lambda, logL = logL(lambda))
datatable(temp)
ggplot(data = temp, aes(x = lambda, y = logL)) + geom_line() + stat_peaks(col = "red")