#lee los datos
pcancer <- read.table(file="cancerp.dat", header=T) 

#Aqui comienza la funcion
MixExp <- function(parametros, tipo, Ti, matriz.d){
n.var <- dim(matriz.d)[2]
J <- max(tipo)
nt <- length(tipo)
C.ij <- matrix(1:(nt * J), ncol = J) * 0
        for(i in 1:J) {
                C.ij[, i][tipo == i] <- 1
        }
uno <- rep(1, J)
ci <- C.ij %*% uno
Tij <- matrix(1:(nt * J), ncol = J) * 0 + 1
T.ij <- Tij * Ti
deltas <- parametros[1:((J-1)*n.var)]
betas <- parametros[((J-1)*n.var+1):((J-1)*n.var+n.var*J)]
dj.u <- matrix(1:(J * nt), ncol = J) * 0
Uij <- matriz.d
Vij <- matriz.d
if(n.var==1){
for(i in 1:(J - 1)) {
                                dj.u[, i] <- c(Uij * deltas[i])
                        }
}
else if(n.var > 1){
for(i in 1:(J - 1)) {
# cat("deltaj:",  deltas[((i-1)*n.var+1):(i*n.var)], "\n")
 dj.u[, i] <- c(Uij %*% deltas[((i-1)*n.var+1):(i*n.var)])
                        }
}
bj.v <- matrix(1:(J * nt), ncol = J)
if(n.var==1){
for(i in 1:J) {
                                bj.v[, i] <- c(Vij * betas[i])
                        }
}
else if(n.var > 1){
for(i in 1:J) {
# cat("betasj:", betas[((i-1)*n.var+1):(i*n.var)], "\n")
         bj.v[, i] <- c(Vij %*% betas[((i-1)*n.var+1):(i*n.var)])
                        }
}
p.ij <- matrix(1:(J * nt), ncol = J) * 0
exp.dj.u <- exp(dj.u)
sum.exps.i <- exp.dj.u %*% uno
for(i in 1:J) {
                        p.ij[, i] <- exp.dj.u[, i]/(sum.exps.i)
                }
l.ij <- exp(bj.v)
A <- sum(C.ij * (log(p.ij) + log(l.ij) - (l.ij * T.ij)))
F0.ij <- 1 - exp( - l.ij * T.ij)
p.F0.ij <- p.ij * F0.ij
sumJ.p.F0.ij <- p.F0.ij %*% uno
B <- sum((1 - ci) * log(1 - sumJ.p.F0.ij))
log.lik <- A + B
-1*log.lik
}
#aqui termina

mle.null <- nlm(MixExp,  c(-0.4421456,rep(0,14), -4.0580595, 
rep(0,14), -3.8259849, rep(0,14)),  
tipo=pcancer$tipo, Ti=pcancer$tiempo, 
matriz.d= model.matrix(~RX*(edad+peso.st+ACTIVO+HX+hg+medida)+sg, 
data=pcancer), hessian=T, gradtol = 1e-3)