require(gtools)
require(tidyverse)
require(ggforce)  
options(kableExtra.latex.load_packages = TRUE)
require(kableExtra)
require(mosaicCalc)
require(gridExtra)
require(numDeriv)
require(DT)
require(ggpmisc) #para o statpeaks

0.1 Prioris conjugadas

Uma família de distribuições a priori é conjugada se as distribuições a posteriori pertencem a esta mesma família de distribuições.

0.1.1 Casos principais de prioris conjugadas

Caso 1: Quando os dados têm distribuição normal com média desconhecida e variância conhecida

  • Atribui-se priori normal para a média \(\mu\);
  • E a distribuição a posteriori para esta média será normal.

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 contínuo (pois \(-\infty < \mu < \infty\));
  • dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(\mu\); \[ f(\mu \vert \boldsymbol{y})\propto f(\mu)L(\mu \vert \boldsymbol{y}) \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Normal), ou seja, temos uma amostra de tamanho \(n\) i.i.d. (independentes e identicamente distribuídos) de uma distribuição de Normal com média igual a \(\mu\) e variância igual a \(\sigma^2\), então \(\boldsymbol{y}\) é um vetor de tamanho \(n\).

Passo I: atribuir priori para \(\mu\), sendo \(-\infty<\mu\infty\),

\[ \begin{array}{ll} & \mu \sim \text{Normal}(m_0,s_0^2):\\ &f(\mu)=\frac{1}{\sqrt{2 \pi} s_0} \exp \left[-\frac{1}{2 s_0^2} (\mu-m_0)^2 \right] \end{array}\\ \text{considere } m_0 \text{ e } s_0^2 \text{ os parâmetros da média e variância da distribuição a priori, respectivamente} \] - com \(m_0\) e \(s_0^2\) conhecidos; - A distribuição a priori nos traz o conhecimento a priori sobre a média \(\mu\); - Se temos pouca informação a respeito de \(\mu\), podemos fixar a média \(m_0\) e atribuir uma variância \(s_0^2\) grande; - Se temos muita informação a respeito de \(\mu\), podemos fixar a média \(m_0\) e atribuir uma variância \(s_0^2\) pequena.

Passo II: atribuir distribuição Normal para os dados \[ Y_1, Y_2,Y_3,\ldots,Y_n \sim \text{Normal}(\mu,\sigma^2) \text{ então }\\ \begin{array}{lll} f_{\boldsymbol{Y}}(\boldsymbol{y})&=&\prod_{i=1}^n\frac{1}{\sqrt{2 \pi} \sigma} \exp\left[-\frac{1}{2\sigma^2} (y_i-\mu)^2 \right]\\ &=& \left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n \exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2\right]\\ &=&L(\mu \vert \boldsymbol{y}) \end{array} \] distribuição dos dados = função de verossimilhança de \(\mu\) dado os dados

  • A função de verossimilhança nos traz toda a informação disponível na amostra (nos dados);
  • com \(\mu\) desconhecido e \(\sigma^2\) conhecido;

Passo III:Aplicar a fórmula - pode-se dispensar os termos constantes \[ \begin{array}{lll} f(\mu \vert \boldsymbol{y}) &\propto& f(\mu) L(\mu \vert \boldsymbol{y})\\ &\propto& \exp \left[-\frac{1}{2} \left( \frac{1}{s_0^2} (\mu-m_0)^2 + \frac{1}{\sigma^2} \sum \limits_{i=1}^n (y_i-\mu)^2 \right) \right],\\ &\text{ Então }& \mu | \boldsymbol{y} \sim N\left(\frac{\frac{m_0}{s_0^2}+\frac{n \overline{y}}{\sigma^2}}{\frac{1}{s_0^2}+\frac{n}{\sigma^2}},\frac{1}{\frac{1}{s_0^2}+\frac{n}{\sigma^2}}\right) \end{array} \] - o símbolo \(\propto\) significa “é proporcional a”, ou seja, todos os termos multiplicativos que não dependem de \(\mu\) podem ser desconsiderados na fórmula. - A demonstração pode ser encontrada em Box & Tiao (1973)

0.2 Prioris Conjugadas - continuação

Exemplo 1 Box & Tiao (1973) Os físicos A e B desejam determinar uma quantidade física \(\mu\). O físico A tem mais experiência nesta área e especifica sua priori como \(\mu \sim N(900,20^2)\). O físico B tem pouca experiência e especifica uma priori muito mais incerta em relação à posição de \(\mu\): \(\mu\sim N(800,80^2)\). Faz-se então uma medição \(X\) de \(\mu\) em laboratório com um aparelho calibrado com distribuição amostral \(X \vert \mu \sim N(\mu,40^2)\) e observa-se \(X=850\). Este exemplo corresponde ao “Caso 1)” de prioris conjugadas.

