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