#
# Inferencia clásica y bayesiana para unos datos normales.
#
rm(list=ls())
#
if (!require("bayess")) install.packages("bayess")
library(bayess)
help(normaldata)
data(normaldata)
y <- normaldata$x2
n <- length(y)
hist(y,probability=TRUE,nclass=10,xlab='y',ylab='f',main="Histograma de y e densidades predictivas")
#
# Calculamos estimadores de mu, sigma^2
#
ybar=mean(y)
ybar # EMV de mu
s2=var(y)
s2 # Estimador insesgado de sigma^2
s <- sqrt(s2)
#
# Intervalo de confianza para mu
#
alpha=0.05
cimu <- qt(1-alpha/2,n-1)*s/sqrt(n)
cimu <- c(ybar-cimu,ybar+cimu)
cimu
#
# Intervalo de confianza para sigma^2
#
cisigma2 <- (n-1)*s2*c(1/qchisq(1-alpha/2,n-1),1/qchisq(alpha/2,n-1))
cisigma2
#
# Incluir la estimación clásica de la densidad de y en el histograma.
#
ygrid <- min(y)+c(0:1000)*(max(y)-min(y))/1000
fclass <- dnorm(ygrid,ybar,s)
lines(ygrid,fclass,type='l',lwd=2,col='blue')
#
# Fijar parámetros de la distribución a priori.
#
# No informativa: resultados de intervalos de credibilidad coinciden con los intervalos clásicos
# m <- 0
# c <- 0
# a <- 0
# b <- 0
#
# A priori propia
# 
m <- 0
c <- 2
a <- 1
b <- 1
#
# Calcular valores a posteriori.
#
cstar <- c+n
mstar <- (c*m+n*ybar)/cstar
astar <- ifelse(c>0,a+n,a+n-1)
bstar <- b+(n-1)*s2+c*n*((m-ybar)^2)/(c+n)
#
# Media a posteriori de mu.
#
mstar
#
# Media a posteriori de sigma2.
#
bstar/(astar-2)
#
# Intervalo de credibilidad para mu.
#
cimubayes <- mstar+sqrt(bstar/(astar*cstar))*c(qt(alpha/2,astar),qt(1-alpha/2,astar))
cimubayes
#
# Intervalo de credibilidad para sigma2.
#
cisigma2bayes <-1/c(qgamma(1-alpha/2,astar/2,bstar/2),qgamma(alpha/2,astar/2,bstar/2))
cisigma2bayes
#
# la densidad predictiva de y.
#
cpred <- cstar/(cstar+1)
fy <- sqrt(astar*cpred/bstar)*dt((ygrid-mstar)/sqrt(bstar/(astar*cpred)),astar)
lines(ygrid,fy,type="l",lwd=2,col="green")
#
# Observamos la diferencia entre la distribución "plug-in" y la densidad predictiva bayesiana cuando utilizamos 
# una a priori no informativa.
#