#---------------------------------------------------------------------
# Lista de Exercícios de Análise de Sobrevivência / MAE0514
#
#   Lista 1 / Entrega: 24/08/2005
#
# Fernando Henrique Ferraz Pereira da Rosa <academic at feferraz.net>
# Vagner Aparecido Pedro Junior <vapjunior at superig.com.br>
#----------------------------------------------------------------------
# 


rexp2 <- function(n,rate=1) {
    return(-log(runif(n,min=0,max=1))/rate)
}

simula.l1 <- function(n,media,p) {
    lambda <- 1/media
    falhas <- rexp2(n,rate=lambda)
    if (p == 0) { censuras <- falhas + 1 }
    else {
        censuras <- rexp2(n,lambda * p/(1-p))
    }
    

    Y <- pmin(falhas,censuras)
    r <- sum(falhas <= censuras)
    
    tempo.medio <- sum(Y)/r
    sd.tempo.medio <- sqrt(tempo.medio^2/r)
    mediana.est <- tempo.medio * log(2)
    sd.mediana.est <- sqrt(log(2)^2*tempo.medio^2/r)

    res <- structure(list(n=n,p=p,media=media,media.est=tempo.medio,
                          sd.media.est=sd.tempo.medio,mediana.est=mediana.est,
                          sd.mediana.est=sd.mediana.est),class='sim.l1.mae514')
    return(res)

}

print.sim.l1.mae514 <- function(x) {
    cat("\n      Simulação de um modelo exponencial com censuras e taxas\n",
        "    de falha independentes.\n\n")
    cat("Parametros de entrada: n = ",x$n," p = ",x$p," media = ",x$media,"\n")
    cat("\nmedia estimada: \t",x$media.est," (",x$sd.media.est,")\n",sep='')
    cat("mediana estimada: \t",x$mediana.est," (",x$sd.mediana.est,")\n",sep='')
}

as.data.frame.sim.l1.mae514 <- function(x) {
    x <- unclass(x)
    data.frame(x)
}

plot.sim.l1.mae514 <- function(x,xlim=c(0,x$media*10),ylim=c(0,1),ylab='Estimativa de máxima verossimilhança da função de sobrevivência',xlab='Tempo de sobrevivência',plot.true=TRUE) {
    func <- function(y) { exp(-1/x$media.est * y) }
    plot(func,xlab=xlab,ylab=ylab,main=paste(" n = ",x$n," p = ",x$p," media = ",x$media),ylim=ylim,xlim=xlim)
    if (plot.true == TRUE) {
        func.true <-  function(y) { exp(-1/x$media * y) }
        plot(func.true,lty=2,col=4,xlim=xlim,add=TRUE)
        legend(x='topright',lty=c(1,2),col=c(1,4),
                legend=c('estimativa','real'),inset=0.05)

    }
}