R/Internal.R

Defines functions .estim.plusTraj .estimNP_CensureFin .estimNP_aucuneCensure .write.results .calcul_NP.q.matrice .calcul_NP.q .p.hat .calcul.Niu .calcul.Niujv .calcul.Niujv.matrice .kernel_param_f .kernel_param_fi .kernel_param_fj .kernel_param_fij .calcul_q .limit.law .stationnary.law .comptage .donneYU .donneSeq .donneJT .decode .code

## __________________________________________________________
##
## .code
##
## __________________________________________________________
##

.code <- function(sequence,E)
{
  seq<-sequence
  s<-length(E)
  for (i in 1:s)
  {
    seq[seq==E[i]]<-i
  }
  seq<-as.numeric(seq)
  return(seq)
}


## __________________________________________________________
##
## .decode
##
## __________________________________________________________
##

.decode <-function(seq,E)
{
  sequence<-seq
  s<-length(seq)
  for (i in 1:s)
  {
    sequence[i]<-E[seq[i]]
  }
  return(sequence)    
}



## __________________________________________________________
##
## .donneJT
##
## __________________________________________________________
##

.donneJT <-function(seq)
{
  J<-NULL
  T<-NULL
  n<-length(seq)
  i1<-1
  i2<-2
  
  while (i1 <= n){
    while (i2<=n && (seq[i2]==seq[i1])){
      i2<-i2+1
    }
    J<-cbind(J,seq[i1])
    T<-cbind(T,i2)
    i1<-i2
    i2<-i1+1
  }
  JT <-rbind(J,T)
  colnames(JT)<-NULL
  return(JT)
}

## __________________________________________________________
##
## .donneSeq
##
## __________________________________________________________
##

.donneSeq <- function(J,T)
{
  seq<-NULL
  i<-1
  n<-length(T)
  t<-1
  while (i<= n)
  {
    k<-T[i]-t
    seq <-c(seq,rep(J[i],k))
    # for (j in 1:k)
    # {
    #   seq<-cbind(seq,J[i])
    # }
    t<-T[i]
    i<- i+1
  }
  
  return(seq)
}

## __________________________________________________________
##
## .donneYU
##
## __________________________________________________________
##

.donneYU<-function(JT, E){
  J<-.code(JT[1,],E)
  T<-as.numeric(JT[2,])
  L<-c(1,T)
  Y<-NULL
  U<-NULL
  n<-length(J)
  
  for (i in 1:n){
    for (j in 0:(L[i+1]-L[i]-1)){
      Y<-cbind(Y,J[i])
      U<-cbind(U,j)
    }
  }
  YU<-rbind(Y,U)
  colnames(YU)<-NULL
  return (YU)
}

## __________________________________________________________
##
## .comptage
##
## __________________________________________________________
##
.comptage <- function(J,L,S,Kmax){
  
  if(!is.list(J) || !is.list(L)){
    stop("J and L should be lists")
  }
  
  nbSeq = length(J)
  sNijk<-array(0,c(S,S,Kmax))
  Nstart = c()
  Nbij = array(0, c(S,S,Kmax))
  Nei = matrix(0, nrow = S, ncol = Kmax)
  Nbj = matrix(0, nrow = S, ncol = Kmax)
  Nb = vector( length = Kmax)
  Ne = vector( length = Kmax)  
  LsNijk<-list()
  LNstart = list()
  LNbij = list()
  LNei = list()
  LNbj = list()
  LNb = list()
  LNe = list()
  LNijk = list()
  for (k in 1:nbSeq){
    Ji = J[[k]]
    Li = L[[k]]
    M = length(Ji)
    Nstart = c(Nstart, Ji[1])
    Nbij[Ji[1],Ji[2],Li[1]] = Nbij[Ji[1],Ji[2],Li[1]] + 1
    LNbij[[k]] = Nbij
    Nei[Ji[M],Li[M]] = Nei[Ji[M],Li[M]] + 1
    LNei[[k]] = Nei
    # Nbj[Ji[2],Li[1]] = Nbj[Ji[2],Li[1]] + 1
    # Nb[Li[1]] = Nb[Li[1]] + 1
    # Ne[Li[M]] = Ne[Li[M]] + 1
    Nijk<-array(0,c(S,S,Kmax))
    for (i in 2:M){
      Nijk[Ji[i-1],Ji[i],Li[i-1]]<-Nijk[Ji[i-1],Ji[i],Li[i-1]]+1
    }
    sNijk <- sNijk + Nijk
    LNijk[[k]] = Nijk
  }
  
  LNij = list()
  Nij<-rowSums(sNijk,dims=2)
  for( k in 1:nbSeq ){
    LNij[[k]] = rowSums(LNijk[[k]], dims = 2)
  }
  
  diago = seq(1, S*S, by = S+1)
  Nij.temp = as.vector(Nij)[-diago]
  if( length(which(Nij.temp==0))!=0) {
    warning("Warning : missing transitions")
  }
  Nbj = colSums(Nbij)
  Ne = colSums(Nei)
  Nb = colSums(Nbj)
  Nbi = apply(Nbij, 3, rowSums)
  Nebi = Nei + Nbi
  Neb = Ne + Nb
  Nstarti = as.vector(count(seq = Nstart, wordsize = 1, alphabet = 1:S))
  LNstarti = list()
  for (k in 1:nbSeq){
    LNstarti[[k]] = as.vector(count(seq = Nstart[k], wordsize = 1, alphabet = 1:S))
  }
  Ni<-rowSums(Nij)
  if( length(which(Ni==0))!=0) {
    warning("Warning : missing letters")
  }
  Nj<-colSums(Nij)
  if( length(which(Nj==0))!=0) {
    warning("Warning : missing letters")
  }
  N<-sum(Nij)
  LNk = list()
  Nk<-colSums(sNijk,dims=2)
  for (k in 1:nbSeq){
    LNk[[k]] = colSums(LNijk[[k]], dims = 2)
  }
  LNik = list()
  Nik<-apply(sNijk,3,rowSums) #  i in raws and k in columns
  for (k in 1:nbSeq){
    LNik[[k]]<-apply(LNijk[[k]],3,rowSums)
  }
  LNjk = list()
  Njk<-colSums(sNijk) #  j in raws and k in columns
  for (k in 1:nbSeq){
    LNjk[[k]]<-colSums(LNijk[[k]])
  }
  return(list(Nijk=sNijk,Nij=Nij, Ni=Ni, Nj=Nj, N=N, Nk=Nk, Nik=Nik, Njk=Njk, Nstarti = Nstarti, Nbij = Nbij, Nei = Nei,
              Nbj = Nbj, Neb = Neb, Nebi = Nebi, Ne = Ne, Nb = Nb, Nbi = Nbi, LNijk = LNijk, LNij = LNij, LNij = LNij, LNik = LNik, LNjk = LNjk, 
              LNk = LNk, Nstart = LNstarti))
  
}

## __________________________________________________________
##
## .strationnary.law
##
## __________________________________________________________
##
# tpm: transition matrix p_ij = P(x_{t+1} = j | x_t = i)
.stationnary.law <- function(Pest){
  ## nb states
  m <- dim(Pest)[1]
  ## computation 
  A <- t(Pest) - diag(1, m, m)
  A[m, ] <- 1
  b <- c(rep(0, (m-1)), 1)
  stat.law <- solve(A, b)
  ## results
  return(stat.law)
}


## __________________________________________________________
##
## .limit.law
##
## __________________________________________________________
##
.limit.law <- function(q, pij){
  
  if ( dim(q)[1] != dim(pij)[1] && dim(pij)[2] != dim(q)[2] ){
    stop("The state number is not valid")
  }
  
  Kmax = dim(q)[3]
  S = dim(pij)[1]
  
  hj = apply(X = q, MARGIN = c(1,3), FUN = sum)
  ksup1 = rep(1:Kmax, S)
  mj_int = matrix(ksup1*as.vector(t(hj)), nrow = S, ncol = Kmax, byrow = TRUE)
  mj = rowSums(mj_int)
  statLaw = .stationnary.law(Pest = pij)
  
  Pi.j = statLaw*mj/sum(statLaw*mj)
  
  return(Pi.j)  
}


################################################################
## Computation of semi-Markov kernel (non-parametric case) ##
################################################################
## __________________________________________________________
##
## .calcul_q
##
## __________________________________________________________
##
.calcul_q = function(p,f,Kmax){
  
  dim.q=dim(f)
  q<-array(0,dim.q)
  for (k in 1:Kmax) {
    q[,,k] = p*f[,,k]
  }
  
  return(q)
  
}


############################################################
## Computation odf semi-Markovian kernel (parametric case)##
############################################################
### CASE fij ###
## __________________________________________________________
##
## .kernel_param_fij
##
## __________________________________________________________
##
# Computation of the kernel with the param just estimated
.kernel_param_fij <- function(distr, param, Kmax, pij, S){
  # get all the distributions in a vector
  lois = as.vector(distr)
  S2 = S*S
  fik = matrix(0, nrow = S2, ncol = Kmax)
  param = as.vector(param)
  diag = seq(from = 1, to = S*S, by = S+1)
  
  # print(S)
  # print(lois)
  # print(param)
  
  if ( "pois" %in% lois ){
    i = which(distr == "pois")
    for (j in i){
      fik[j,] = dpois(0:(Kmax-1), lambda = param[j])  
    }
  }
  if ( "geom" %in% lois ){
    i = which(distr == "geom")
    for (j in i){
      fik[j,] = dgeom(0:(Kmax-1), prob = param[j])
    }
  }
  if ( "nbinom" %in% lois ){
    i = which(distr == "nbinom")
    for (j in i){
      # print(j)
      # print(param[j])
      # print(param[j+S*S])
      fik[j,] = dnbinom(0:(Kmax-1), size = param[j], mu = param[j+S*S])
      if (fik[j,1] == 1){ fik[j,-1] = 0}else if (j %in% diag){ fik[j,] = 0}
      # print(fik)
      
    }
  }
  if ( "dweibull" %in% lois ){
    i = which(distr == "dweibull")
    for (j in i){
      fik[j,] = ddweibull(1:Kmax, q = param[j], beta = param[j+S*S])
      if (fik[j,1] == 1){ fik[j,-1] = 0}else if (j %in% diag){ fik[j,] = 0}
    }
  }
  fijk = array(fik, c(S,S,Kmax))
  q = array(pij, c(S,S,Kmax))*fijk
  
  return(list(q = q, f = fijk))
}