Explorando o exemplo e obtendo a posteriori:

  • para o físico A:

    • distribuição a priori: \(\mu \sim N(900,40^2)\);
    • \(P(860<\mu<940)\approx 0,95\): o intervalo que abrange 95% dos valores é mais estreito.
    • distribuição a posteriori: \(\mu \vert \boldsymbol{x} \sim N(890,320)\).
    • a variância do nosso parâmetro diminuiu \(\rightarrow\) significa que ganhamos informação com os dados observados;

  • para o físico B:

    • distribuição a priori: \(\mu \sim N(800,80^2)\);
    • \(P(640<\mu<960)\approx 0,95\): o intervalo que abrange 95% dos valores é mais largo.
    • distribuição a posteriori: \(\mu \vert \boldsymbol{x} \sim N(840,1280)\).
    • a variância de \(\mu\) era igual a \(6400\), e passou a ser igual a \(1280\) \(\rightarrow\) agregamos informação da amostra;

  • A distribuição a posteriori representa um compromisso entre a distribuição a priori e a verossimilhança. Além disso, como as incertezas iniciais são bem diferentes, o mesmo experimento fornece muito pouca informação adicional para o físico A enquanto que a incerteza do físico B foi bastante reduzida.

  • Cálculos para este exemplo:

require(DT)
#Informações a priori
m_0=c(900,800)
s_0=c(20,80)

quantis_A=qnorm(c(0.025,0.975),mean=m_0[1],sd=s_0[1]) 
quantis_B=qnorm(c(0.025,0.975),mean=m_0[2],sd=s_0[2]) 
quantis=data.frame(prob=c(0.025,0.975),quantis_A,quantis_B)
datatable(quantis,caption="Distribuição a priori")
#Amostra
x=850
n=1
sigma=40
sigma2=sigma^2
x_bar=mean(x)

#Informações a posteriori
media=(m_0/s_0^2+n*x_bar/sigma2)/(1/s_0^2+n/sigma2)
variancia=1/(1/s_0^2+n/sigma2)
posteriori=data.frame(rbind(media,variancia))
names(posteriori)=c("Fisico A","Fisico B")
datatable(posteriori,caption="Distribuição a posteriori")
  • Representação Gráfica:
mu=seq(500,1000,1) 
priori_A=dnorm(mu,mean=m_0[1],sd=s_0[1])
priori_B=dnorm(mu,mean=m_0[2],sd=s_0[2])
L_mu=dnorm(x,mean=mu,sd=40)

#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=1
L_mu_normalizado=L_mu/sum(L_mu*intervalos)
posteriori_A=dnorm(mu,mean=media[1],sd=sqrt(variancia[1]))
posteriori_B=dnorm(mu,mean=media[2],sd=sqrt(variancia[2]))
par(mfrow=c(1,2)) 
plot(mu,L_mu_normalizado,type='l',xlab=expression(mu),ylab=expression(f(mu)),main="Fisico A",ylim=c(0,0.03),col=2)
lines(mu,priori_A,type='l',col=1)
lines(mu,posteriori_A,type='l',col=3)
legend("topleft",c(expression("priori: "*mu*"~"*N(900,40^2)),expression("Veross.: "*X*"~"*N(mu,40^2)),expression("poster.: "*mu*"~"*N(890,320))),col=c(1,2,3),lty=c(1,1,1),bty = "n")

plot(mu,L_mu_normalizado,type='l',xlab=expression(mu),ylab=expression(f(mu)),main="Fisico B",ylim=c(0,0.03),col=2)
lines(mu,priori_B,type='l',col=1)
lines(mu,posteriori_B,type='l',col=3)
legend("topleft",c(expression("priori: "*mu*"~"*N(800,80^2)),expression("Veross.: "*X*"~"*N(mu,40^2)),expression("poster.: "*mu*"~"*N(840,1280))),col=c(1,2,3),lty=c(1,1,1),bty = "n")

Figura 1: Gráficos das três funções conjuntamente: priori, verossimilhança e posteriori

Caso 2: Quando os dados têm distribuição de Poisson

  • Atribui-se priori Gmma para a taxa \(\lambda\);
  • E a distribuição a posteriori para esta taxa será Gamma.

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 contínuo (pois \(0<\lambda<\infty\));
  • dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(\lambda\); \[ f(\lambda \vert \boldsymbol{y})\propto f(\lambda)L(\lambda \vert \boldsymbol{y}) \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Poisson), ou seja, temos uma amostra de tamanho \(n\) i.i.d. (independentes e identicamente distribuídos) de uma distribuição de Poisson com taxa \(\lambda\), então \(\boldsymbol{y}\) é um vetor de tamanho \(n\).

