Teoria das probabilidades e axiomas

Definição 1.1: Partição de um Espaço Amostral Dizemos que os eventos \(A_1,A_2,\ldots,A_n\) formam uma partição do espaço amostral \(\Omega\) se as seguintes propriedades são satisfeitas:

• \(A_i\neq\emptyset,i=1,\ldots,n\): significa que nenhum evento pode ser igual ao conjunto vazio;
• \(A_i\) \(\cap\) \(A_j\) = \(\emptyset\), para \(i\) \(\neq\) \(j\): significa que os eventos são disjuntos;
• \(\cup^{n}_{i=1}A_i=\Omega\): significa que a união (ou reunião) de todos os eventos totaliza o espaço amostral.

Figura 1: Representação Gráfica de partição de um espaço amostral

Definição 1.2: Classe de eventos do espaço amostral A classe de eventos do espaço amostral \(\Omega\), também chamada de classe de subconjuntos do espaço amostral \(\Omega\), ou Conjunto das partes de \(\Omega\), é o conjunto que contém todos os subconjuntos de \(\Omega\) e é representado por \(\mathcal{P}(\Omega)\).

Conceito de probabilidade A probabilidade é definida numa classe de eventos do espaço amostral, com certas propriedades.

Definição 1.3: Probabilidade é uma função \(P(.)\) que associa a cada evento de \(\mathcal{P}(\Omega)\) (ou subconjunto de \(\Omega\)) um número real pertencente ao intervalo \([0,1]\), satisfazendo aos Axiomas de Kolmogorov

• Axioma 1: \(P(A) \geq 0\) para todo evento \(A\), \(A \subset \Omega\): significa que a probabilidade é sempre um número real não negativo;
• Axioma 2: \(P(\Omega)=1\): significa que \(\Omega\) é um evento certo pois reúne todas as possibilidades;
• Axioma 3: \(P(\cup^{\infty}_{i=1}A_i)=\sum_{i=1}^{\infty}P(A_{i})\), se \(A_1, A_2, \ldots\) forem, dois a dois, mutuamente exclusivos:

significa que a probabilidade da união de dois ou mais eventos é igual à soma de suas respectivas probabilidades, se os eventos forem mutuamente exclusivos aos pares.

Teorema 1.4: Se os eventos \(A_1, A_2, \ldots, A_n\) formam uma partição do espaço amostral, então: \[ \sum_{i=1}^n P(A_{i})=1. \]

Demonstração: A demonstração vem da definição de partição e dos Axiomas 2 e 3 de Kolmogorov.

Definição 1.5: Probabilidade condicional A probabilidade condicional de \(B\) dado \(A\) é dada pela fórmula: \[ P(B|A) = \frac{P(A \cap B)}{P(A)}, \] sendo que \(P(A)\) deve ser maior do que zero.

Teorema 1.6: Teorema do produto \[ P(A \cap B) = P(A) P(B|A) \text{ e também } P(A \cap B) = P(B) P(A|B). \]

Demonstração: A demonstração vem da definição de probabilidade condicional.

Proposição 1.7: Generalização do Teorema do Produto Sejam \(A_1\), \(A_2\),\(\ldots\), \(A_{n-1}\), \(A_n\) eventos do espaço amostral \(\Omega\), onde está definida a probabilidade \(P\), temos: \[ P(A_1\cap A_2 \ldots \cap A_{n-1}\cap A_n) = P(A_1)P(A_2|A_1)P(A_3|A_1\cap A_2)P(A_n|A1\cap A2 \ldots A_{n-1}). \] Demonstração: A demonstração é através do Princípio da Indução Finita.

Teorema 1.8: Teorema da Probabilidade Total (ou Fórmula da Probabilidade Total): Sejam \(A_1,A_2,\ldots,A_n\) eventos que formam uma partição do espaço amostral. Seja \(B\) um evento deste espaço. \[ P(B) = \sum_{i=1}^n P(A_i)P(B|A_i), \] onde \(A_1, A_2, \ldots A_n\) formam uma partição no espaço amostral.

A fórmula da Probabilidade Total permite calcular a probabilidade de um evento \(B\) a partir das probabilidades de um conjunto de eventos disjuntos cuja reunião é o espaço amostral; e as probabilidades condicionais de \(B\) dado cada um destes eventos são fornecidas.