## __________________________________________________________
##
## .kernel_param_fj
##
## __________________________________________________________
##
### CASE fj ###
# Computation of the kernel with the param just estimated
.kernel_param_fj <- function(distr, param, Kmax, pij, S){
  
  lois = as.vector(distr)
  fik = matrix(data = 0, nrow = S, ncol = Kmax)
  
  if ( "pois" %in% lois ){
    i = which(distr == "pois")
    for (j in i){
      fik[j,] = dpois(0:(Kmax-1), lambda = param[j,1])
    }
    # i = which(distr == "pois")
    # fik[i,] = dpois(0:(Kmax-1), lambda = param[i])
  }
  if ( "geom" %in% lois ){
    i = which(distr == "geom")
    for (j in i){
      fik[j,] = dgeom(0:(Kmax-1), prob = param[j,1])
    }
  }
  if ( "nbinom" %in% lois ){
    i = which(distr == "nbinom")
    for (j in i){
      fik[j,] = dnbinom(0:(Kmax-1), size = param[j,1], mu = param[j,2])
    }
  }
  if ( "dweibull" %in% lois ){
    i = which(distr == "dweibull")
    for (j in i){
      fik[j,] = ddweibull(1:Kmax, q = param[j,1], beta = param[j,2])  
    }
  }
  
  f = rep(as.vector(t(fik)), each = S)
  fmat = matrix(f, nrow = Kmax, ncol =S*S, byrow = T)
  fk = array(as.vector(t(fmat)), c(S,S,Kmax))
  # fi.k = matrix(rep(as.vector(t(fik)), S), nrow = S*S, ncol = Kmax, byrow = TRUE)
  
  q = array(pij, c(S,S,Kmax))*fk
  q[which(is.na(q))]<-0
  
  return(list(q = q, f = fk))
}

## __________________________________________________________
##
## .kernel_param_fi
##
## __________________________________________________________
##
### CASE fi ###
# Computation of the kernel with the param just estimated
.kernel_param_fi <- function(distr, param, Kmax, pij, S){
  
  lois = as.vector(distr)
  fik = matrix(data = 0, nrow = S, ncol = Kmax)
  
  if ( "pois" %in% lois ){
    i = which(distr == "pois")
    for (j in i){
      fik[j,] = dpois(0:(Kmax-1), lambda = param[j,1])
    }
    # i = which(distr == "pois")
    # fik[i,] = dpois(0:(Kmax-1), lambda = param[i])
  }
  if ( "geom" %in% lois ){
    i = which(distr == "geom")
    for (j in i){
      fik[j,] = dgeom(0:(Kmax-1), prob = param[j,1])
    }
  }
  if ( "nbinom" %in% lois ){
    i = which(distr == "nbinom")
    for (j in i){
      fik[j,] = dnbinom(0:(Kmax-1), size = param[j,1], mu = param[j,2])
    }
  }
  if ( "dweibull" %in% lois ){
    i = which(distr == "dweibull")
    for (j in i){
      fik[j,] = ddweibull(1:Kmax, q = param[j,1], beta = param[j,2])  
    }
  }
  
  f = rep(as.vector(t(fik)), each = S)
  fmat = matrix(f, nrow = Kmax, ncol =S*S, byrow = TRUE)
  # fmat2 = matrix(as.vector(t(fmat)))
  fk = array(as.vector(t(fmat)), c(S,S,Kmax))
  fi.k = apply(X = fk,MARGIN =  c(1,3),FUN =  t)
  # fi.k = matrix(rep(as.vector(t(fik)), S), nrow = S*S, ncol = Kmax, byrow = TRUE)
  
  q = array(pij, c(S,S,Kmax))*fi.k
  q[which(is.na(q))]<-0
  
  return(list(q = q, f = fk))
}


## __________________________________________________________
##
## .kernel_param_f
##
## __________________________________________________________
##
### CASE f ###
.kernel_param_f <- function(distr, param, Kmax, pij, S){
  if ( distr == "pois" ) {
    f = dpois(0:(Kmax-1), lambda = param[1])
  }
  if ( distr == "geom" ) {
    f = dgeom(0:(Kmax-1), prob = param[1])
  }
  if ( distr == "dweibull" ) {
    f = ddweibull(1:Kmax, q = param[1], beta = param[2])
  }
  if ( distr == "nbinom" ){
    f = dnbinom(0:(Kmax-1), size = param[1], mu = param[2])
  }
  
  f = rep(f, each = S*S)
  fmat = matrix(f, nrow = Kmax, ncol =S*S, byrow = TRUE)
  fk = array(as.vector(t(fmat)), c(S,S,Kmax))
  # fk = array(fmat, c(S,S,Kmax)) 
  q = array(pij, c(S,S,Kmax))*fk
  q[which(is.na(q))]<-0
  
  return(list(q = q, f = fk))
}



##################################
## NON-PARAMETRIC (couple MC) ##
##################################
## __________________________________________________________
##
## .calcul.Niujv.matrice 
##
## __________________________________________________________
##
##
.calcul.Niujv.matrice <- function(Y,U,S,Kmax){
  
  if(!is.list(Y) || !is.list(U)){
    stop("Y and U should be lists")
  }
  
  nbSeq = length(Y)
  sNiujv<-matrix(0, nrow=S*Kmax, ncol=S*Kmax)
  for (k in 1:nbSeq){
    Yi = Y[[k]]
    Ui = U[[k]]
    M = length(Yi)
    Niujv=matrix(0, nrow=S*Kmax, ncol=S*Kmax)
    for (i in 2:M){
      if((Yi[i-1] != Yi[i] && Ui[i] == 0) || (Yi[i-1] == Yi[i] && Ui[i] == Ui[i-1] + 1)){
        Niujv[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+Ui[i-1],Yi[i]+(Yi[i]-1)*(Kmax-1)+Ui[i]]<-
          Niujv[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+Ui[i-1],Yi[i]+(Yi[i]-1)*(Kmax-1)+Ui[i]]+1
      }
    }
    
    #     print(Niujv)
    sNiujv <- sNiujv + Niujv
  }
  
  return(sNiujv)
  
}

## __________________________________________________________
##
## .calcul.Niujv 
##
## __________________________________________________________
##
.calcul.Niujv <- function(Y,U,S,Kmax){
  
  if(!is.list(Y) || !is.list(U)){
    stop("Y and U should be lists")
  }
  
  nbSeq = length(Y)
  for (k in 1:nbSeq){
    ## Shift from aphabet {1,2,3,4,...} to {0,1,2,3,...} 
    Yi = Y[[k]]-1
    Ui = U[[k]]
    M = length(Yi)
    Niujv=rep(0, (Kmax+1)*S*(S-1)+S*(Kmax))
    for (i in 2:M){
#       if((Yi[i-1] != Yi[i] && Ui[i] == 0) || (Yi[i-1] == Yi[i] && Ui[i] == Ui[i-1] + 1)){
#         Niujv[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+Ui[i-1],Yi[i]+(Yi[i]-1)*(Kmax-1)+Ui[i]]<-
#           Niujv[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+Ui[i-1],Yi[i]+(Yi[i]-1)*(Kmax-1)+Ui[i]]+1
#       }
      if( Yi[i-1] != Yi[i] ){
        if( Ui[i] != 0 ){
          stop("erreur")
        }
        if( Yi[i] > Yi[i-1] ){
          Niujv[Yi[i-1]*(Kmax+1)*(S-1)+(Kmax+1)*(Yi[i]-1)+Ui[i-1]+1] =
            Niujv[Yi[i-1]*(Kmax+1)*(S-1)+(Kmax+1)*(Yi[i]-1)+Ui[i-1]+1] + 1
        }else{
          Niujv[Yi[i-1]*(Kmax+1)*(S-1)+(Kmax+1)*Yi[i]+Ui[i-1]+1] = 
            Niujv[Yi[i-1]*(Kmax+1)*(S-1)+(Kmax+1)*Yi[i]+Ui[i-1]+1] + 1
        } 
      }else{
        if( Ui[i] != Ui[i-1] + 1 ){
          stop("erreur")
        }
        Niujv[(Kmax+1)*S*(S-1)+Yi[i-1]*Kmax+Ui[i-1]+1] = Niujv[(Kmax+1)*S*(S-1)+Yi[i-1]*Kmax+Ui[i-1]+1] + 1 
      }
    }
    
    #     print(Niujv)
#     sNiujv <- sNiujv + Niujv
  }
  
  return(Niujv)
  
}

## __________________________________________________________
##
## .calcul.Niu
##
## __________________________________________________________
##
.calcul.Niu = function(Kmax,S,Niujv){

  Niu<-matrix(0, nrow=S, ncol=Kmax+1)
  for (i in 0:(S-1)){
    for (u in 0:(Kmax-1)){
      Niu[i+1,u+1] = sum(Niujv[i*(Kmax+1)*(S-1)+(Kmax+1)*(0:(S-2))+u+1]) + 
        Niujv[(Kmax+1)*S*(S-1)+i*Kmax+u+1]
          
    }
    Niu[i+1,Kmax] = sum(Niujv[i*(Kmax+1)*(S-1)+(Kmax+1)*(0:(S-2))+u+1])
  }  
    
  return (Niu)
  
}

## __________________________________________________________
##
## .p.hat 
##
## __________________________________________________________
##
.p.hat = function(Kmax,S,Niujv,Niu){
  
  p = rep(0, (Kmax+1)*S*(S-1)+S*(Kmax))
  for (i in 0:(S-1)){
    for (u in 0:(Kmax-1)){
      for (j in 0:(S-2)){
        #         for (v in 0:(Kmax-1)){
        if( i != j ){
          if( i < j ){
            if( Niu[i+1,u+1] == 0 ){
              break
            }
            p[(Kmax+1)*i*(S-1)+(j-1)*(Kmax+1)+u+1] = 
              Niujv[(Kmax+1)*i*(S-1)+(j-1)*(Kmax+1)+u+1]/Niu[i+1,u+1]
            
          }else{
            if( Niu[i+1,u+1] == 0 ){
              break
            }
            p[(Kmax+1)*i*(S-1)+j*(Kmax+1)+u+1] = 
              Niujv[(Kmax+1)*i*(S-1)+j*(Kmax+1)+u+1]/Niu[i+1,u+1]
            
          }
        }else{
          if( Niu[i+1,u+1] == 0 ){
            break
          }
          #             if( v == u + 1 ){
          
          p[(Kmax+1)*S*(S-1)+i*Kmax+u+1] = Niujv[(Kmax+1)*S*(S-1)+i*Kmax+u+1]/Niu[i+1,u+1]
          
          #             }
        }
        
        #           }
      }
    }      
    
  }
  return (p)
}  