Passo I: atribuir priori para \(\lambda\), sendo \(0<\lambda<\infty\),

\[ \begin{array}{lll} \lambda \sim \text{Gamma}(a,r) \text{ então}\\ f(\lambda)=\underbrace{\frac{a^r}{\Gamma[r]} \lambda^{r-1} e^{-a \lambda} I_{(0,\infty)}(\lambda)}_{\text{considere } a \text{ o primeiro parâmetro da Gamma, para não confundir com }\lambda \text{ na fórmula}} \end{array} \] Passo II: atribuir distribuição Poisson para os dados \[ Y_1, Y_2,Y_3,\ldots,Y_n \sim \text{Poisson}(\lambda) \text{ então }\\ \begin{array}{lll} \underbrace{P(\boldsymbol{Y}=\boldsymbol{y})}_{\text{vetor }\boldsymbol{Y}\text{ ser igual a "vetorzinho" }\boldsymbol{y}\text{, de valores observados}}&=&\prod_{i=1}^n \frac{e^{-\lambda}\lambda^{y_i}}{y_i!}\times \underbrace{\prod_{i=1}^n I_{\{0,1,\ldots\}} (y_i)}_{\text{produto de funções indicadoras}}\\ &=& \frac{e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i}}{\prod_{i=1}^n y_i!} \times I_{\{0,1,\ldots\}}(\prod_{i=1}^n y_i)\\ &=&L(\lambda \vert \boldsymbol{y}) \end{array}\\ \text{distribuição dos dados = função de verossimilhança de } \lambda \text{ dado os dados} \] Passo III:Aplicar a fórmula - pode-se dispensar os termos constantes \[ \begin{array}{lll} f(\lambda \vert \boldsymbol{y})&\propto& f(\lambda)L(\lambda \vert \boldsymbol{y})\\ &\propto& \frac{a^r}{\Gamma[r]} \lambda^{r-1} e^{-a \lambda} I_{(0,\infty)}(\lambda) \times \frac{e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i}}{\prod_{i=1}^n y_i!}\\ &\propto& \lambda^{r+\sum_{i=1}^n y_i-1} e^{-(a+n)\lambda} I_{(0,\infty)}(\lambda)\\ &\text{ Então }& \lambda \vert \boldsymbol{y} \sim \text{Gamma}(a+n,r+\sum_{i=1}^n y_i) \text{ é só incluir as constantes na fórmula:} \\ f(\lambda \vert \boldsymbol{y})&=& \frac{(a+n)^{r+\sum_{i=1}^n y_i}}{\Gamma[r+\sum_{i=1}^n y_i]} \lambda^{r+\sum_{i=1}^n y_i-1} e^{-(a+n) \lambda} I_{(0,\infty)}(\lambda) \end{array} \]

  • Na inferência clássica, nós procedemos com o método usual de maximizar \(\log \left(L(\lambda|\boldsymbol{y})\right)\) com respeito a \(\lambda\), pois estamos tratando de um parâmetro fixo e desconhecido;
  • Na inferência Bayesiana, ao invés de encontrar a estimativa de máxima verossimilhança, nós estamos interessados na distribuição a posteriori, pois estamos tratando de um parâmetro desconhecido cuja distribuição a priori é a nossa crença a priori antes de observar os dados.

Caso 3: Quando os dados têm distribuição Binomial

  • Atribui-se priori Beta para a proporção \(p\);
  • E a distribuição a posteriori para esta proporção será Beta.

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 contínuo (pois 0<p<1);
  • dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(p\); \[ f(p \vert \boldsymbol{y})\propto f(p)L(p \vert \boldsymbol{y}) \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Binomial), ou seja, quando realizados de \(n\) ensaios de Bernoulli independentes com probabilidade igual a \(p\), verificou-se \(y\) sucessos.

Passo I: atribuir priori para \(p\), sendo \(0<p<1\), \[ \begin{array}{lll} p \sim \text{Beta}(a,b) \text{ então }\\ f(p)=\frac{1}{B(a,b)} \times p^{a-1} (1-p)^{b-1} I_{(0,1)}(p) \end{array} \] Passo II: atribuir distribuição Binomial para os dados \[ \begin{array}{lll} Y \sim \text{Binominal}(n,p) \text{ então }\\ P(Y=y)=\underbrace{ \left(\begin{array}{ll} n\\y \end{array}\right) \times p^y (1-p)^{n-y}=L(p \vert \boldsymbol{y})}_{\text{distribuição dos dados = função de verossimilhança de }p\text{ dado os dados}} \end{array} \] Passo III: Aplicar a fórmula - pode-se dispensar os termos constantes \[ \begin{array}{lll} f(p \vert \boldsymbol{y})&\propto& f(p)L(p \vert \boldsymbol{y})\\ &\propto& \frac{1}{B(a,b)} \times p^{a-1} (1-p)^{b-1} I_{(0,1)}(p) \times \left(\begin{array}{ll} n\\y \end{array}\right) \times p^y (1-p)^{n-y}\\ &\propto& p^{a+y-1} \times (1-p)^{b+n-y-1} \times I_{(0,1)}(p)\\ &\text{ Então }& p \vert \boldsymbol{y} \sim \text{Beta}(a+y,b+n-y) \text{ é só incluir as constantes na fórmula:}\\ 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) \end{array}\\ \] - Também pode-se mostrar que sendo \(Y_1,Y_2,\ldots,Y_N \sim\) Binomial \((n,p)\) a priori conjugada é Beta - com parâmetros atualizados para N experimentos de Bernoulli.

  • Veja abaixo um exemplo no R para este cenário - olhe o código