Demonstração - Passo 1: Como \(A_1, \ldots, A_n\) formam uma partição, então podemos escrever \(B\) da seguinte forma: \[ B = (B \cap A_1) \cup (B \cap A_2) \cup (B \cap A_3) \cup \ldots \cup (B \cap A_n). \]

Figura 2: Representação gráfica do teorema da probabilidade total

• Passo 2: Como \(A_1,A_2,\ldots,A_n\) são disjuntos, então pelo axioma 3 de Kolmogorov, temos: \[ P(B)= P(B \cap A_1) + P(B \cap A_2) + P(B \cap A_3) \ldots + P(B \cap A_n). \]
• Passo 3: Utilizando o Teorema do Produto, podemos escrever: \[ P(B)=P(A_1)P(B|A_1)+ P(A_2)P(B|A_2)+P(A_3)P(B|A_3)+ \ldots +P(A_n)P(B|A_n). \]

Teorema 1.9: Teorema de Bayes (ou fórmula de Bayes):

\[ P(A|B)=\frac{P(B|A)P(A)}{P(B)}. \] Demonstração: A demonstração vem da Definição de probabilidade condicional, Teorema do Produto e Teorema da Probabilidade Total.

Teorema de Bayes - Caso geral \[ P(A_i|B)=\frac{P(A_i)P(B|A_i)}{\sum_{i=1}^n P(A_i)P(B|A_i)}, \] onde \(A_1,A_2,\ldots A_n\) formam uma partição no espaço amostral

Exemplo 1: Exames de diagnóstico não são infalíveis, mas deseja-se que tenha probabilidade pequena de erro. Um exame detecta uma certa doença, caso ela exista, com probabilidade 0,9. Se a doença não existir, o exame corretamente aponta isso com probabilidade 0,8. Considere que estamos aplicando esses exames em uma população com \(10\%\) de prevalência dessa doença. Para um indivíduo escolhido ao acaso, pergunta-se:

• 1. A probabilidade de ser realmente doente se o exame indicou que era.
• 2. A probabilidade de acerto do exame.

Resolução:

• Sejam os eventos A: o indivíduo ter a doença e B: o teste dar positivo.
• Podemos construir o diagrama da árvore

Figura 3: Diagrama da árvore

• item a) esta probabilidade é denotada por \(P(doente \vert +)= P(A|B)\) não é mesma coisa que \(P(B \vert A)\).

Passo I: Calcular a probabilidade Total - que vai no denominador da probabilidade que queremos:

\[ \begin{array}{lll} P(B)&=&P(B|A) \times P(A)+P(B|A^c) \times P(A^c)\\ &=&0,9 \times 0,1+0,2 \times 0,9 \\ &=&0,27 \end{array} \]

Passo II: Calcular a probabilidade condicional \[ \begin{array}{lll} P(A|B)&=& \frac{P(A \cap B)}{P(B)}\\ &=& \frac{P(B| A) \times P(A)}{P(B)} \text{ , o numerador vem da fórmula do produto} \\ &=& \frac{0,9 \times 0,1}{0,27} = \approx 0,33 = 33\% \end{array} \] - item b) esta probabilidade é dada pela soma: P(+ e o indivíduo ser doente ) + P(- e o indivíduo não ser doente): \[ \begin{array}{lll} P(A \cap B) + P(A^c \cap B^c) &=& P(B|A) \times P(A) + P( B^c|A^c) \times P(A^c) \text{ pela fórmula do produto} \\ &=& 0,9 \times 0,1 + 0,8 \times 0,9 = 0,81 = 81 \% \end{array} \]

Exercícios

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,
• 2. 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

Exercícios - continuação

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

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

• 2. 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")

• 2. 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 4: Priori, Verossimilhança e Posteriori

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 = 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 = 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)

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

• 2. 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")

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.

Casos principais de prioris conjugadas

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 5: Gráficos das três funções conjuntamente: priori, verossimilhança e posteriori

• 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.

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

• 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\);
• 2. 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\);
• 3. Segundo o item b), qual é a média a posteriori para \(\tau\)?
• 4. 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).

• 2. 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"

• 3. 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 6: Gráficos do exemplo 2

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

require(DT)
require(tidyverse)
require(ggpmisc)

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))

Síntese de prioris conjugadas

Exercícios

Princípio da Verossimilhança

Exemplo: Suponha que desejamos estimar \(\theta\), a probabilidade de observar cara (C) no lançamento de uma moeda e que, para um determinado experimento, observou-se: \[ \{C,\bar{C},\bar{C},C,C,\bar{C},\bar{C},\bar{C},\bar{C},C\} \] Entre outras possibilidades, os dados acima podem ter sido gerados a partir dos seguintes experimentos:

• Seja \(X\) o número de caras em \(10\) lançamentos da moeda: \[ X \sim \text{ Binomial }(10,\theta), \text{ e os resultados poderiam ser: } X = 0,1,2,\ldots,10. \]

• Seja \(Y\) o número de lançamentos da moeda até a obtenção de \(4\) caras: \[ Y \sim \text{ Binomial Negativa }(4,\theta), \text{ e os resultados poderiam ser: } Y = 4,5,6,\ldots \]

• Considerando os resultados do experimento no modelo Binomial:

\[ \begin{array}{lll} &&P(X = 4 \mid \theta) = \left(\begin{array}{l l}10\\4\end{array}\right) \theta^4 (1 - \theta)^{10-4}, \\ \end{array} \]

de modo que a função de verossimilhança será: \(L(\theta \mid \text{ x} ) \propto \theta^4 (1 - \theta)^6\);

• E no modelo Binomial Negativa:

\[ \begin{array}{lll} && P(Y = 10 \mid \theta) = \left(\begin{array}{l l}10-1\\4-1\end{array}\right) \theta^4 (1 -\theta)^{10-4},\\ && \text{ de modo que a função de verossimilhança será: } L(\theta \mid \text{ y} ) \propto \theta^4 (1 - \theta)^6; \end{array} \]

• \(X\) Sob a mesma priori para \(\theta\), a posteriori obtida a partir do modelo Binomial é igual à posteriori obtida a partir do modelo Binomial-Negativa.

• Porém, as estimativas de máxima verossimilhança sob cada um dos modelos são diferentes. Tarefa: justificar esta afirmação;

• Formalmente: Se temos dois vetores aleatórios pertencentes a um mesmo espaço amostral, que dependem do mesmo parâmetro \(\theta\) e que possuem verossimilhancas distintas, diferindo apenas por uma constante que não depende de \(\theta\), Então as posterioris obtidas a partir destes dois vetores são iguais.

• Em outras palavras, a inferência Bayesiana é a mesma quando a condição de proporcionalidade das verossimilhanças é satisfeita.

prioris não informativas

• As prioris não informativas estão presentes quando se espera que a informação dos dados seja dominante, significa que a informação a priori é vaga, então temos o conceito de “conhecimento vago”, “não informação” ou “ignorância a priori”.
• Referências sobre prioris não informativas estão em , e .

priori Uniforme

É uma priori intuitiva porque todos os possíveis valores do parâmetro \(\theta\) são igualmente prováveis: \[ f(\theta)\propto k, \] com \(\theta\) variando em um subconjunto da reta de modo que nenhum valor particular tem preferência (Bayes, 1763).

A priori uniforme, no entanto, apresenta algumas dificuldades:

• Se o intervalo de variação de \(\theta\) for a reta real então a distribuição é imprópria: \[ \int \limits_{-\infty}^{\infty} f(\theta )d\theta =\infty, \] mas este não chega a ser um impedimento para a escolha de prioris, como veremos mais adiante.
• Se \(\phi =g(\theta )\) é uma reparametrização não linear monótona de \(\theta\) então a priori para o parâmetro \(\phi\) será: \[ f(\phi)=f(\theta(\phi))\left\vert\frac{d\theta}{d\phi}\right\vert\propto \left\vert\frac{d\theta}{d\phi}\right\vert, \] e vemos pelo Teorema de transformação de variáveis que a priori para \(\phi\) não é uniforme.

priori de Jeffreys

É uma priori construída a partir da medida de informação esperada de Fisher, proposta por Jeffreys (1961).

• é uma priori imprópria;
• é invariante a transformações \(1\) a \(1\).

Definição: Medida de informação esperada de Fisher Considere uma única observação \(X\) com f.d.p. indexada pelo parâmetro \(\theta\): \(f(x \vert \theta)\). A medida de informação esperada de Fisher de \(\theta\) através de \(X\) é definida como \[ \displaystyle I(\theta)=E \left[-\frac{\partial^2\log f(x \vert \theta)}{\partial\theta^2}\right], \] em que a esperança matemática é tomada em relação à distribuição amostral \(f(x \vert \theta)\) (a esperança é com respeito a \(X\) e não com respeito a \(\theta\)). - A informação esperada de Fisher \(I(\theta)\) é uma medida de informação global.