## __________________________________________________________
##
## .calcul_NP.q
##
## __________________________________________________________
##
.calcul_NP.q = function(S,Kmax,p){
  
  q = array(0, c(S,S,Kmax))
  Pi = 1
  for (i in 0:(S-1)){
    for (j in 0:(S-1)){
      if (i!=j){
      for (k in 1:Kmax){
        Pr = 1        
       
        if ( i > j ){
          Pi = p[(Kmax+1)*i*(S-1)+j*(Kmax+1)+k] #(k-1)+1
          #           cat("Pi :")
          #           print(Pi)
          #           cat("Pr :")
          #           print(Pr)
          #           cat("==================================")
          #           cat("\n")

        }
        if ( i < j ){
          Pi = p[(Kmax+1)*i*(S-1)+(j-1)*(Kmax+1)+k] #(k-1)+1
          #           cat("Pi :")
          #           print(Pi)
          #           cat("Pr :")
          #           print(Pr)
          #           cat("==================================")
          #           cat("\n")

        }
        Pr=Pi*Pr
        if( k >= 2 ){
          for (t in 0:(k-2)){
            
            Pr = Pr * p[(Kmax+1)*S*(S-1)+ i*Kmax + t+1]
            
          }
        }
        
        q[i+1,j+1,k]=Pr
        
      }
      }
    }
  }
  return (q)
}

## __________________________________________________________
##
## .calcul.NP.q.matrice
##
## __________________________________________________________
##
.calcul_NP.q.matrice = function(Y,U,p,S,Kmax){
  if(!is.list(Y) || !is.list(U)){
    stop("Y and U should be lists")
  }
  
  nbSeq = length(Y) 
  q=array(0, c(S,S,Kmax))
  for (l in 1:nbSeq){
    Yi = Y[[l]]
    Ui = U[[l]]
    M = length(Yi)
    
    for (i in 2:M){
      if((Yi[i-1] != Yi[i] && Ui[i] == 0) || (Yi[i-1] == Yi[i] && Ui[i] == Ui[i-1] + 1)){
        for (k in 1:Kmax){
          Pr=1
          if( k >= 2){
            for (j in 1:(k-1)){
              Pr = Pr*p[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+j-1,Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+j]
            }
          }
          q[Yi[i-1],Yi[i],k]=Pr*p[Yi[i-1]+(Yi[i-1]-1)*(Kmax-1)+k-1,Yi[i]+(Yi[i]-1)*(Kmax-1)]
        }
      }
    }
    
  }
  return (q)
}

## __________________________________________________________
##
## .write.results
##
## __________________________________________________________
##
.write.results = function(TypeSojournTime, cens.beg, cens.end, p, param){
  date.res = paste("results", Sys.Date(), sep="_")
  File.out = paste(date.res, "txt", sep=".")
  i = 1
  while(file.exists(File.out)){
    # print(File.out)
    ## get the file name before extension
    File.txt = unlist(strsplit(File.out, split = "[.]"))[1]
    if (i>1) {
      ## get the file name without version number
      File.1 = paste(unlist(strsplit(File.txt, split = "_"))[1], unlist(strsplit(File.txt, split = "_"))[2], sep="_")
      Filei = paste(File.1, i, sep="_")
      # Filei = paste(File.int, "txt", sep=".")
    }else{
      Filei = paste(File.txt, i, sep="_")
    }
    File.out = paste(Filei, "txt", sep=".")
    i = i + 1
  }
  sink(File.out)
  cat("The Results of estimation for one or several sequences \n")
  cat("\n")
  cat("The sojourn time chosen : ")
  cat("\n")
  print(TypeSojournTime)
  # cat("The matrix of the distributions")
  # cat("\n")
  # print(distr)
  cat("Censoring at the beginnig")
  cat("\n")
  if(cens.beg == 1){print("yes")}else{print("no")}
  cat("Censoring at the end")
  cat("\n")
  if(cens.end == 1){print("yes")}else{print("no")}
  cat("The transition probabilities : ")
  cat("\n")
  print(p)
  cat("The estimated parameters : ")
  cat("\n")
  print(param)
  sink()
} 


## __________________________________________________________
##
## .estimNP_aucuneCensure
##
## __________________________________________________________
##
####################
##### INPUT: 
## seq: list of sequences
## E: state space
## TypeSojournTime: Type of Sojourn time ("fij","fi","fj","f")
## distr = "NP"
## No censoring
## All trajectories with the same type of censoring
#####################
##### OUTPUT: 
## p: transition matrix
## f: array
## q: array
.estimNP_aucuneCensure<-function(file = NULL, seq, E, TypeSojournTime="fij", distr="NP"){
  
  
  if(is.null(file) && is.null(seq)){
    stop("One of the two parameters file or seq should be given")
  }
  
  if(!is.null(file)){
    fasta = read.fasta(file = file)
    seq = getSequence(fasta)
    attr = getAnnot(fasta)
  }
  
  
  if(!is.list(seq)){
    stop("The parameter seq should be a list")
  }
  
  ## States
  TYPES<-c("fij", "fi", "fj", "f")
  type<-pmatch(TypeSojournTime,TYPES)
  
  ## Length of the state space
  S<-length(E)
  
  ## All sequences in a vector
  nbSeq<-length(seq)
  #  seq1<-c()
  #  for (k in 1:nbSeq){
  #    seq1<-c(seq1, unlist(seq[[k]]))
  #  }
  
  Kmax = 0
  J<-list()
  T<-list()
  L<-list()
  M<-list()
  
  for (k in 1:nbSeq){
    ## Manipulations
    JT<-.donneJT(unlist(seq[[k]]))
    J[[k]]<-.code(JT[1,],E)
    T[[k]]<-as.numeric(JT[2,])
    L[[k]]<-T[[k]] - c(1,T[[k]][-length(T[[k]]-1)]) ## sojourn time
    Kmax <- max(Kmax, max(L[[k]])) 
  }
  
  
  ## J: list
  ## L: list
  ## S: state space size
  ## Kmax: maximal sojourn time
  ################
  ## get the counts
  res<-.comptage(J,L,S,Kmax)
  
  Nij = res$Nij
  Ni = res$Ni
  Nj = res$Nj
  N= res$N
  Nk = res$Nk
  Nik = res$Nik
  Njk = res$Njk
  Nijk = res$Nijk
  Nstart = res$Nstarti
  Nbij = res$Nbij
  Nei = res$Nei
  Nbj = res$Nbj
  Neb = res$Neb
  Nebi = res$Nebi
  Nbi = res$Nbi
  Nb = res$Nb
  Ne = res$Ne
  #########
  
  
  
  if(is.na(distr)){
    stop("invalid distribution")
    
    ## Non-parametric
  }else if(distr == "NP"){
    
    if( type == 4 ){
      # Nij<-rowSums(Nijk,dims=2)
      # Ni<-rowSums(Nij)
      # N<-sum(Nij)
      mat.Nk<-matrix(1,nrow=S,ncol=S)
      diag(mat.Nk)<-0
      Nk<-array(mat.Nk,c(S,S,Kmax))
      
      ## Estimators
      if(length(which(Ni==0))!=0){ # Presence of all states
        
        print("Warning : All the states are not observed")
        cat("\n")
        
        p<-Nij/(as.vector(Ni))
        p[which(is.na(p))]<-0
        
        f <- Nk/N
        f[which(is.na(f))]<-0
        
        q<-Nijk/array(Ni,c(S,S,Kmax))
        q[which(is.na(q))]<-0
        
      } else {
        diago = seq(1, S*S, by = S+1)
        Nij.temp = as.vector(Nij)[-diago]
        if( length(which(Nij.temp==0))!=0) # Presence of all transitions
        {
          print("Warning: The transitions are not all observed")
          cat("\n")
          
          p<-Nij/(as.vector(Ni))
          
          f <- Nk/N
          f[which(is.na(f))]<-0
          
          q<-Nijk/array(Ni,c(S,S,Kmax))                    
        } else {
          p<-Nij/(as.vector(Ni))
          f <- Nk/N
          q<-Nijk/array(Ni,c(S,S,Kmax))
        }
      }
      # return (list(p,f,q))
      
    }else {
      
      if( type == 2 ){
        # Nij<-rowSums(Nijk,dims=2)
        Nik<-vector("list", Kmax)
        N<-vector("list", Kmax)
        list.Nijk<-lapply(seq(dim(Nijk)[3]), function(x) Nijk[ , , x])
        for (k in 1:Kmax){
          N[[k]]<-apply(list.Nijk[[k]],1,sum)
          Nik[[k]]<-matrix(unlist(N[[k]]),4,4)
          diag(Nik[[k]])<-0
        }
        Nik<-array(unlist(Nik), c(4,4,Kmax))
        Ni<-rowSums(Nij)
        
        ## Estimators
        if(length(which(Ni==0))!=0){ # Presence of all states
          
          print("Warning: Not all the states are observed")
          cat("\n")
          
          p<-Nij/(as.vector(Ni))
          p[which(is.na(p))]<-0
          
          f<-Nik/Ni
          f[which(is.na(f))]<-0
          
          q<-Nijk/array(Ni,c(S,S,Kmax))
          q[which(is.na(q))]<-0
          
        } else {
          diago = seq(1, S*S, by = S+1)
          Nij.temp = as.vector(Nij)[-diago]
          if( length(which(Nij.temp==0))!=0) # Presence of all transitions
          {
            print("Warning: The transitions are not all observed")
            cat("\n")
            
            p<-Nij/(as.vector(Ni))
            
            f<-Nik/Ni
            f[which(is.na(f))]<-0
            
            q<-Nijk/array(Ni,c(S,S,Kmax))                    
          } else {
            p<-Nij/(as.vector(Ni))
            f<-Nik/Ni
            q<-Nijk/array(Ni,c(S,S,Kmax))
          }
        }
        # return (list(p,f,q))
      }
      
      if( type == 1 ){
        # Nij<-rowSums(Nijk,dims=2)
        # Ni<-rowSums(Nij)
        ## Estimators
        
        if(length(which(Ni==0))!=0){ # Presence of all states
          
          print("Warning: The transitions are not all observed")
          cat("\n")
          
          p<-Nij/(as.vector(Ni))
          p[which(is.na(p))]<-0
          
          f <- Nijk/array(Nij,c(S,S,Kmax))
          f[which(is.na(f))]<-0
          
          q<-Nijk/array(Ni,c(S,S,Kmax))
          q[which(is.na(q))]<-0
          
        } else {
          diago = seq(1, S*S, by = S+1)
          Nij.temp = as.vector(Nij)[-diago]
          if( length(which(Nij.temp==0))!=0) # Presence of all transitions
          {
            print("Warning: The transitions are not all observed")
            cat("\n")
            
            p<-Nij/(as.vector(Ni))
            
            f <- Nijk/array(Nij,c(S,S,Kmax))
            f[which(is.na(f))]<-0
            
            q<-Nijk/array(Ni,c(S,S,Kmax))                    
          } else {
            p<-Nij/(as.vector(Ni))
            f <- Nijk/array(Nij,c(S,S,Kmax))
            q<-Nijk/array(Ni,c(S,S,Kmax))
          }
        }
        # return(list(p,f,q))
      }
      
      if( type == 3 ){   
        # Nij<-rowSums(Nijk,dims=2)
        # Ni<-rowSums(Nij)
        Njk<-vector("list", Kmax)
        N<-vector("list", Kmax)
        list.Nijk<-lapply(seq(dim(Nijk)[3]), function(x) Nijk[ , , x])
        for (k in 1:Kmax){
          N[[k]]<-apply(list.Nijk[[k]],2,sum)
          Njk[[k]]<-matrix(unlist(N[[k]]),4,4)
          diag(Njk[[k]])<-0
        }
        Njk<-array(unlist(Njk), c(4,4,Kmax))
        Nj<-colSums(Nij)
        ## Estimators
        
        
        if(length(which(Ni==0))!=0){ # Presence of all states
          
          print("Warning: The transitions are not all observed")
          cat("\n")
          
          p<-Nij/(as.vector(Ni))
          p[which(is.na(p))]<-0
          
          f<-Njk/Nj
          f<-lapply(seq(dim(f)[3]), function(x) f[ , , x])
          Lf<-vector("list",5)
          for (k in 1:Kmax){
            Lf[[k]]<-matrix(f[[k]],4,4)
            f[[k]]<-t(Lf[[k]])
          }
          f<-array(unlist(f), c(4,4,Kmax))
          f[which(is.na(f))]<-0
          
          q<-Nijk/array(Ni,c(S,S,Kmax))
          q[which(is.na(q))]<-0
          
        } else {
          diago = seq(1, S*S, by = S+1)
          Nij.temp = as.vector(Nij)[-diago]
          if( length(which(Nij.temp==0))!=0) # Presence of all transitions
          {
            print("Warning: The transitions are not all observed")
            cat("\n")
            
            p<-Nij/(as.vector(Ni))
            
            f<-Njk/Nj
            f<-lapply(seq(dim(f)[3]), function(x) f[ , , x])
            Lf<-vector("list",5)
            for (k in 1:Kmax){
              Lf[[k]]<-matrix(f[[k]],4,4)
              f[[k]]<-t(Lf[[k]])
            }
            f<-array(unlist(f), c(4,4,Kmax))
            f[which(is.na(f))]<-0
            
            q<-Nijk/array(Ni,c(S,S,Kmax))                    
          } else {
            p<-Nij/(as.vector(Ni))
            f<-Njk/Nj
            f<-lapply(seq(dim(f)[3]), function(x) f[ , , x])
            Lf<-vector("list",5)
            for (k in 1:Kmax){
              Lf[[k]]<-matrix(f[[k]],4,4)
              f[[k]]<-t(Lf[[k]])
            }
            f<-array(unlist(f), c(4,4,Kmax))
            q<-Nijk/array(Ni,c(S,S,Kmax))
          }
          
        }
        
    
        
      }
      
    }
    
  }
  
  ## Inital law
  if(nbSeq > 1){
    init1 = Nstart/sum(Nstart)
  }else{
    ## Computation of the limit distribution
    init1 = .limit.law(q = q, pij = p)
  }
  
  f[which(is.na(f))]<-0
  return(list(init = init1, Ptrans = p, laws = f, q = q))
  
}