N = 20 #Número de experimentos
n = 10 #Número de ensaios de Bernoulli para cada experimento
p = 0.4 #probabilidade de sucesso para cada ensaio de Bernoulli
set.seed(18102020)
y = rbinom(N, n, p) #vetor y
y 
##  [1] 4 4 5 5 7 5 5 3 5 5 4 6 6 4 4 5 3 5 7 5

Caso 4: Quando os dados têm distribuição normal com média conhecida e variância desconhecida

  • Atribui-se priori Gamma invertida para a variância \(\sigma^2\);
  • E a distribuição a posteriori para esta variância será Gamma invertida.

Resolução: Você vai precisar de

  • Apêndice A do material do prof. Ricardo Ehlers - distribuições ;
  • Fórmula de Bayes (da seção 1.3) para o caso contínuo (pois \(0<\sigma^2<\infty\));
  • dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(\sigma^2\); \[ f(\sigma^2 \vert \boldsymbol{y})\propto f(\sigma^2)L(\sigma^2 \vert \boldsymbol{y}) \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Normal), ou seja, temos uma amostra de tamanho \(n\) i.i.d. (independentes e identicamente distribuídos) de uma distribuição de Normal com média igual a \(\mu\) e variância igual a \(\sigma^2\), então \(\boldsymbol{y}\) é um vetor de tamanho \(n\).

Passo I: atribuir priori para \(\sigma^2\), sendo \(0<\sigma^2 < \infty\), \[ \begin{array}{lll} \sigma^2 \sim \text{Gamma Invertida}(\alpha,\beta) \text{ então}\\ f(\sigma^2)=\frac{\beta^\alpha}{\Gamma (\alpha)} \left(\sigma^2\right)^{-\alpha+1} \exp\left[-\frac{\beta}{\sigma^2}\right] \end{array} \] \(\text{ com } \alpha>0 \text{ e } \beta>0\) parâmetros conhecidos.

Passo II: atribuir distribuição Normal para os dados \[ Y_1, Y_2,Y_3,\ldots,Y_n \sim \text{Normal}(\mu,\sigma^2) \text{ então }\\ \begin{array}{lll} f_{\boldsymbol{Y}}(\boldsymbol{y})&=&\prod_{i=1}^n\frac{1}{\sqrt{2 \pi} \sigma} \exp \left[-\frac{1}{2 \sigma^2} (y_i-\mu)^2 \right]\\ &=& \left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n \exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2\right]\\ &=&L(\sigma^2 \vert \boldsymbol{y}) \end{array}\\ \text{distribuição dos dados = função de verossimilhança de } \sigma^2 \text{ dado os dados}\\ \] - A função de verossimilhança nos traz toda a informação disponível na amostra (nos dados); - com \(\mu\) conhecido e \(\sigma^2\) desconhecido;

Passo III: Aplicar a fórmula - pode-se dispensar os termos constantes \[ \begin{array}{llll} f(\sigma^2 \vert \boldsymbol{y}) &\propto& f(\sigma^2) L(\sigma^2 \vert \boldsymbol{y})\\ &\propto& \frac{\beta^\alpha}{\Gamma (\alpha)} \left(\sigma^2\right)^{-\alpha+1} \exp\left[-\frac{\beta}{\sigma^2}\right] \times \left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n \exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2\right]\\ &\propto& \left(\sigma^2 \right)^{-\left(\alpha+\frac{n}{2}\right)+1} \exp\left[-\frac{\beta+\frac{1}{2}\sum_{i=1}^n (y_i-\mu)^2}{\sigma^2}\right]\\ & \text{Então}& \sigma^2 \vert \boldsymbol{y} \sim \text{Gamma Invertida} \left(\alpha + \frac{n}{2},\beta+\frac{1}{2} \sum \limits_{i=1}^n (y_i-\mu)^2\right), \text{ é só incluir as constantes na fórmula} \\ f(\sigma^2 \vert \boldsymbol{y}) &=& \frac{\left(\beta+\frac{1}{2} \sum_{i=1}^n (y_i-\mu)^2\right)^{\alpha + \frac{n}{2}}}{\Gamma[\alpha + \frac{n}{2}]} \left(\sigma^2 \right)^{-\left(\alpha+\frac{n}{2}\right)+1} \exp\left[-\frac{\beta+\frac{1}{2}\sum_{i=1}^n (y_i-\mu)^2}{\sigma^2}\right] \end{array} \]