Extendendo esta definição para uma amostra i.i.d. \(X_1,X_2,\ldots,X_n\), temos: \(f(\boldsymbol{x} \vert \theta)=\prod \limits_{i=1}^n f(x_i \vert \theta)\) e \[ \displaystyle I(\theta)=E \left[-\frac{\partial^2\log f(\boldsymbol{x} \vert \theta)}{\partial\theta^2}\right], \] é a informação esperada de Fisher de \(\theta\) através do vetor \(\boldsymbol{x}\).

Definição: priori de Jeffreys A priori de Jeffreys é dada por: \[ \sqrt{I(\theta)}. \]

No caso multiparamétrico (mais de um parâmetro), a medida de Informação de Fisher é dada de forma matricial, então temos: \[ \sqrt{\left| \text{det}\left[I(\boldsymbol{\theta})\right]\right|}. \]

Exemplo: Sejam \(X_1,\ldots,X_n \sim \text{Poisson} (\theta)\).

\[ \begin{array}{lll} \log f(\boldsymbol{x} \vert \theta) &=&-n\theta + \sum_{i=1}^n x_i \log (\theta) - \log\left(\prod_{i=1}^n x_i!\right) \\ \frac{\partial \log f(\boldsymbol{x} \vert\theta)}{\partial\theta}&=&-n+\frac{\sum_{i=1}^n x_i}{\theta}\\ \frac{\partial^2 \log f(\boldsymbol{x} \vert \theta)}{\partial\theta^2}&=&-\frac{\sum_{i=1}^n x_i}{\theta^2}\\ I(\theta) &=& \frac{n}{\theta}\\ &\propto& \frac{1}{\theta}\\ \end{array} \]

• A priori de Jeffreys para \(\theta\) no modelo Poisson é \(f(\theta)\propto\theta^{-1/2}\);
• Esta priori também pode ser obtida tomando-se a priori conjugada \(\text{Gamma}(\alpha,\beta)\) com \(\alpha=\frac{1}{2}\) e \(\beta \rightarrow 0\). Note que o parâmetro \(\beta\) é sempre positivo, por isso a noção de “tender a zero”. Tarefa: verificar;
• Em geral, quando o modelo admite priori conjugada, basta fixar um dos parâmetros da priori conjugada e o outro parâmetro “tender a zero”, resultando na priori de Jeffreys;
• A priori de Jeffreys não satisfaz o princípio da verossimilhança, pois a informação esperada de Fisher depende da distribuição amostral (o cálculo das esperanças matemáticas podem ser diferentes se os modelos forem diferentes como no exemplo modelos Binomial e Binomial-Negativa).
• A priori de Jeffreys apresenta algumas particularidades nos modelos de locação-escala, como veremos a seguir.

Modelos de locação-escala

Modelo de Locação \(X\) tem um modelo de locação se existem uma função \(g\) e uma quantidade \(\mu\) tais que: \[ f(x \vert \mu)=g(x-\mu), \] logo \(\mu\) é o parâmetro de locação.

• A definição se extende para o caso multiparamétrico;
• Exemplo: distribuição normal com variância conhecida é um modelo de locação:

mu=c(-1,0,1)  #cria um vetor de médias de -1, 0 e 1
sigma=2 #fixa o desvio padrao em 2, ou seja, a variancia é igual a 4
x=seq(-10,10,0.1)
y_1=dnorm(x,mean=mu[1],sd=sigma)
y_2=dnorm(x,mean=mu[2],sd=sigma)
y_3=dnorm(x,mean=mu[3],sd=sigma)
dados=data.frame(x,y_1,y_2,y_3)

colors <- c("y_1"="blue", "y_2"="red", "y_3"="green")

ggplot(dados, aes(x = x,y = y_1, color = "y_1")) +
    geom_line() +
    geom_line(aes(x = x,y = y_2, color = "y_2")) +
    geom_line(aes(x = x,y = y_3, color = "y_3")) +
    labs(x="x",
         y="f(x)",
         title=expression(N(mu,4)),
         color = "Valores da média") +
    scale_color_manual(labels=c(expression(mu==-1),expression(mu==0),expression(mu==1)),values = colors)

• È possivel calcular valores de uma normal não central a partir da normal centrada no zero

x=seq(4,6,.5)
f=  dnorm(x,mean=5,sd=2)
g=  dnorm(x-5,mean=0,sd=2)
temp=data.frame(x,f,g)

