Segunda Lista de Exercícios:

Exercícios

Escolher 1 caso de prioris conjugadas - olhando no quadro Síntese da unidade de prioris conjugadas: “https://rpubs.com/liamorita/bayes_un2_conjugadas” ou nos exercícios 1 & 2 ao final da unidade, ou outro caso nos materiais de referência da disciplina.

  • Fazer a implementação numérica no R:

  • Gerar os dados;

  • Atribuir o número do RGA como semente para geração dos dados;

  • Fazer o gráfico da função de verossimilhança;

  • Encontrar a estimativa de máxima verossimilhança e exibir graficamente;

  • Exibir a tabela com os valores do parâmetro, juntamente com os valores da verossimilhança;

  • Atribuir uma distribuição a priori para o parâmetro de interesse (escolha livre para os parâmetros da priori);

  • Fazer o gráfico das três funções conjuntamente: priori, verossimilhança & posteriori;

Forma de entrega: Relatório em R markdown - pdf ou html.

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

Fazer a implementação numérica no R:

  • Gerar os dados;

  • Atribuir o número do RGA como semente para geração dos dados.

set.seed(2016113190) #cria uma semente única com com RGA
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]  5.6958  3.7998  2.8419  0.6255  1.9669 -2.0207  4.8482  1.1834  2.7007
[10]  4.8945

Estimativa Verossimilhança

sigma2_hat=sum((y-mean(y))^2)/n
print(paste0("Modo alternativo: sigma2_hat= ",sigma2_hat))
[1] "Modo alternativo: sigma2_hat= 4.903964758"

Função 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] 4.903945

$objective
[1] -22.13961

Exibir a tabela com os valores do parâmetro, juntamente com os valores da verossimilhança

sigma2=seq(0.01,20,0.01) #só assume valores positivos
temp <- data.frame(sigma2=sigma2, logL = logL(sigma2))
datatable(temp,options = list(
  columnDefs = list(list(className = 'dt-center', targets = 4)),
  pageLength = 20,
  lengthMenu = c(5, 10, 15, 20)
))

Gráfico da função de verossimilhança

ggplot(data = temp, aes(x = sigma2, y = logL)) + 
  geom_line() +
  theme_stata()+ # pacote ggthemes
  labs(x = "Sigma2",
       y = "logL",
       title = "Função de Verossimilhança",
       subtitle = "Estatística Bayesiana",
       caption = "Autor: Mateus Elias")+
  stat_peaks(col = "red")

Fazendo os cálculos no R: τ|y∼Gamma(27.6557,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= 27.65578859 e par_2= 5.5"

Atribuir uma distribuição a priori para o parâmetro de interesse

(escolha livre para os parâmetros da priori):

  • Fazer o gráfico das três funções conjuntamente: priori, verossimilhança & posteriori.
tau=seq(0,1.3,0.01)     
lambda=1
r=1/2
mu=3
constante=exp(25)

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