## __________________________________________________________
##
## .estimNP_CensureFin
##
## __________________________________________________________
##
## Estimation with the associated MC (Y,U) for several trajectories
.estimNP_CensureFin<-function(file = NULL, seq, E, TypeSojournTime){

  
  if(is.null(file) && is.null(seq)){
    stop("One of the two parameters file or seq must be given")
  }
  
  if(!is.null(file)){
    fasta = read.fasta(file = file)
    seq = getSequence(fasta)
    attr = getAnnot(fasta)
  }
  
  
  if(!is.list(seq)){
    stop("The parameter seq must be a list")
  }
  
  ## States
  TYPES<-c("fij", "fi", "fj", "f")
  type<-pmatch(TypeSojournTime,TYPES)
  
  ## Size of the state space
  S<-length(E)
  
  ## All sequences 
  nbSeq<-length(seq)
  KmaxStart = 0
  Kmax = 0
  Y<-list()
  U<-list()
  for (k in 1:nbSeq){
    ## Manipulations
    JT<-.donneJT(unlist(seq[[k]]))
    J[[k]]<-.code(JT[1,],E)
    T[[k]]<-as.numeric(JT[2,])
    L[[k]]<-T[[k]] - c(1,T[[k]][-length(T[[k]]-1)]) ## sojourn times
    KmaxStart <- max(KmaxStart, max(L[[k]])) 
    YU<-.donneYU(JT, E)
    Y[[k]]<-YU[1,]
    U[[k]]<-YU[2,]
    Kmax <- max(Kmax, max(U[[k]])+1) 
  }
  
  ## get the counts in first position for each state 
  res<-.comptage(J,L,S,KmaxStart)
  Nstart = res$Nstarti
  
  ## Computation of Niujv
  ## with array ?
  Niujv = .calcul.Niujv(Y,U,S,Kmax)
  
  #   Niujv.mat = .calcul.Niujv.matrice(Y,U,S,Kmax)
  
  #   Niu.mat<-rowSums(Niujv.mat)
  Niu = .calcul.Niu(Kmax,S,Niujv)
  
  #   p.mat<-Niujv.mat/(as.vector(Niu.mat))
  #   p.mat[which(is.na(p.mat))]<-0
  p = .p.hat(Kmax = Kmax, S = S, Niujv = Niujv, Niu = Niu)
  
  ##computation of qs  
  ## q.old = .calcul_NP.q.matrice(Y,U,p.mat,S,Kmax)
  q = .calcul_NP.q(S = S, Kmax = Kmax, p = p)
  
  
  pij = rowSums(q,dims=2)

  # fij
  if( type == 1 ){
    f=q/array(pij, c(S,S,Kmax))
    f[which(is.na(f))]<-0   
    
    # fi
  }else if( type == 2 ){
    # sum over the raws of each matrix of the array
    a = apply(q, c(1,3), sum)
    # transformation of the matrix to an array
    f = array(apply(a,2,function(x) matrix(x, nrow = S, ncol = S)), c(S,S,Kmax))
    # fj
  }else if( type == 3 ){
    a = apply(q, c(2,3),sum)
    sum_qi = array(0, c(S,S,Kmax))
    sum_qi = array(apply(a,2,function(x) matrix(x, nrow = S, ncol = S)), c(S,S,Kmax))
    mat_p = colSums(pij)
    calcul.f = array(apply(sum_qi,c(2,3),function(x) x/mat_p), c(S,S,Kmax))
    calcul.f[which(is.na(calcul.f))]<-0
    f = array(apply(calcul.f, 3, function(x) t(x)), c(S,S,Kmax))
    #     for (i in 1:Kmax){
    #       for( j in 1:S ){
    #         sum_qi[,,i] = matrix(colSums(q)[,i])
    #         sum_qi[,,i] = t(sum_qi[,,i])
    #         f[,j,i]  = sum_qi[,j,i]/colSums(p)[j]
    #       }
    #     }
    
    # f
  }else if ( type == 4 ){
    ## warning S = 3 (The states are not all observed)
    sum_qij = apply(q,3,sum)/S
    f = array(apply(as.matrix(sum_qij),1,function(x) matrix(x, nrow = S, ncol = S)),
              c(S,S,Kmax))
    
  }else{
    stop("The parameter \"TypeSojournTime\" is not valid")
  }
  
  ##Initial law
  if(nbSeq > 1){
    init = Nstart/sum(Nstart)
  }else{
    ## Computation of limit law
    init = .limit.law(q = q, pij = pij)
  }
  
  return (list(init = init, Ptrans = pij, laws = f, q = q))
  
}