Caso 5: Quando os dados têm distribuição normal com média conhecida e variância desconhecida

  • Atribui-se priori Gamma para a precisão \(\tau=\frac{1}{\sigma^2}\);
  • E a distribuição a posteriori para esta precisão será Gamma.

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 contínuo (pois \(0< \tau <\infty\));
  • dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(\tau\); \[ f(\tau \vert \boldsymbol{y})\propto f(\tau)L(\tau \vert \boldsymbol{y}) \]
  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Normal), ou seja, temos uma amostra de tamanho \(n\) i.i.d. (independentes e identicamente distribuídos) de uma distribuição de Normal com média igual a \(\mu\) e precisão igual a \(\tau\), então \(\boldsymbol{y}\) é um vetor de tamanho \(n\).

Passo I: atribuir priori para \(\tau\), sendo \(0<\tau < \infty\), \[ \begin{array}{lll} \tau \sim \text{Gamma}(\lambda,r) \text{ então}\\ f(\tau)=\frac{\lambda^r}{\Gamma (r)} \tau^{r-1} e^{-\lambda \tau} I_{(0,\infty)}(\tau) \end{array} \] \(\text{ com } \lambda>0 \text{ e } r>0\) parâmetros conhecidos (chamados de hiperparâmetros no contexto Bayesiano)

Passo II: atribuir distribuição Normal para os dados \[ Y_1, Y_2,Y_3,\ldots,Y_n \sim \text{Normal}(\mu,\sigma^2) \text{ então }\\ \begin{array}{lll} f_{\boldsymbol{Y}}(\boldsymbol{y})&=&\prod_{i=1}^n\frac{\sqrt{\tau}}{\sqrt{2 \pi}} \exp \left[-\frac{\tau}{2} (y_i-\mu)^2 \right]\\ &=& \left(\frac{\tau}{2\pi}\right)^{\frac{n}{2}} \exp\left[-\frac{\tau}{2}\sum_{i=1}^n(y_i-\mu)^2\right]\\ &=&L(\tau \vert \boldsymbol{y}) \end{array}\\ \text{distribuição dos dados = função de verossimilhança de } \tau \text{ dado os dados}\\ \] - A função de verossimilhança nos traz toda a informação disponível na amostra (nos dados); - com \(\mu\) conhecido e \(\tau\) desconhecido;

Passo III: Aplicar a fórmula - pode-se dispensar os termos constantes

\[ \begin{array}{lll} f(\tau \vert \boldsymbol{y}) &\propto& f(\tau) L(\tau \vert \boldsymbol{y})\\ &\propto& \frac{\lambda^r}{\Gamma (r)} \tau^{r-1} e^{-\lambda \tau} I_{(0,\infty)}(\tau) \times \left(\frac{\tau}{2\pi}\right)^{\frac{n}{2}} \exp\left[-\frac{\tau}{2}\sum_{i=1}^n(y_i-\mu)^2\right] \\ &\propto& \tau^{r+\frac{n}{2}-1} \exp \left[-\left(\lambda+\sum_{i=1}^n(y_i-\mu)^2\right) \tau \right] I_{(0,\infty)}(\tau)\\ & \text{Então}& \tau \vert \boldsymbol{y} \sim \text{Gamma} \left(\lambda+\frac{1}{2} \sum \limits_{i=1}^n (y_i-\mu)^2,r+\frac{n}{2}\right), \text{ é só incluir as constantes na fórmula} \\ f(\tau \vert \boldsymbol{y}) &=& \frac{\left(\lambda+\frac{1}{2} \sum \limits_{i=1}^n (y_i-\mu)^2 \right)^{r+\frac{n}{2}}}{\Gamma[r+\frac{n}{2}]}\tau^{r+\frac{n}{2}-1} \exp \left[-\left(\lambda+\sum_{i=1}^n(y_i-\mu)^2\right) \tau \right] I_{(0,\infty)}(\tau) \end{array} \]

  • E se dizemos que a distribuição a posteriori para \(\tau\) é gamma:
  • É equivalente a dizer que a distribuição a posteriori para \(\sigma^2\) é gamma invertida:

Exemplo 2 Abaixo temos \(10\) valores provenientes de uma distribuição normal com média \(\mu\) e variância \(\sigma^2\).

set.seed(05062017) #cria uma semente única
n=10
mu=2 #este é o valor da média mu
sigma2=3 #este é o valor da variancia sigma2 para a criação dos dados
y=rnorm(n,mean=mu,sd=sqrt(sigma2))
y=round(y,4)
y
##  [1]  1.2156  1.2000  2.1362  2.1139  2.6546  0.0135 -0.0007  0.2131  3.3849
## [10]  4.9196
    1. Obtenha a estimativa de máxima verossimilhança para \(\sigma^2\);
    1. Considere a média conhecida e igual a \(2\), e a variância desconhecida. Considere a distribuição a priori Gamma com média \(0.5\) e variância \(0.5\) para a precisão \(\tau=\frac{1}{\sigma^2}\), obtenha a distribuição a posteriori para \(\tau\);
    1. Segundo o item b), qual é a média a posteriori para \(\tau\)?
    1. Segundo o item b), qual é a varância a posteriori para \(\tau\)?

Solução: Você vai precisar de

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

    1. A estimativa de máxima verossimilhança para \(\sigma^2\) é \(\widehat{\sigma^2}=2.2837\), pois \[ Y_1,Y_2,\ldots,Y_n \text{ variáveis aleatórias i.i.d.} \sim N(\mu,\sigma^2), \text{ então}\\ f_{Y_i}(y)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{1}{2\sigma^2} (y_i-\mu)^2\right\} \text{ e aplicando o produtório}\\ L(\mu,\sigma^2 \vert \boldsymbol{y})= \left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n \exp\left\{-\frac{1}{2\sigma^2} \sum_{i=1}^n(y_i-\mu)^2\right\}\\ \frac{\partial L(\mu,\sigma^2 \vert \boldsymbol{y})}{\partial \sigma^2}=0 \Rightarrow \widehat{\sigma^2}=\frac{\sum_{i=1}^n \left(y_i-\bar{y}\right)^2}{n}, \text{ ou, equivalentemente, }\\ \widehat{\sigma^2}=\frac{s^2 (n-1)}{n}, \text{onde }s^2=\frac{\sum_{i=1}^n \left(y_i-\bar{y}\right)^2}{n-1} \text{ é a variância amostral}, \] e também pode ser verificado que a segunda derivada avaliada na e.m.v. tem sinal negativo.

  • Passo a passo na página 27 ;

  • Fazendo os cálculos no R:

s2=var(y)
sigma2_hat=s2*(n-1)/n
print(paste0("sigma2_hat= ",sigma2_hat))
## [1] "sigma2_hat= 2.2837365641"
  • Outro modo:
sigma2_hat=sum((y-mean(y))^2)/n
print(paste0("Modo alternativo: sigma2_hat= ",sigma2_hat))
## [1] "Modo alternativo: sigma2_hat= 2.2837365641"
  • Visualizando a log-verossimilhança:
logL=function(sigma2){
L=(1/sqrt(2*pi*sigma2))^n*exp(-1/(2*sigma2)*sum((y-mean(y))^2))
log(L)
}
optimize(logL, interval=c(0.01,10), maximum=TRUE) #encontra o ponto de máximo
## $maximum
## [1] 2.283739
## 
## $objective
## [1] -18.31845
sigma2=seq(0.01,20,0.01) #só assume valores positivos
temp <- data.frame(sigma2=sigma2, logL = logL(sigma2))
datatable(temp)
ggplot(data = temp, aes(x = sigma2, y = logL)) + geom_line() + stat_peaks(col = "red")
## Warning: Removed 1 rows containing non-finite values (stat_peaks).

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

\[ \tau \vert \boldsymbol{y} \sim \text{Gamma}\left(\lambda+\frac{1}{2}\sum_{i=1}^n (y_i-\mu)^2,r+\frac{n}{2}\right), \text{ ou seja, }\\ \text{Como } \lambda=1, r=\frac{1}{2} \text{ e } \mu=2:\\ \tau \vert \boldsymbol{y} \sim \text{Gamma}\left(\frac{2+\sum_{i=1}^n (y_i-2)^2}{2},\frac{1+n}{2}\right) \] - Fazendo os cálculos no R: \(\tau | \boldsymbol{y} \sim \text{Gamma} (12.6496,5.5)\)

par1=(2+sum((y-2)^2))/2
par2=(1+n)/2
print(paste0("par_1= ",par1," e par_2= ",par2))
## [1] "par_1= 12.649657345 e par_2= 5.5"
    1. e d) Média e variância a posteriori: \[ \begin{array}{lll} \text{E}(\tau| \boldsymbol{y})&=&\frac{5.5}{12.6497} \approx 0.4351, \text{ próxima da e.m.v: } \widehat{\tau}= \frac{1}{2.2837}=0.4379,\\ \text{VAR}(\tau| \boldsymbol{y})&=&\frac{5.5}{12.6497^2}=0.0344 \end{array} \]
  • Graficamente:
