#
# El muestreo Gibbs para datos censurados.  Los datos son tiempos de vida modelados con una distribución log-normal.
# Las maquinas con cens = 1 han fallado  y las maquinas con cens = 0 siguen funcionando.
#
rm(list=ls())
#
y <- c(3.4,2.9,1.2,1.4,3.2,1.8,4.6,1.7,2.0,1.4,2.8,0.6)
cens <- c(0,0,1,0,0,0,0,1,1,1,0,1)
#
# Dividimos los datos en censurados y no censurados.
#
ysincens <- y[cens==0]
m <- length(ysincens)
n <- length(y)
ycens <- y[cens==1]
#
# Fijar las iteraciones para el muestreador de Gibbs.
#
burnin <- 10000
iters <- 10000
totits <- burnin+iters
#
# Definir los valores a posteriori.
#
postmu <- rep(NA,iters)
postphi <- rep(NA,iters)
#
# Valores iniciales: usar estimaciones de la media y varianza basadas en los datos sin censurar. Inicializamos los
# datos censurados como si no lo fuesen.
#
mu <- mean(log(ysincens))
phi <- 1/var(log(ysincens))
ycompleto <- c(ysincens,ycens)
#
# En el muestreados, utilizamos el paquete truncdist para muestrear de una log-normal truncada.
#
if (!require("truncdist")) install.packages("truncdist")
library(truncdist)
#
# Muestreador.
#
for (i in 1:totits){
  #
  # Generar los tiempos de vida completos para los datos censurados
  #
  for (j in (m+1):n){
    ycompleto[j] <- rtrunc(1,"lnorm",a = ycens[j-m],mean=mu,sd=1/sqrt(phi))
  }
  #
  # Estadísticos suficientes dada la muestra completa.
  #
  mediacompleto <- mean(log(ycompleto))
  #
  # Generar de la distribución a posteriori condicional de mu|phi, datos completos.
  #
  mu <- rnorm(1,mediacompleto,sqrt(1/(n*phi)))
  #
  # Generar de la distribución a posteriori condicional de phi|mu, datos completos.
  #
  phi <- rgamma(1,0.5*n,0.5*sum((log(ycompleto)-mu)^2))
  if (i > burnin){
    postmu[i-burnin] <- mu
    postphi[i-burnin] <- phi
  }
}
#
# Comprobar la convergencia.
#
ts.plot(postmu,xlab='iters',ylab='mu')
ts.plot(cumsum(postmu)/c(1:iters),xlab='iters',ylab='E[mu|datos]')
acf(postmu,main='ACF de mu')
ts.plot(postphi,xlab='iters',ylab='phi')
ts.plot(cumsum(postphi)/c(1:iters),xlab='iters',ylab='E[phi|datos]')
acf(postphi,main='ACF de phi')
#
# (Los acf sugieren que debemos hacer algo de thinning pero no es demasiado serio).  
#
# Densidades a posteriori y resumenes
#
plot(density(postmu),xlab=expression(mu),ylab='f',main='Densidad a posteriori de mu')
mean(postmu)
quantile(postmu,c(0.025,0.5,0.975))
plot(density(postphi),xlab=expression(phi),ylab='f',main='Densidad a posteriori de phi')
mean(postphi)
quantile(postphi,c(0.025,0.5,0.975))