## __________________________________________________________
##
## .estim.plusTraj
##
## __________________________________________________________
##
####################
##### INPUT: 
## file: path to the file
## seq: list of sequences
## E: state space
## TypeSojournTime: Type of Sojour time ("fij","fi","fj","f")
## distr:  - matrix of distributions if fij
##         - vector of distributions if  fi or fj
## cens.end = 1 if censoring at the end
##            0 if not
## cens.beg = 1 if censoring at the beginning 
##            0 if not
#####################
##### OUTPUT : 
## p: transition matrix
## param:  - array (in param[,,1] the first parameter of the distribution and in param[,,2] the second) if fij
##         - matrix (in param[,1] the first parameter of the distribution and in param[,2] the second) if fi ou fj
##         - vector (in param[1] the first parameter of the distribution and in param[2] the second) if f
.estim.plusTraj<-function(seq, E, TypeSojournTime = "fij", distr = "NP", cens.end = 0, cens.beg = 0){
  
  # require(seqinr)
  # require(DiscreteWeibull)
  # 

  
  if(!is.list(seq)){
    stop("the parameter seq should be a list")
  }
  # 
  ## States
  TYPES<-c("fij", "fi", "fj", "f")
  type<-pmatch(TypeSojournTime,TYPES)
  
  ## Length of state space
  S<-length(E)
  
  ## All sequences 
  nbSeq<-length(seq)
  Kmax = 0
  J<-list()
  T<-list()
  L<-list()
  M<-list()
  
  for (k in 1:nbSeq){
    ## Manipulations
    JT<-.donneJT(unlist(seq[[k]]))
    J[[k]]<-.code(JT[1,],E)
    T[[k]]<-as.numeric(JT[2,])
    L[[k]]<-T[[k]] - c(1,T[[k]][-length(T[[k]]-1)]) ## sojourn times
    Kmax <- max(Kmax, max(L[[k]])) 
  }
  
  
  ## J: list
  ## L: list
  ## S: alphabet size
  ## Kmax: maximal sojourn time
  ################
  ## get the counts
  res<-.comptage(J,L,S,Kmax)
  
  Nij = res$Nij
  Ni = res$Ni
  Nj = res$Nj
  N= res$N
  Nk = res$Nk
  Nik = res$Nik
  Njk = res$Njk
  Nijk = res$Nijk
  Nstart = res$Nstarti
  Nbij = res$Nbij
  Nei = res$Nei
  Nbj = res$Nbj
  Neb = res$Neb
  Nebi = res$Nebi
  Nbi = res$Nbi
  Nb = res$Nb
  Ne = res$Ne
  #########
  
  
  ## Initialization of P
  P <- matrix(0,S,S)
  
  ## fij
  if ( type == 1 ){
    
    if ( !is.matrix(distr) ){
      stop("The parameter distr should be a matrix")
    }
    
    param = array(0, c(S,S,2))
    ## Warning unif infinite ?
    for (i in 1:S){
      for ( j in 1:S ){
        
        if(i != j){
           if ( distr[i,j] == "unif" ){
             
            ## No censoring
            if( cens.end == 0 && cens.beg == 0 ){
              P = Nij/Ni
              
              maxUnif = which.max(Nijk[i,j,])
              param[i,j,1] = maxUnif
              
              ## Censoring at the end
            }else if ( cens.end == 1 && cens.beg == 0 ){
              ## Particular case S = 2 because index 0 for P
              if ( S == 2 ){
                #############################
                ## Estimation
                #############################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ## the Nijk[i,j,] correspond for the vector Nijk to vect.Nijk[i+S*(j-1)+plus]
                ## estimation of p_{jt+1v} et theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute the log-likeihood
                Nij1 = Nij.vect[1]
                Nij2 = Nij.vect[2]
                lois = distr[i,]
                lois = lois[-i]
                
                ## Function to estimate
                logvrais = function(par){
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  u = NULL
                  
                  if("unif" %in% lois){
                    u = 2*which(lois == "unif")-1
                    u = c(u, u+1)
                  }
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw,
                                 dunif(1:Kmax, min = 0, max = par[u[1]])), nrow = Kmax)
                  
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  } 
                  if("unif" %in% lois){    
                    nr = which(lois == "unif")
                    fuv[nr,] = MGM[,5]
                  }
                  
                  puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                  puv.vect0 = as.vector(t(puv))
                  puv.vect = puv.vect0[-El.diag]
                  
                  -( sum(Nij.vect*log(puv.vect)) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dunif(1:Kmax, min = 0, max = par[u[1]], log = TRUE))
                     + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                  )  
                  
                  
                }
                
                ### Vector for the initialisation of optim
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                if("unif" %in% lois){ 
                  u = 2*which(lois == "unif")-1
                  u = c(u, u+1)
                  init[u] = c(Kmax, Kmax)
                }
                
                
                ### Likelihood Optimisation
                CO2<-optim(par = init,fn = logvrais, method ="Nelder-Mead")
                ### Get parameters
                parC = CO2$par
                maxUnif = parC[u[1]]
                param[i,j,1] = maxUnif
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                P = puv
                
              }else{
                #############################
                ## Estimation
                #############################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk est le vecteur de l'array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] corresponds to the Nijk vector for vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} et theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute log-likeihood
                Nij1 = c()
                Nij2 = c()
                a = 1
                while (a < length(Nij.vect)){
                  Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## corresponds to the (S-2) first values of Nij on each line
                  Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## corresponds to the last value of Nij on each line
                  a = a + (S-1)
                }
                
                lois = distr[i,]
                lois = lois[-i]
                
                ## Function to estimate
                logvrais = function(par){
                  m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  p = cbind(m, 1-rowSums(m))
                  puv = p[i,]
                  # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  u = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  if("unif" %in% lois){
                    u = 2*which(lois == "unif")-1
                    u = c(u, u+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw, 
                                 dunif(1:Kmax, min = 0, max = par[S*(S-2)+u[1]])), nrow = Kmax)
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  if("unif" %in% lois){    
                    nr = which(lois == "unif")
                    fuv[nr,] = MGM[,5]
                  }
                  
                  
                  
                  -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dunif(1:(Kmax), min = 0, max = par[(S*(S-2)+u[1])], log = TRUE))
                     + sum( Nei[i,]* log(1-puv%*%fuv))
                  )  
                  
                  
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                diag(u0) <- 1
                c0 = rep(0,(S*(S-2))+2*(S-1))
                
                # puv < 1
                u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                diag(u1) <- -1
                c1 = c(rep(-1,S*(S-2)))
                
                u2 = rbind(u0,u1)
                c2 = c(c0,c1)
           
                
                ###  Vector for the initialisation of optim
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                if("unif" %in% lois){ 
                  u = 2*which(lois == "unif")-1
                  u = c(u, u+1)
                  init[u] = c(Kmax, Kmax)
                }
                
                
                ### log-likelihood optimisation
                CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
                ### get the parameters
                parC = CO2$par
                parL = parC[1:(S*(S-2))]
                maxUnif = parC[(S*(S-2)+u[1])]
                
                ### Transformation of parameters
                ## Add parameters of type: pat = 1 - pac - pag
                ## in par vector of optim
                b = 1
                parP = c()
                while( b <= length(parL) ){
                  parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                  b = b +S-2
                  
                }
                
                ## Add zero in diagonal 
                c = 1 
                parP0 = c()
                while( c < length(parP) ){
                  parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                  c = c + S
                }
                parP0 = c(parP0, 0)
                P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = maxUnif
                
              }
              ## censoring at the beginning
            }else if ( cens.end == 0 && cens.beg == 1 ){
              ########################
              ## Estimation of transition matrix
              #######################
              P = Nij/Ni
              
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              ## In the same way, we create the vector  vect.Nbij : 
              ## vect.Nbij is the vector of the array Nbij
              vect.Nbij = as.vector((Nbij))
              ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
              ########################
              ## Estimation of parameters
              ########################
              logvrais_theta = function(par){
                -(sum(vect.Nijk[i+S*(j-1)+plus]*dunif(1:(Kmax), min = 0, max = par, log = TRUE))
                  + sum(vect.Nbij[i+S*(j-1)+plus]*log(punif(1:Kmax-1, min = 0, max = par, lower.tail = FALSE)))
                )
              }
              
              out<-optim(par = Kmax, fn = logvrais_theta, method = "Brent", lower = 0, upper = Kmax+1)
              param[i,j,1] = out$par
              
              
              ## censoring at the end and at the beginning
            }else if (cens.end == 1 && cens.beg == 1 ){
              if ( S == 2 ){
                #######################################
                ## Estimation of parameters  ##########
                #######################################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} et theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                ## vect.Nbij is the vector of the array Nbij
                vect.Nbij = as.vector((Nbij))
                ## the Nbij[i,j,] corresponds for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
                
                ### Separation of parameters in order to compute log-likelihood
                Nij1 = Nij.vect[1]
                Nij2 = Nij.vect[2]
                lois = distr[i,]
                lois = lois[-i]
                
                ## Function to estimate
                logvrais = function(par){
                  # m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  # p = cbind(m, 1-rowSums(m))
                  # puv = p[i,]
                  # # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  u = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  if("unif" %in% lois){
                    u = 2*which(lois == "unif")-1
                    u = c(u, u+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw,
                                 dunif(1:Kmax, min = 0, max = par[(S*(S-2)+u[1])])), nrow = Kmax)
                  
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  if("unif" %in% lois){    
                    nr = which(lois == "unif")
                    fuv[nr,] = MGM[,5]
                  }
                  
                  puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                  puv.vect0 = as.vector(t(puv))
                  puv.vect = puv.vect0[-El.diag]
                  
                  -( sum(Nij.vect*log(puv.vect))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dunif(1:(Kmax), min = 0, max = par[(S*(S-2)+u[1])], log = TRUE))
                     + sum(vect.Nbij[i+S*(j-1)+plus]*log(punif(1:Kmax-1, min = 0, max = par[(S*(S-2)+u[1])], lower.tail = FALSE)))
                     + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                  )
                }
                
              
                p = c(0,0)
                g = c(0,0)
                nb = c(0,0)
                dw = c(0,0)
                u = c(0,0)
                
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                if("unif" %in% lois){ 
                  u = 2*which(lois == "unif")-1
                  u = c(u, u+1)
                  init[u] = c(Kmax, Kmax+1)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-optim(par=init,logvrais, method ="Nelder-Mead")
                ### Get the parameterss
                parC = CO2$par
                # parL = parC[1:(S*(S-2))]
                theta = parC[u[1]]
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                P = puv
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
                
                
              }else{
                #######################################
                ##Estimation of parameters ##########
                #######################################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## In the same way, we create the vector  vect.Nbij : 
                ## vect.Nbij is the vector of the array Nbij
                vect.Nbij = as.vector((Nbij))
                ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute the log-likelihood
                Nij1 = c()
                Nij2 = c()
                a = 1
                while (a < length(Nij.vect)){
                  Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## corresponds to the (S-2) first values of Nij on each line
                  Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## corresponds to the last value of Nij on each line
                  a = a + (S-1)
                }
                
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  p = cbind(m, 1-rowSums(m))
                  puv = p[i,]
                  
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  u = NULL
                  
                  # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  if("unif" %in% lois){
                    u = 2*which(lois == "unif")-1
                    u = c(u, u+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw,
                                 dunif(1:Kmax, min = 0, max = par[S*(S-2)+u[1]])), nrow = Kmax)
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  if("unif" %in% lois){    
                    nr = which(lois == "unif")
                    fuv[nr,] = MGM[,5]
                  }
                  
                  
                  
                  -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dunif(1:(Kmax), min = 0, max = par[(S*(S-2)+u[1])], log = TRUE))
                     + sum(vect.Nbij[i+S*(j-1)+plus]*log(punif(1:Kmax-1, min = 0, max = par[(S*(S-2)+u[1])], lower.tail = FALSE)))
                     + sum( Nei[i,]* log(1-puv%*%fuv))
                  )
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                diag(u0) <- 1
                c0 = rep(0,(S*(S-2))+2*(S-1))
                
                # puv < 1
                u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                diag(u1) <- -1
                c1 = c(rep(-1,S*(S-2)))
                
                u2 = rbind(u0,u1)
                c2 = c(c0,c1)

                p = c(0,0)
                g = c(0,0)
                nb = c(0,0)
                dw = c(0,0)
                u = c(0,0)
                
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                if("unif" %in% lois){ 
                  u = 2*which(lois == "unif")-1
                  u = c(u, u+1)
                  init[u] = c(Kmax, Kmax + 1)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                parL = parC[1:(S*(S-2))]
                theta = parC[(S*(S-2)+u[1])]
                
                ### Transformation of parameters
                ## Add parameters of type  : pat = 1 - pac - pag
                ## in par vector of optim
                b = 1
                parP = c()
                while( b <= length(parL) ){
                  parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                  b = b +S-2
                  
                }
                
               ## Add zeros on diagonal
                c = 1 
                parP0 = c()
                while( c < length(parP) ){
                  parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                  c = c + S
                }
                parP0 = c(parP0, 0)
                P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
                
              }
              
              
            }
          }else if ( distr[i,j] == "pois" ){
            ## censoring at the end
            if( cens.end == 1 && cens.beg == 0){
              
              ## Particular case S = 2 
               if ( S == 2 ){
                #############################
                ## Estimation of parameters 
                #############################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} et theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute log-likelihood
                Nij1 = Nij.vect[1]
                Nij2 = Nij.vect[2]
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  
                  puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                  puv.vect0 = as.vector(t(puv))
                  puv.vect = puv.vect0[-El.diag]
                  
                  -( sum(Nij.vect*log(puv.vect)) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par[p[1]], log = TRUE))
                     + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                  )  
                  
                  
                }
                
                
                ### Vector to initialise in optim
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-optim(par = init,fn = logvrais, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                theta = parC[p[1]]
                param[i,j,1] = theta
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                P = puv
                
              }else{
                #############################
                ## Estimation of parameters 
                #############################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} et theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute the log-likelihood
                Nij1 = c()
                Nij2 = c()
                a = 1
                while (a < length(Nij.vect)){
                  Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## correspond to the (S-2) first values of Nij on each line
                  Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## correspond to the last value of Nij on each line
                  a = a + (S-1)
                }
                
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  p = cbind(m, 1-rowSums(m))
                  puv = p[i,]
                  # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  
                  
                  
                  -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])], log = TRUE))
                     + sum( Nei[i,]* log(1-puv%*%fuv))
                  )  
                  
                  
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                diag(u0) <- 1
                c0 = rep(0,(S*(S-2))+2*(S-1))
                
                # puv < 1
                u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                diag(u1) <- -1
                c1 = c(rep(-1,S*(S-2)))
                
                u2 = rbind(u0,u1)
                c2 = c(c0,c1)
                
            
                ### Vector to initialise par in optim
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                parL = parC[1:(S*(S-2))]
                theta = parC[(S*(S-2)+p[1])]
                              
                ### Transformation of parameters
                ## Add parameters of type: pat = 1 - pac - pag
                ## in par vector of optim
                b = 1
                parP = c()
                while( b <= length(parL) ){
                  parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                  b = b +S-2
                  
                }
                
               ## Add zeros in diagonale 
                c = 1 
                parP0 = c()
                while( c < length(parP) ){
                  parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                  c = c + S
                }
                parP0 = c(parP0, 0)
                P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
            
              }    
                ##  censoring at the beginning 
            }else if( cens.end == 0 && cens.beg == 1 ){   
              
              ########################
              ## Estimation of the transition matrix
              #######################
              P = Nij/Ni
              
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] corresponds for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              ## In the same way, we create the vector  vect.Nbij : 
              ## vect.Nbij is the vector of the array Nbij
              vect.Nbij = as.vector((Nbij))
              ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
              ########################
              ## Estimation of parameters  of the distribution 
              ########################
              logvrais_theta = function(par){
                -(sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par, log = TRUE))
                  + sum(vect.Nbij[i+S*(j-1)+plus]*log(ppois(1:Kmax-1, lambda = par, lower.tail = FALSE)))
                )
              }
              
              out<-optim(par = 1, fn = logvrais_theta, method = "Brent", lower = 0, upper = Kmax)
              param[i,j,1] = out$par
              
              ## censoring at the end and at the beginning
            }else if( cens.end == 1 && cens.beg == 1 ){
              
              if ( S == 2 ){
                #######################################
                ##Estimation of parameters ##########
                #######################################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} and theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                ## vect.Nbij is the vector of the array Nbij
                vect.Nbij = as.vector((Nbij))
                ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
                
                
                ### Separation of parameters in order to compute the log-likelihood
                Nij1 = Nij.vect[1]
                Nij2 = Nij.vect[2]
                # a = 1
                # while (a < length(Nij.vect)){
                #   Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## correspond to the (S-2) first values of Nij on each line
                #   Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## correspond to the last value of Nij on each line
                #   a = a + (S-1)
                # }
                # 
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  # m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  # p = cbind(m, 1-rowSums(m))
                  # puv = p[i,]
                  # # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  
                  puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                  puv.vect0 = as.vector(t(puv))
                  puv.vect = puv.vect0[-El.diag]
                  
                  -( sum(Nij.vect*log(puv.vect))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])], log = TRUE))
                     + sum(vect.Nbij[i+S*(j-1)+plus]*log(ppois(1:Kmax-1, lambda = par[(S*(S-2)+p[1])], lower.tail = FALSE)))
                     + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                  )
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                # u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                # diag(u0) <- 1
                # c0 = rep(0,(S*(S-2))+2*(S-1))
                # 
                # # puv < 1
                # u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                # diag(u1) <- -1
                # c1 = c(rep(-1,S*(S-2)))
                # 
                # u2 = rbind(u0,u1)
                # c2 = c(c0,c1)
                # 
                ## verification of constraints
                # k2 <- length(c2)
                # n <- dim(u2)[2]
                # for(i in seq_len(k2)) {
                #   j <- which( u2[i,] != 0 )
                #   cat(paste( u2[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
                #   cat(" >= " )
                #   cat( c2[i], "\n" )
                # }
                p = c(0,0)
                g = c(0,0)
                nb = c(0,0)
                dw = c(0,0)
                
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-optim(par=init,logvrais, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                # parL = parC[1:(S*(S-2))]
                theta = parC[p[1]]
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                P = puv
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
                
                
              }else{
                #######################################
                ##Estimation of parameters ##########
                #######################################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## In the same way, we create the vector  vect.Nbij : 
                ## vect.Nbij is the vector of the array Nbij
                vect.Nbij = as.vector((Nbij))
                ##  Nbij[i,j,] corresponds for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute log-likelihood
                Nij1 = c()
                Nij2 = c()
                a = 1
                while (a < length(Nij.vect)){
                  Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## correspond to the (S-2) first values of Nij on each line
                  Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## correspond to the last value of Nij on each line
                  a = a + (S-1)
                }
                
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  p = cbind(m, 1-rowSums(m))
                  puv = p[i,]
                  
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  
                  
                  
                  -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])], log = TRUE))
                     + sum(vect.Nbij[i+S*(j-1)+plus]*log(ppois(1:Kmax-1, lambda = par[(S*(S-2)+p[1])], lower.tail = FALSE)))
                     + sum( Nei[i,]* log(1-puv%*%fuv))
                  )
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                diag(u0) <- 1
                c0 = rep(0,(S*(S-2))+2*(S-1))
                
                # puv < 1
                u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                diag(u1) <- -1
                c1 = c(rep(-1,S*(S-2)))
                
                u2 = rbind(u0,u1)
                c2 = c(c0,c1)
                
                ## verification of constraints
                # k2 <- length(c2)
                # n <- dim(u2)[2]
                # for(i in seq_len(k2)) {
                #   j <- which( u2[i,] != 0 )
                #   cat(paste( u2[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
                #   cat(" >= " )
                #   cat( c2[i], "\n" )
                # }
                p = c(0,0)
                g = c(0,0)
                nb = c(0,0)
                dw = c(0,0)
                
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                parL = parC[1:(S*(S-2))]
                theta = parC[(S*(S-2)+p[1])]
                
                ### Transformation of parameters
                ## Add parameters of type  : pat = 1 - pac - pag
                ## in par vector of optim
                b = 1
                parP = c()
                while( b <= length(parL) ){
                  parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                  b = b +S-2
                  
                }
                
               ## Add zeros in diagonale 
                c = 1 
                parP0 = c()
                while( c < length(parP) ){
                  parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                  c = c + S
                }
                parP0 = c(parP0, 0)
                P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
                
              }
                            
              ## No censoring  
            }else{        
              
              ##############
              ## Estimation of transition matrix
              ##############
              P = Nij/Ni
              P[which(is.na(P))] = 0
              
              ##############
              ## Estimation of theta : 
              #############
              ## vector plus for getting all the Nij for each k
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] corresponds for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              logvrais_theta = function(par){
                -(sum(vect.Nijk[i+S*(j-1)+plus]*dpois(0:(Kmax-1), lambda = par, log = TRUE)))
              }
              
              if (i != j ){ 
                ## Warning for warnings due to par initialization
                # out<-optim(par = 1, logvrais_theta, method = "BFGS")
                # when maximization is posssible withot optim, do without
                x = rep(0:(Kmax-1), vect.Nijk[i+S*(j-1)+plus])
                param[i,j,1] = (sum(x)/sum(vect.Nijk[i+S*(j-1)+plus]))
                # param[i,j,1] = out$par
              }
              
            } 
            
            
            
          }else if ( distr[i,j] == "geom" ){
            
            ## censoring at the end
            if( cens.end == 1 && cens.beg == 0){
              if (S == 2){
                #############################
                ## Estimation of parameters 
                #############################
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} and theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                
                ### Separation of parameters in order to compute the log-likelihood
                Nij1 = Nij.vect[1]
                Nij2 = Nij.vect[2]
                # a = 1
                # while (a < length(Nij.vect)){
                #   Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## corresponds to the (S-2) first values of Nij on each line
                #   Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## corresponds to the last value of Nij on each line
                #   a = a + (S-1)
                # }
                # 
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  # m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  # p = cbind(m, 1-rowSums(m))
                  # puv = p[i,]
                  # # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  
                  puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                  puv.vect0 = as.vector(t(puv))
                  puv.vect = puv.vect0[-El.diag]
                  
                  -( sum(Nij.vect*log(puv.vect)) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par[g[1]], log = TRUE))
                     + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                  )  
                  
                  
                }
                
                
                ### CONSTRAINTS
                # # puv > 0
                # u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                # diag(u0) <- 1
                # c0 = rep(0,(S*(S-2))+2*(S-1))
                # 
                # # puv < 1
                # u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                # diag(u1) <- -1
                # c1 = c(rep(-1,S*(S-2)))
                # 
                # u2 = rbind(u0,u1)
                # c2 = c(c0,c1)
                
                ## verification of constraints
                # k2 <- length(c2)
                # n <- dim(u2)[2]
                # for(i in seq_len(k2)) {
                #   j <- which( u2[i,] != 0 )
                #   cat(paste( u2[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
                #   cat(" >= " )
                #   cat( c2[i], "\n" )
                # }
                
                
                ### Vector to initialise par in optim
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                
                ### Log-likelihood optimisation
                CO2<-optim(par = init,fn = logvrais, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                theta = parC[g[1]]
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                P = puv
                param[i,j,1] = theta
                
              }else{
                #############################
                ## Estimation of parameters 
                #############################
                # tNijk = aperm(Nijk, c(2,1,3))i
                # Nijk.vect = as.vector(tNijk)
                # n = 1
                # d = 1
                # Diag = c()
                # while ( n < S*S*Kmax ){
                #   Diag = c(Diag, seq(from = n, to = d*S*S, by = S+1))
                #   d=d+1
                #   n=n+S*S
                # }
                # Nijk.vect = Nijk.vect[-Diag]
                plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
                ## vect.Nijk is the vector of the array Nijk
                vect.Nijk = as.vector(Nijk)
                ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
                ## estimation function of  p_{jt+1v} and theta_{jt+1v}
                Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
                El.diag = seq(from= 1, to = S*S, by = S+1)
                Nij.vect = Nij.vect0[-El.diag]
                # Nij0 = Nij0[c(etat.cens),]
                
                # 
                ### Separation of parameters in order to compute log-likelihood
                Nij1 = c()
                Nij2 = c()
                a = 1
                while (a < length(Nij.vect)){
                  Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## corresponds to the (S-2) first values of Nij on each line
                  Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## corresponds to the last value of Nij on each line
                  a = a + (S-1)
                }
                
                lois = distr[i,]
                lois = lois[-i]
                
               ## Functions to estimate
                logvrais = function(par){
                  m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                  p = cbind(m, 1-rowSums(m))
                  puv = p[i,]
                  p = NULL
                  g = NULL
                  nb = NULL
                  dw = NULL
                  
                  # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                  # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                  if("pois" %in% lois){
                    p = 2*which(lois == "pois")-1
                    p = c(p, p+1)
                  }
                  if("geom" %in% lois){
                    g = 2*which(lois == "geom")-1
                    g = c(g, g+1)
                  }
                  if("nbinom" %in% lois){
                    nb = 2*which(lois == "nbinom")-1
                    nb = c(nb, nb+1)
                  }
                  if("dweibull" %in% lois){
                    dw = 2*which(lois == "dweibull")-1
                    dw = c(dw, dw+1)
                  }
                  
                  if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                  
                  MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                                 dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                                 dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                                 ddw), nrow = Kmax)
                  
                  # print(MGM)
                  fuv = matrix(nrow = S-1, ncol = Kmax)
                  if("pois" %in% lois){
                    nr = which(lois == "pois")
                    fuv[nr,] = MGM[,1]
                  }
                  if("geom" %in% lois){
                    nr = which(lois == "geom")
                    fuv[nr,] = MGM[,2]
                  }
                  if("nbinom" %in% lois){
                    nr = which(lois == "nbinom")
                    fuv[nr,] = MGM[,3]
                  }
                  if("dweibull" %in% lois){    
                    nr = which(lois == "dweibull")
                    fuv[nr,] = MGM[,4]
                  }
                  -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                     # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                     +  sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])], log = TRUE))
                     + sum( Nei[i,]* log(1-puv%*%fuv))
                  )
                }
                
                
                ### CONSTRAINTS
                # puv > 0
                u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
                diag(u0) <- 1
                c0 = rep(0,(S*(S-2))+2*(S-1))
                
                # puv < 1
                u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
                diag(u1) <- -1
                c1 = c(rep(-1,S*(S-2)))
                
                u2 = rbind(u0,u1)
                c2 = c(c0,c1)
                
                ## verification of constraints
                # k2 <- length(c2)
                # n <- dim(u2)[2]
                # for(i in seq_len(k2)) {
                #   j <- which( u2[i,] != 0 )
                #   cat(paste( u2[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
                #   cat(" >= " )
                #   cat( c2[i], "\n" )
                # }
                p = c(0,0)
                g = c(0,0)
                nb = c(0,0)
                dw = c(0,0)
                
                init = c()
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                  init[p] = c(2, 1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                  init[g] = c(0.5, 0.4)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                  init[nb] = c(2, 4)
                }
                if("dweibull" %in% lois){ 
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                  init[dw] = c(0.6, 0.8)
                }
                
                ### Log-likelihood optimisation
                CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
                ### Get the parameters
                parC = CO2$par
                parL = parC[1:(S*(S-2))]
                theta = parC[(S*(S-2)+g[1])]
                
                ### Transformation of parameters
                ## Add parameters of type: pat = 1 - pac - pag
                ## in par vector of optim
                b = 1
                parP = c()
                while( b <= length(parL) ){
                  parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                  b = b +S-2
                  
                }
                
               ## Add zeros in diagonale 
                c = 1 
                parP0 = c()
                while( c < length(parP) ){
                  parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                  c = c + S
                }
                parP0 = c(parP0, 0)
                P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
                P[which(is.na(P))] = 0
                # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
                param[i,j,1] = theta
              }
            ##  censoring at the beginning 
          }else if( cens.end == 0 && cens.beg == 1 ){   
            
            ########################
            ## Estimation of the transition matrix
            #######################
            P = Nij/Ni
            P[which(is.na(P))] = 0
            
            plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
            ## vect.Nijk is the vector of the array Nijk
            vect.Nijk = as.vector(Nijk)
            ##  Nijk[i,j,] correspondS for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
            ## In the same way, we create the vector  vect.Nbij : 
            ## vect.Nbij is the vector of the array Nbij
            vect.Nbij = as.vector((Nbij))
            ## the Nbij[i,j,] corresponds for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
            ########################
            ## Estimation of parameters  of the distribution 
            ########################
            logvrais_theta = function(par){
              -(sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par, log = TRUE))
                + sum(vect.Nbij[i+S*(j-1)+plus]*log(pgeom(1:Kmax-1, prob = par, lower.tail = FALSE)))
              )
            }
            
            out<-optim(par = 0.5, logvrais_theta, method = "Brent", lower = 0, upper = 1)
            param[i,j,1] = out$par
            
            ## censoring at the end and at the beginning
          }else if( cens.end == 1 && cens.beg == 1 ){
            if ( S == 2 ){
              #######################################
              ##Estimation of parameters ##########
              #######################################
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              ## estimation function of  p_{jt+1v} and theta_{jt+1v}
              Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
              El.diag = seq(from= 1, to = S*S, by = S+1)
              Nij.vect = Nij.vect0[-El.diag]
              # Nij0 = Nij0[c(etat.cens),]
              ## vect.Nbij is the vector of the array Nbij
              vect.Nbij = as.vector((Nbij))
              ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
              
              
              ### Separation of parameters for log-likelihood estimation
              Nij1 = Nij.vect[1]
              Nij2 = Nij.vect[2]
              # a = 1
              # while (a < length(Nij.vect)){
              #   Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## corresponds to the (S-2) first values of Nij on each line
              #   Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## corresponds to the last value of Nij on each line
              #   a = a + (S-1)
              # }
              # 
              lois = distr[i,]
              lois = lois[-i]
              
             ## Functions to estimate
              logvrais = function(par){
                # m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                # p = cbind(m, 1-rowSums(m))
                # puv = p[i,]
                # # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                p = NULL
                g = NULL
                nb = NULL
                dw = NULL
                
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                }
                if("dweibull" %in% lois){
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                }
                
                if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                
                MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                               dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                               dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                               ddw), nrow = Kmax)
                
                
                # print(MGM)
                fuv = matrix(nrow = S-1, ncol = Kmax)
                if("pois" %in% lois){
                  nr = which(lois == "pois")
                  fuv[nr,] = MGM[,1]
                }
                if("geom" %in% lois){
                  nr = which(lois == "geom")
                  fuv[nr,] = MGM[,2]
                }
                if("nbinom" %in% lois){
                  nr = which(lois == "nbinom")
                  fuv[nr,] = MGM[,3]
                }
                if("dweibull" %in% lois){    
                  nr = which(lois == "dweibull")
                  fuv[nr,] = MGM[,4]
                }
                
                puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
                puv.vect0 = as.vector(t(puv))
                puv.vect = puv.vect0[-El.diag]
                
                -( sum(Nij.vect*log(puv.vect))
                   # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                   +  sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par[g[1]], log = TRUE))
                   + sum(vect.Nbij[i+S*(j-1)+plus]*log(pgeom(1:Kmax-1, prob = par[g[1]], lower.tail = FALSE)))
                   + sum( Nei[i,]* log(1-puv.vect%*%fuv))
                )
              }
              
              
              ### CONSTRAINTS
               p = c(0,0)
              g = c(0,0)
              nb = c(0,0)
              dw = c(0,0)
              
              init = c()
              if("pois" %in% lois){
                p = 2*which(lois == "pois")-1
                p = c(p, p+1)
                init[p] = c(2, 1)
              }
              if("geom" %in% lois){
                g = 2*which(lois == "geom")-1
                g = c(g, g+1)
                init[g] = c(0.5, 0.4)
              }
              if("nbinom" %in% lois){
                nb = 2*which(lois == "nbinom")-1
                nb = c(nb, nb+1)
                init[nb] = c(2, 4)
              }
              if("dweibull" %in% lois){ 
                dw = 2*which(lois == "dweibull")-1
                dw = c(dw, dw+1)
                init[dw] = c(0.6, 0.8)
              }
              
              
              ### Log-likelihood optimisation
              CO2<-optim(par=init,logvrais, method ="Nelder-Mead")
              ### Get the parameters
              parC = CO2$par
              # parL = parC[1:(S*(S-2))]
              theta = parC[g[1]]
              
              puv = matrix(c(0,1,1,0), nrow = S, byrow = TRUE)
              P = puv
              param[i,j,1] = theta
              
              
            }else{
              #############################
              ## Estimation of parameters 
              #############################
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              ## In the same way, we create the vector  vect.Nbij : 
              ## vect.Nbij is the vector of the array Nbij
              vect.Nbij = as.vector((Nbij))
              ## the Nbij[i,j,] correspond for the vector Nbij to vect.Nbij[i+S*(j-1)+plus]
              Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
              El.diag = seq(from= 1, to = S*S, by = S+1)
              Nij.vect = Nij.vect0[-El.diag]
              # Nij0 = Nij0[c(etat.cens),]

                
              ### Separation of parameters for log-likelihood estimation
              Nij1 = c()
              Nij2 = c()
              a = 1
              while (a < length(Nij.vect)){
                Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## correspond to the (S-2) first values of Nij on each line
                Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## correspond to the last value of Nij on each line
                a = a + (S-1)
              }
              
              lois = distr[i,]
              lois = lois[-i]
              
             ## Functions to estimate
              logvrais = function(par){
                m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                p = cbind(m, 1-rowSums(m))
                puv = p[i,]
                p = NULL
                g = NULL
                nb = NULL
                dw = NULL
                
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")
                  p = c(p, p+1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                }
                if("dweibull" %in% lois){
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                }
                
                if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                
                MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[(S*(S-2)+p[1])]), 
                               dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])]),
                               dnbinom(0:(Kmax-1), size = par[(S*(S-2)+nb[1])], mu = par[(S*(S-2)+nb[2])]), 
                               ddw), nrow = Kmax)
                
                # print(MGM)
                fuv = matrix(nrow = S-1, ncol = Kmax)
                if("pois" %in% lois){
                  nr = which(lois == "pois")
                  fuv[nr,] = MGM[,1]
                }
                if("geom" %in% lois){
                  nr = which(lois == "geom")
                  fuv[nr,] = MGM[,2]
                }
                if("nbinom" %in% lois){
                  nr = which(lois == "nbinom")
                  fuv[nr,] = MGM[,3]
                }
                if("dweibull" %in% lois){    
                  nr = which(lois == "dweibull")
                  fuv[nr,] = MGM[,4]
                }
                
                -( sum(Nij1*log(par[1:(S*(S-2))])) + sum(Nij2*log(1-rowSums(m))) #sum(par[1:(S-2)])))
                   # + sum(Nijk.vect * dpois(rep(0:(Kmax-1), each = (S-1)*S), lambda = par[(S*(S-2)+1)], log = TRUE))
                   +  sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par[(S*(S-2)+g[1])], log = TRUE))
                   + sum(vect.Nbij[i+S*(j-1)+plus]*log(pgeom(1:Kmax-1, prob = par[(S*(S-2)+g[1])], lower.tail = FALSE)))
                   + sum( Nei[i,]* log(1-puv%*%fuv))
                )
              }
              
              
              ### CONSTRAINTS
              # puv > 0
              u0 = matrix(0, nrow = (S*(S-2))+2*(S-1), ncol = (S*(S-2))+2*(S-1))
              diag(u0) <- 1
              c0 = rep(0,(S*(S-2))+2*(S-1))
              
              # puv < 1
              u1 = matrix(0, nrow = S*(S-2), ncol = (S*(S-2))+2*(S-1))
              diag(u1) <- -1
              c1 = c(rep(-1,S*(S-2)))
              
              u2 = rbind(u0,u1)
              c2 = c(c0,c1)
              
              ## verification of constraints
              # k2 <- length(c2)
              # n <- dim(u2)[2]
              # for(i in seq_len(k2)) {
              #   j <- which( u2[i,] != 0 )
              #   cat(paste( u2[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
              #   cat(" >= " )
              #   cat( c2[i], "\n" )
              # }
              p = c(0,0)
              g = c(0,0)
              nb = c(0,0)
              dw = c(0,0)
              
              init = c()
              if("pois" %in% lois){
                p = 2*which(lois == "pois")-1
                p = c(p, p+1)
                init[p] = c(2, 1)
              }
              if("geom" %in% lois){
                g = 2*which(lois == "geom")-1
                g = c(g, g+1)
                init[g] = c(0.5, 0.4)
              }
              if("nbinom" %in% lois){
                nb = 2*which(lois == "nbinom")-1
                nb = c(nb, nb+1)
                init[nb] = c(2, 4)
              }
              if("dweibull" %in% lois){ 
                dw = 2*which(lois == "dweibull")-1
                dw = c(dw, dw+1)
                init[dw] = c(0.6, 0.8)
              }
              
              ### Log-likelihood optimisation
              CO2<-constrOptim(theta=c(rep(1/(S-1),S*(S-2)), init),logvrais, ui=u2,ci=c2, method ="Nelder-Mead")
              ### Get the parameters
              parC = CO2$par
              parL = parC[1:(S*(S-2))]
              theta = parC[(S*(S-2)+g[1])]
              
              ### Transformation of parameters
              ## Add parameters of type: pat = 1 - pac - pag
              ## in par vector of optim
              b = 1
              parP = c()
              while( b <= length(parL) ){
                parP = c(parP, parL[b:(b+(S-3))], 1-sum(parL[b:(b+(S-3))]))
                b = b +S-2
                
              }
              
             ## Add zeros in diagonale 
              c = 1 
              parP0 = c()
              while( c < length(parP) ){
                parP0 = c(parP0, 0,parP[c:(c+(S-1))])
                c = c + S
              }
              parP0 = c(parP0, 0)
              P = matrix(parP0, nrow = S, ncol = S, byrow = TRUE)
              P[which(is.na(P))] = 0
              # param = array(c(matrix(theta0, nrow = S, ncol = S, byrow = TRUE),matrix(0, nrow = S, ncol = S)), c(S,S,2))
              param[i,j,1] = theta
            }  
            ## No censoring  
          }else{        
            
            ##############
            ## Estimation of the transition matrix
            ##############
            P = Nij/Ni
            P[which(is.na(P))] = 0
            
            ##############
            ## Estimation of theta: 
            #############
            ## vecteur plus for getting all the Nij for each k
            plus = seq(from = 0, to = Kmax*S*S-1, by = S*S )
            ## vect.Nijk is the vector of the array Nijk
            vect.Nijk = as.vector(Nijk)
            ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
            logvrais_theta = function(par){
              -(sum(vect.Nijk[i+S*(j-1)+plus]*dgeom(0:(Kmax-1), prob = par, log = TRUE)))
            }
            
            if (i != j ){ 
              ## Warning for warnings due to par initialization
              # out<-optim(par = 1, logvrais_theta, method = "BFGS")
              x = rep(1:(Kmax), vect.Nijk[i+S*(j-1)+plus])
              param[i,j,1] = 1/(sum(x)/sum(vect.Nijk[i+S*(j-1)+plus]))
              # param[i,j,1] = out$par
            }
            
          } 
          
          
        }else if ( distr[i,j] == "nbinom"){
          
          ## censoring at the end
          if( cens.end == 1 && cens.beg == 0){
            if (S == 2){
              #############################
              ## Estimation of parameters 
              #############################
              plus = seq(from = 0, to = Kmax*S*S-1, by = S*S)
              ## vect.Nijk is the vector of the array Nijk
              vect.Nijk = as.vector(Nijk)
              ##  Nijk[i,j,] correspond for Nijk vector to vect.Nijk[i+S*(j-1)+plus]
              ## estimation function of  p_{jt+1v} and theta_{jt+1v}
              Nij.vect0<-as.vector(t(Nij)) #[as.vector(t(Nij))!=0]
              El.diag = seq(from= 1, to = S*S, by = S+1)
              Nij.vect = Nij.vect0[-El.diag]
              # Nij0 = Nij0[c(etat.cens),]
              
              
              ### Separation of parameters for log-likelihood estimation
              Nij1 = Nij.vect[1]
              Nij2 = Nij.vect[2]
              # a = 1
              # while (a < length(Nij.vect)){
              #   Nij1 = c(Nij1,Nij.vect[a:(a+(S-3))]) ## correspond to the (S-2) first values of Nij on each line
              #   Nij2 = c(Nij2, Nij.vect[a+(S-2)]) ## correspond to the last value of Nij on each line
              #   a = a + (S-1)
              # }
              # 
              lois = distr[i,]
              lois = lois[-i]
              
             ## Functions to estimate
              logvrais = function(par){
                # m = matrix(par[1:(S*(S-2))], ncol = S-2, byrow = TRUE)
                # p = cbind(m, 1-rowSums(m))
                # puv = p[i,]
                # # ind1 = seq(from = (S*(S-2))+1, to =(S*(S-2)+5), by=2)
                # ind2 = seq(from = (S*(S-2))+2, to =(S*(S-2)+6), by=2)
                p = NULL
                g = NULL
                nb = NULL
                dw = NULL
                
                if("pois" %in% lois){
                  p = 2*which(lois == "pois")-1
                  p = c(p, p+1)
                }
                if("geom" %in% lois){
                  g = 2*which(lois == "geom")-1
                  g = c(g, g+1)
                }
                if("nbinom" %in% lois){
                  nb = 2*which(lois == "nbinom")-1
                  nb = c(nb, nb+1)
                }
                if("dweibull" %in% lois){
                  dw = 2*which(lois == "dweibull")-1
                  dw = c(dw, dw+1)
                }
                
                if(is.null(dw)){ddw = rep(0, Kmax)}else{ddw = ddweibull(1:Kmax, q = par[(S*(S-2)+dw[1])], beta = par[(S*(S-2)+dw[2])])}
                
                MGM = matrix(c(dpois(0:(Kmax-1), lambda = par[p[1]]), 
                               dgeom(0:(Kmax-1), prob = par[g[1]]),
                               dnbinom(0:(Kmax-1), size = par[nb[1]], mu = par[nb[2]]), 
                               ddw), nrow = Kmax)
                
                
                # print(MGM)
                fuv = matrix(nrow = S-1, ncol = Kmax)
                if("pois" %in% lois){
                  nr = which(lois == "pois")
                  fuv[nr,] = MGM[,1]
                }
                if("geom" %in% lois){
                  nr = which(lois == "geom")
                  fuv[nr,] = MGM[,2]
                }
                if("nbinom" %in% lois){
                  nr = which(lois == "nbinom")
                  fuv[nr,] = MGM[,3]
                }