tau=seq(0,1.5,0.01)     
lambda=1
r=1/2
mu=2
constante=exp(19)

priori=dgamma(tau,shape=r, scale=1/lambda) 
L_tau=(1/sqrt(2*pi))^n*tau^(n/2)*exp(-tau/2*sum((y-mu)^2))

#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
L_tau_normalizado=L_tau/sum(L_tau*intervalos)   
posteriori=dgamma(tau,shape=par2, scale=1/par1)

plot(tau,priori,type='l',xlab=expression(tau),ylab=expression(f(tau)))
lines(tau,L_tau_normalizado,type='l',col=2)  
lines(tau,posteriori,type='l',col=3)
legend("topright",c(expression("priori "*f(tau)),expression("verossimilhança "*L(tau*"|"*bold(y))),expression("posteriori "*f(tau*"|"*bold(y)))),col=c(1,2,3),bty = "n",lty=c(1,1,1))

Figura 2: Gráficos do exemplo 2

Caso 6: Vemos que a família de distribuições Beta é conjugada a distribuição Geométrica:

Sejam \(Y_1,Y_2,\ldots,Y_n\) variáveis aleatórias i.i.d. com distribuição geométrica com probabilidade de sucesso igual a \(p\).

  • Distribuição geométrica (segundo apêndice do Mood): Cada \(Y_i\) é igual ao número de tentativas anteriores ao primeiro sucesso em um experimento com ensaios independentes de Bernoulli com probabilidade de sucesso igual a \(p\).
  • Atribui-se priori Beta para a proporção \(p\);
  • E a distribuição a posteriori para esta proporção será Beta.

Resolução:

  • Pela fórmula de Bayes, dispensando o termo do denominador (pois é uma constante com relação ao parâmetro \(p\); \[ f(p \vert \boldsymbol{y})\propto f(p)L(p \vert \boldsymbol{y}) \]

  • o termo \(\boldsymbol{y}\) se refere aos dados observados (Modelo Geométrico), ou seja, quando realizados ensaios de Bernoulli independentes com probabilidade igual a \(p\), verificou-se os valores observados \(y_1,y_1,\ldots,y_n\) das tentativas anteriores ao primeiro sucesso.

Passo I: atribuir priori para \(p\), sendo \(0<p<1\), \[ \begin{array}{lll} p \sim \text{Beta}(a,b) \text{ então }\\ f(p)=\frac{1}{B(a,b)} \times p^{a-1} (1-p)^{b-1} I_{(0,1)}(p) \end{array} \]

Passo II: atribuir distribuição geométrica para os dados \[ \begin{array}{lll} Y_1, Y_2,\ldots,Y_n \sim \text{Geom}(p) \text{ então }\\ P(\boldsymbol{Y}=\boldsymbol{y})&=&\underbrace{\prod_{i=1}^n p (1-p)^{y_i}}_{\text{os experimentos são independentes}} \\ &=& p^n (1-p)^{\sum_{i=1}^n y_i}\\ &=& L(p \vert \boldsymbol{y}) \end{array}\\ \text{distribuição dos dados = função de verossimilhança de }p\text{ dado os dados} \] Passo III: Aplicar a fórmula - pode-se dispensar os termos constantes \[ \begin{array}{lll} f(p \vert \boldsymbol{y})&\propto& f(p)L(p \vert \boldsymbol{y})\\ &\propto& \frac{1}{B(a,b)} \times p^{a-1} (1-p)^{b-1} I_{(0,1)}(p) \times p^n (1-p)^{\sum_{i=1}^n y_i}\\ &\propto& p^{a+n-1} \times (1-p)^{b+\sum_{i=1}^n y_i-1} \times I_{(0,1)}(p)\\ &\text{ Então }& p \vert \boldsymbol{y} \sim \text{Beta}(a+n,b+\sum_{i=1}^n y_i) \text{ é só incluir as constantes na fórmula:}\\ f(p \vert \boldsymbol{y})&=& \frac{1}{B\left(a+n,b+\sum_{i=1}^n y_i\right)}\times p^{a+n} \times (1-p)^{b+\sum_{i=1}^n y_i-1} \times I_{(0,1)}(p) \end{array}\\ \] Implementação numérica do Caso 6:

set.seed(10102020) #cria uma semente única
p0=0.3 #este é o valor de p utilizado para gerar os dados
n=10 #tamanho da amostra
y=rgeom(n,prob=p0)
y
##  [1] 4 6 0 1 0 1 7 7 0 0
logL=function(p){
L=p^n* (1-p)^sum(y)
log(L)
}
optimize(logL, interval=c(0.001, 0.999), maximum=TRUE) #encontra o ponto de máximo
## $maximum
## [1] 0.2777811
## 
## $objective
## [1] -21.27032
p=seq(0.001,0.999,0.001)
temp <- data.frame(p, logL = logL(p))
datatable(temp)
ggplot(data = temp, aes(x = p, y = logL)) + geom_line() + stat_peaks(col = "red")

a=0.4
b=0.6
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.001
vero=p^n* (1-p)^sum(y)
vero_normalizado=vero/sum(vero*intervalos)
pos=dbeta(p,shape1=a+n,shape2=b+sum(y))
plot(p,pri,type='l',xlab=expression(p),ylab=expression(f(p)))
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))

