#
# Análisis bayesiano de datos normales con una a priori semi-conjugada y el muestreo de Gibbs.  Los datos son los 
# famosos datos de Illingworth y miden las diferencias en velocidades de dos rayos de luz yendo la misma distancia
# en dos direcciones distintas.
#
rm(list=ls())
if (!require("bayess")) install.packages("bayess")
library(bayess)
data(normaldata)
y <- normaldata$x2
hist(y,probability=TRUE,xlab='y',ylab='f',main='Histograma de los datos de Illingworth con la densidad ajustada')
#
# Estadísticos suficientes.
#
n <- length(y)
ybar <- mean(y)
s2 <- var(y)
#
# A prioris.
#
m <- 0
c <- 0.001
a <- 0.5
b <- 0.5
#
# Cálculo de la verdadera distribución marginal a posteriori de mu mediante integración numérica.
#
fmu <- function(mu,m,c,a,b,n,ybar,s2){
smu <- sqrt((n-1)*s2/(n*(a+n-1)))
  f <- dnorm(mu,m,1/sqrt(c))*dt((mu-ybar)/smu,a+n-1)*smu
  return(f)
}
cmu <- integrate(fmu,lower=-Inf,upper=Inf,m=m,c=c,a=a,b=b,n=n,ybar=ybar,s2=s2)$value
#
# Cálculo de la verdadera distribución marginal a posteriori de phi mediante integración numérica.
#
fphi <- function(phi,m,c,a,b,n,ybar,s2){
  f <- dgamma(phi,0.5*(a+n),0.5*(b+(n-1)*s2))*dnorm(m,ybar,sqrt((n*phi+c)/(n*phi*c)))/sqrt(c*n*phi)
  return(f)
}
cphi <- integrate(fphi,lower=0,upper=Inf,m=m,c=c,a=a,b=b,n,ybar=ybar,s2=s2)$value
#
# Cálculo de la verdadera densidad conjunta, f(mu, phi| datos) = f(mu|datos) f(phi|mu, datos). 
#
fconjunta <- function(mu,phi,m,c,a,b,n,ybar,s2,cmu){
  f <- (fmu(mu,m,c,a,b,n,ybar,s2)/cmu)*dgamma(phi,0.5*(a+n),0.5*(b+(n-1)*s2+(n*(ybar-mu)^2)))
}
#
# Iteraciones para correr el muestreador.
#
burnin <- 1000
iters <- 10000
totits <- burnin+iters
#
# Valores a posteriori.
#
mupost <- rep(NA,iters)
phipost <- rep(NA,iters)
rejillay <- min(y)+c(0:1000)*(max(y)-min(y))/1000
fy <- rep(0,1001)
#
# Valores iniciales para el muestreador.
#
mu <- ybar
#
# Muestreador
#
for (i in 1:totits){
  phi <- rgamma(1,0.5*(a+n),0.5*(b+(n-1)*s2+n*(mu-ybar)*(mu-ybar)))
  mu <- rnorm(1,(c*m+n*phi*ybar)/(c+n*phi),1/sqrt(c+n*phi))
  if (i>burnin){
    mupost[i-burnin] <- mu
    phipost[i-burnin] <- phi
    fy <- fy+dnorm(rejillay,mupost[i-burnin],1/sqrt(phipost[i-burnin]))
  }
}
#
# Normalizar la densidad predictiva de y.
#
fy <- fy/iters
#
# Para calcular valores de sigma, transformamos.
#
sigmapost <- 1/sqrt(phipost)
#
# Dibujar la distribución predictiva estimada.
#
lines(rejillay,fy,type='l',lwd=2,col='orange')
#
# Distribuciones marginales a posteriori. El histograma muestra los valores generados por el MCMC.
# Las lineas sólidas son suavizaciones del histograma y las lineas de color magenta son las densidades
# estimadas a través de integración numérica.
#
hist(mupost,probability=TRUE,nclass=30,xlab=expression(mu),ylab='f',main='')
lines(density(mupost),lwd=2,col='red')
rejillamu = min(mupost)+c(0:1000)*(max(mupost)-min(mupost))/1000
lines(rejillamu,fmu(rejillamu,m,c,a,b,n,ybar,s2)/cmu,type='l',lwd=2,lty=2,col='magenta')
#
hist(phipost,probability=TRUE,nclass=30,xlab=expression(phi),ylab='f',main='')
lines(density(phipost),lwd=2,col='blue')
rejillaphi = min(phipost)+c(0:1000)*(max(phipost)-min(phipost))/1000
lines(rejillaphi,fphi(rejillaphi,m,c,a,b,n,ybar,s2)/cphi,type='l',lwd=2,lty=2,col='magenta')
#
# También miramos la distribución a posteriori de sigma.
#
hist(sigmapost,probability=TRUE,nclass=30,xlab=expression(sigma),ylab='f',main='')
lines(density(sigmapost),lwd=2,col='green')
#
# Distribución conjunta a posteriori de mu y phi.
#
if (!require("MASS")) install.packages("MASS")
library(MASS)
f <- kde2d(x=mupost,y=phipost,n=50)
plot(mupost,phipost,xlab=expression(mu),ylab=expression(phi))
#
# Mirar los contornos estimados ... 
#
contour(f,nlevels=10,lwd=2,col='brown',add=TRUE)
#
# ... y añadir los contornos calculados a través de integración numérica.
#
z <- matrix(rep(NA,1001*1001),nrow=1001)
for (i in 1:1001){
  for (j in 1:1001){
    z[i,j] <- fconjunta(mu=rejillamu[i],phi=rejillaphi[j],m,c,a,b,n,ybar,s2,cmu)
  }
}
contour(rejillamu,rejillaphi,z,nlevels=10,lwd=2,col='magenta',lty=2,add=TRUE)
#
# Pruebas de convergencia.
#
plot(mupost,type='l',lwd=2,col='red',ylab=expression(mu),xlab='iters')
plot(cumsum(mupost)/c(1:iters),type='l',lwd=2,col='red',ylab=expression(bar(mu)),xlab='iters')
acf(mupost,lwd=2,col='red',main='')
plot(sigmapost,type='l',lwd=2,col='red',ylab=expression(sigma),xlab='iters')
plot(cumsum(sigmapost)/c(1:iters),type='l',lwd=2,col='red',ylab=expression(bar(sigma)),xlab='iters')
acf(sigmapost,lwd=2,col='red',main='')