require(kableExtra)
names(temp) = c("$x$","$f(x)$","$g(x-5$)")
temp %>%
   kbl(escape = FALSE) %>%
 kable_classic(full_width = F, html_font = "Cambria") %>%
kable_styling(bootstrap_options = c("striped","hold_position")) %>%
 add_header_above(c("","$f(.)$ vem da N(5,4)","$g(.)$ vem da N(0,4)"))

• Propriedade: A priori de Jeffreys para o parâmetro de locação \(\mu\) é: \[ f(\mu)\propto k, \] onde \(k\) é uma constante.

Modelo de Escala \(X\) tem um modelo de escala se existem uma função \(g\) e uma quantidade \(\sigma\) tais que: \[ f(x \vert \sigma)=\frac{1}{\sigma} g\left(\frac{x}{\sigma}\right), \] logo \(\sigma\) é o parâmetro de escala.

• Exemplos: Na distribuição \(\text{Exp}(\theta)\) o parâmetro de escala é \(\sigma=\frac{1}{\theta}\), e na distribuição \(\text{N}(\mu,\sigma^2)\) com média conhecida o parâmetro de escala é \(\sigma\);

• Propriedade: A priori de Jeffreys para o parâmetro de escala \(\sigma\) é: \[ f(\sigma) \propto \frac{1}{\sigma}. \]

• Mostrando que o Modelo normal com media conhecida é modelo de escala

sigma=c(1,2,3)  #cria um vetor de desvios padroes de 1, 2 e 3
media=0 #fixa a média em zero
x=seq(-10,10,0.1)
y_1=dnorm(x,mean=media,sd=sigma[1])
y_2=dnorm(x,mean=media,sd=sigma[2])
y_3=dnorm(x,mean=media,sd=sigma[3])
dados=data.frame(x,y_1,y_2,y_3)

colors <- c("y_1"="blue", "y_2"="red", "y_3"="green")

ggplot(dados, aes(x = x,y = y_1, color = "y_1")) +
    geom_line() +
    geom_line(aes(x = x,y = y_2, color = "y_2")) +
    geom_line(aes(x = x,y = y_3, color = "y_3")) +
    labs(x="x",
         y="f(x)",
         title=expression(N(0,sigma^2)),
         color = "Valores da variância") +
    scale_color_manual(labels=c(expression(sigma^2==1),expression(sigma^2==4),expression(sigma^2==9)),values = colors)

• È possivel calcular valores de uma normal com desvio padrão diferente de 1 a partir de uma normal com desvio padrão igual a 1:

x=seq(4,6,.5)
mu=5
sigma=2
f=dnorm(x,mean=5,sd=sigma)
g=1/sigma*dnorm(x/sigma,mean=mu/sigma,sd=1)

temp=data.frame(x,f,g)
names(temp) = c("$x$","$f(x)$","$\\frac{1}{2}g(\\frac{x}{2})$")
temp %>%
   kbl(escape = FALSE) %>%
 kable_classic(full_width = F, html_font = "Cambria") %>%
kable_styling(bootstrap_options = c("striped","hold_position")) %>%
 add_header_above(c("","$f(.)$ vem da N(5,4)","$g(.)$ vem da N(5,1)"))

Definição: Modelo de Locação-escala \(X\) tem um modelo de locação-escala se existem uma função \(g\) e as quantidades \(\mu\) e \(\sigma\) tais que \[ f(x\vert\mu,\sigma)=\frac{1}{\sigma}g\left(\frac{x-\mu}{\sigma}\right), \] logo \(\mu\) é o parâmetro de locação e \(\sigma\) é o parâmetro de escala.

• Exemplos: Na distribuição \(\text{N}(\mu,\sigma^2)\) o parâmetro de locação é \(\mu\) e o parâmetro de escala é \(\sigma\), e a distribuição de Cauchy também é um modelo de locação-escala.

    #é modelo de locação-escala
media=c(-1,1) #a média assume os valores -1 ou 1
sigma=c(1,2)  #o desvio padrão assume os valores 1 ou 2
x=seq(-10,10,0.1)
y_1= dnorm(x,media[1],sigma[1])
y_2= dnorm(x,media[1],sigma[2])
y_3= dnorm(x,media[2],sigma[1])
y_4= dnorm(x,media[2],sigma[2])
dados=data.frame(x,y_1,y_2,y_3,y_4)
    