Cenários Distribuição a priori Distribuição dos dados Distribuição a posteriori
Caso 1 \(\mu \sim\) N\((m_0,\sigma_0^2\)) \(X_1,\ldots,X_n \sim\) N\((\mu,\sigma^2)\), com \(\sigma^2\) conhecido \(\mu \vert \boldsymbol{x} \sim\) N\(\left(\frac{\frac{m_0}{s_0^2}+\frac{n \overline{y}}{\sigma^2}}{\frac{1}{s_0^2}+\frac{n}{\sigma^2}},\frac{1}{\frac{1}{s_0^2}+\frac{n}{\sigma^2}}\right)\)
Caso 2 \(\lambda \sim\) Gamma \((a,r)\) \(X_1,\ldots,X_n \sim\) Poisson \((\lambda)\) \(\lambda \vert \boldsymbol{x} \sim\) Gamma\((a+n,r+\sum_{i=1}^n x_i)\)
Caso 3 \(p \sim\) Beta\((a,b)\) \(X \sim\) Binomial\((n,p)\) \(p \vert \boldsymbol{x} \sim\) Beta\((a+y,b+n-x)\)
Caso 4 \(\sigma^2 \sim\) Gamma Inv. \((\alpha,\beta)\) \(X_1,\ldots,X_n \sim\) N\((\mu,\sigma^2)\), com \(\mu\) conhecido \(\sigma^2 \vert \boldsymbol{x} \sim\) Gamma Inv.\((\alpha+\frac{n}{2},\beta+\frac{n}{2} \sum_{i=1}^n (x_i-\mu)^2)\)
Caso 5 \(\tau \sim\) Gamma\((\lambda,r)\) \(X_1,\ldots,X_n \sim\) N\((\mu,\frac{1}{\tau})\), com \(\mu\) conhecido \(\tau \vert \boldsymbol{x} \sim\) Gamma\((\lambda+\frac{1}{2}\sum_{i=1}^n (x_i-\mu)^2,r+\frac{n}{2})\)
Caso 6 \(p \sim\) Beta\((a,b)\) \(X_1,\ldots,X_n \sim\) Geom\((p)\) \(p \vert \boldsymbol{x} \sim\) Beta\((a+n,b+\sum_{i=1}^n x_i)\)

0.3 Exercícios

    1. Mostre que a família de distribuições Beta é conjugada em relação à binomial, geométrica e binomial negativa.
    1. Para uma amostra aleatória \(X_1,\dots,X_n\) tomada da distribuição \(U(0,\theta)\), mostre que a família de distribuições de Pareto com parâmetros \(a\) e \(b\), cuja função de densidade é \(f(\theta)=\frac{ab^a}{\theta^{a+1}}\), é conjugada à uniforme.
    1. Suponha que o tempo, em minutos, para atendimento a clientes segue uma distribuição exponencial com parâmetro \(\theta\) desconhecido. Com base na experiência anterior assume-se uma distribuição a priori Gamma com média \(0.2\) e desvio-padrão \(1\) para \(\theta\). Se o tempo médio para atender uma amostra aleatória de \(20\) clientes foi de \(3.8\) minutos, determine a distribuição a posteriori de \(\theta\).
    1. Seja \(X_1,\dots,X_n\) uma amostra aleatória da distribuição de Poisson com parâmetro \(\theta\). Determine os parâmetros da priori conjugada de \(\theta\) sabendo que E\((\theta)=4\) e o coeficiente de variação a priori é igual a \(0.5\).
    1. O número médio de defeitos por \(100\) metros de uma fita magnética é desconhecido e denotado por \(\theta\). Atribui-se uma distribuição a priori \(\text{Gamma }(2,10)\) para \(\theta\). Se um rolo de \(1200\) metros desta fita foi inspecionado e encontrou-se \(4\) defeitos, qual é a distribuição a posteriori de \(\theta\)?