colors <- c("y_1"="blue", "y_2"="red", "y_3"="green","y_4"="orange")

ggplot(dados, aes(x = x,y = y_1, color = "y_1")) +
    geom_line() +
    geom_line(aes(x = x,y = y_2, color = "y_2")) +
    geom_line(aes(x = x,y = y_3, color = "y_3")) +
    geom_line(aes(x = x,y = y_4, color = "y_4")) +
    labs(x="x",
         y="f(x)",
         title=expression(N(mu,sigma^2)),
         color = "Valores da média e variância") +
    scale_color_manual(labels=c(
      expression(mu==-1~"e"~sigma^2==1),
      expression(mu==-1~"e"~sigma^2==4),
      expression(mu==1~"e"~sigma^2==1),
      expression(mu==1~"e"~sigma^2==4)),
      values = colors)

• È possivel calcular valores de uma normal genérica a partir de uma normal padrão:

x=seq(4,6,.5)
mu=5
sigma=2
f=dnorm(x,mean=5,sd=sigma)
g=1/sigma*dnorm((x-mu)/sigma,mean=0,sd=1)

temp=data.frame(x,f,g)
names(temp) = c("$x$","$f(x)$","$\\frac{1}{2}g(\\frac{x-5}{2})$")
temp %>%
   kbl(escape = FALSE) %>%
 kable_classic(full_width = F, html_font = "Cambria") %>%
kable_styling(bootstrap_options = c("striped","hold_position")) %>%
 add_header_above(c("","$f(.)$ vem da N(5,4)","$g(.)$ vem da N(0,1)"))

• Propriedade A priori conjunta de Jeffreys para os parâmetros de locação \(\mu\) e escala \(\sigma\) é: \[ f(\mu,\sigma)=f(\mu)f(\sigma) \propto \frac{1}{\sigma}, \] onde nós assumimos independência (a priori conjunta é o produto das prioris).

Exemplo: Sejam \(X_1,\dots,X_n \sim N(\mu,\sigma^2)\) com \(\mu\) e \(\sigma^2\) desconhecidos, temos: \[ f\left(x \vert \mu, \sigma^2 \right)=\frac{1}{\sigma} \left\{\frac{1}{\sqrt{2\pi}} \exp\left[-\frac{1}{2} \left(\frac{ x-\mu}{\sigma}\right)^2\right] \right\}, \] logo \(\mu\) é o parâmetro de locação e \(\sigma\) é o parâmetro de escala.

• A priori não informativa de Jeffreys para o vetor \((\mu,\sigma)\) é: \[ f(\mu,\sigma)\propto\frac{1}{\sigma} \]
• Pela propriedade da invariância, a priori não informativa de Jeffreys para o vetor \((\mu,\sigma^2)\) é: \[ f(\mu,\sigma^2)\propto\frac{1}{\sigma^2} \]

Síntese de prioris de Jeffreys

tabela <- data.frame(
      col0 =  c("Cenários","Caso 1","Caso 2","Caso 3","Caso 4","Caso 5","Caso 6"),
      col1 =  c("priori de Jeffreys (é não informativa)","$\\mu \\sim$ N$(m_0,\\sigma_0^2$)","$\\lambda \\sim$ Gamma $(a,r)$","$p \\sim$ Beta$(a,b)$","$\\sigma^2 \\sim$ Gamma Inv. $(\\alpha,\\beta)$","$\\tau \\sim$ Gamma$(\\lambda,r)$","$p \\sim$ Beta$(a,b)$"),
      col2 =  c("Distribuição dos dados","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\sigma^2)$, com $\\sigma^2$ conhecido", "$X_1,\\ldots,X_n \\sim$ Poisson $(\\lambda)$", "$X \\sim$ Binomial$(n,p)$","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\sigma^2)$, com $\\mu$ conhecido","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\frac{1}{\\tau})$, com $\\mu$ conhecido","$X_1,\\ldots,X_n \\sim$ Geom$(p)$"),
      col3 =  c("Distribuição dos dados","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\sigma^2)$, com $\\sigma^2$ conhecido", "$X_1,\\ldots,X_n \\sim$ Poisson $(\\lambda)$", "$X \\sim$ Binomial$(n,p)$","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\sigma^2)$, com $\\mu$ conhecido","$X_1,\\ldots,X_n \\sim$ N$(\\mu,\\frac{1}{\\tau})$, com $\\mu$ conhecido","$X_1,\\ldots,X_n \\sim$ Geom$(p)$"))

      
 tabela %>%
  kbl(col.names = NULL,escape = FALSE) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
     kable_styling(bootstrap_options = c("striped","hold_position"),
                full_width = F)

Exemplo 3.1: Regressão linear simples

O problema envolve as variáveis \(X\): a dose de um medicamento anti-alérgico em estudo, e \(Y\): tempo de duração do efeito (alívio dos sintomas alérgicos). Abaixo temos a representação gráfica dos dados observados Pelo gráfico, nós concluímos que uma relação linear (reta) é satisfatória para os dados. Também iremos supor que \(X\) é uma variável controlada pelo pesquisador (sem a presença de erros).

x=c(3,3,4,5,6,6,7,8,8,9)
y=c(9,5,12,9,14,16,22,18,24,22)
a=cbind(x,y)
tab_nums <- captioner(prefix = "Tabela")
tab_cap = tab_nums("tabela1","Dados do Exemplo")
datatable(as.data.frame(a),caption=tab_cap)

plot(x,y,lwd=2,xlim=c(0,10),ylim=c(0,25))

Figura 7: Diagrama de dispersão dos dados

Modelo Estatístico O modelo estatístico será o modelo de regressão simples com erros i.i.d. normais: \[ y_i=\beta_0+\beta_1 x_i + \epsilon_i, i=1,\ldots,n, \] onde \(\beta_0\): intercepto da linha de regressão com o eixo \(y\);
\(\beta_1\): coeficiente de inclinação da reta: é o nº de unidades em \(y\) que mudam para cada unidade da variável independente \(x\).
\(\epsilon_i\): erros aleatórios com distribuição normal: \(\epsilon_i \sim N(0,\sigma^2)\).

Estimadores de mínimos quadrados da regressão Encontrar \(\widehat{\beta_0}\) e \(\widehat{\beta_1}\) que minimizam a soma de quadrados dos erros:
\(S(\beta_0,\beta_1)=\sum \limits_{i=1}^{n} \epsilon_i^2 = \sum \limits_{i=1}^{n} (y_i-\beta_0-\beta_1 x_i)^2\).
Então temos as equações normais:

\[ \frac{\partial S} {\partial \beta_0} =0 \text{ e } \frac{\partial S} {\partial \beta_1} =0. \]

A solução é dada por:

\[ \begin{array}{lll} \widehat{\beta_1}&=&\frac{\sum \limits_{i=1}^{n} (x_i-\bar{x})(y_i-\bar{y})}{\sum \limits_{i=1}^{n} (x_i-\bar{x})^2}\\ \widehat{\beta_0}&=&\bar{y}-\widehat{\beta_1}\bar{x}. \end{array} \]

Com respeito a variância \(\sigma_2\): \[ \widehat{\sigma_2}=\frac{SQR}{n-2}, \] ou seja, a estimativa da variância é igual à soma dos quadrados dos resíduos sobre o número graus de liberdade do modelo. Os intervalos de confiança e testes de hipóteses para \(\beta_0\) e \(\beta_1\) são baseados na distribuição t-student.

Voltando ao R:
- Método dos mínimos quadrados: cálculo das estimativas passo a passo

soma_xx=sum((x-mean(x))^2)
soma_xy=sum((x-mean(x))*(y-mean(y)))
beta1_hat=soma_xy/soma_xx
beta0_hat=mean(y)-beta1_hat*mean(x)
print(paste0("beta0_hat=",beta0_hat," & beta1_hat=",beta1_hat))

## [1] "beta0_hat=-1.07090464547677 & beta1_hat=2.74083129584352"

• Método dos mínimos quadrados: utilizando a função lm:linear model
  • Resíduos aproximadamente normais com variância constante.

a=lm(y ~ x)
summary(a)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6333 -2.0128 -0.3741  2.0428  3.8851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.0709     2.7509  -0.389 0.707219    
## x             2.7408     0.4411   6.214 0.000255 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.821 on 8 degrees of freedom
## Multiple R-squared:  0.8284, Adjusted R-squared:  0.8069 
## F-statistic: 38.62 on 1 and 8 DF,  p-value: 0.0002555

anova(a)

## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x          1 307.247 307.247  38.615 0.0002555 ***
## Residuals  8  63.653   7.957                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(a)

shapiro.test(residuals(a))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(a)
## W = 0.93227, p-value = 0.4706