R/InternalTrajeR.R

In trajeR: Group Based Modeling Trajectory

#################################################################################
# find parameters for a CNORM model
#################################################################################
#' Internal function to fit CNORM Model
#'
#' @inheritParams trajeR
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @param ng Integer. The number of groups.
#' @param nx Integer. The number of covariates.
#' @param n Integer. Number of individuals.
#' @param nbeta Vector of integers. Number of beta parameters for each group.
#' @param nw Integer. Number of time dependent covariate.
#' @param ntheta Vector of integers. Number of theta parameters for each group.
#' @param period Integer.
#' @param theta Vector of real. The parameter for calculated the group membership probability.
#' @param beta Vector of real. The beta parameter.
#' @param sigma Vector of real. The sigma parameter.
#' @param delta Vector of real. The delta parameter.
#' @param pi Vector of real. The group membership probability.
#' @param EMIRLS Boolean. True if we use EMIRLS method.
#' @return  return a object of class Trajectory.CNORM
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item sigma - vector of the parameters sigma.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape.
#'   \item min - a real with the minimum value for censored data.
#'   \item max - a real with the maximum value for censored data.
#' }

trajeR.CNORM <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nw, ntheta, period, degre, theta, beta, sigma, delta, pi, Method, ssigma,
                         ymax, ymin, hessian, itermax, paraminit, EMIRLS, refgr){
  set_tour(1)
  set_storelik(10**100)
  nsigma = ng
  theta = theta - theta[1:nx]
  if (Method == "L"){
    theta = theta[-c(1:nx)]
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      sigma = rep(stats::sd(Y, na.rm = TRUE),ng)
      paraminit = c(theta, unlist(beta), sigma, unlist(delta))
    }else{
      paraminit = paraminit[-c(1:nx)]
    }
    if (ssigma == FALSE){
      # different sigma
      paraminitL = c(paraminit[1:((ng -1)*nx+sum(nbeta))], log(paraminit[((ng - 1)*nx+sum(nbeta)+1):((ng - 1)*nx+sum(nbeta)+ng)]),
                     paraminit[-c(1: ((ng - 1)*nx+sum(nbeta)+ng))])
      if (!hessian){
        newparam = stats::optim(par = paraminitL, fn = Likelihoodalpha_cpp, gr=difLalpha_cpp,
                                method = "BFGS",
                                hessian = hessian,
                                control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                                ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta,
                                ymin = ymin, ymax = ymax, nw = nw, TCOV = TCOV, ssigma = ssigma)
      } else {
        newparam = ucminf::ucminf(par = paraminitL, 
                                  fn = LCNORM, 
                                  gr = difLCNORM,
                                  hessian = 2,
                                  control = list(stepmax=10**(-6), maxeval = itermax),
                                  ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta,
                                  ymin = ymin, ymax = ymax, nw = nw, TCOV = TCOV, ssigma = ssigma)     
      }
    } else {
      # same sigma
      paraminitL = c(paraminit[1:((ng -1)*nx+sum(nbeta))], log(paraminit[(ng - 1)*nx+sum(nbeta)+1]),
                     paraminit[-c(1: ((ng - 1)*nx+sum(nbeta)+ng))])
      if (!hessian){
        newparam = stats::optim(par = paraminitL, fn = Likelihoodalpha_cpp, gr=difLalphaunique_cpp,
                                method = "BFGS",
                                hessian = hessian,
                                control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                                ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta,
                                ymin = ymin, ymax = ymax, nw = nw, TCOV = TCOV, ssigma = ssigma)
      }else {
        newparam = ucminf::ucminf(par = paraminitL, 
                                  fn = LCNORM, 
                                  gr = difLCNORMss,
                                  hessian = 2,
                                  control = list(stepmax=10**(-5), maxeval = itermax),
                                  ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta,
                                  ymin = ymin, ymax = ymax, nw = nw, TCOV = TCOV, ssigma = ssigma)
      }
    }
    param = newparam$par
    indtheta = 1:((ng-1)*nx)
    indbeta = ((ng-1)*nx+1):((ng-1)*nx+sum(nbeta))
    if (ssigma){
      indsigma = ((ng-1)*nx+sum(nbeta)+1)
      if (nw!=0){
        inddelta = ((ng-1)*nx+sum(nbeta)+2):((ng-1)*nx+sum(nbeta)+1+nw*ng) 
      }
    }else{
      indsigma = ((ng-1)*nx+sum(nbeta)+1):((ng-1)*nx+sum(nbeta)+ng)
      if (nw!=0){
        inddelta = ((ng-1)*nx+sum(nbeta)+ng+1):((ng-1)*nx+sum(nbeta)+ng+nw*ng)  
      }
    }
    if (ssigma){
      if (nw!=0){
        param = c(param[c(indtheta, indbeta)], rep(exp(param[indsigma]), ng), param[inddelta])  
      }else{
        param = c(param[c(indtheta, indbeta)], rep(exp(param[indsigma]), ng))
      }
    }else{
      sigma = exp(param[indsigma])  
      param[indsigma] = sigma
    }
    
    
    theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
    param = c(param[-c(1:((ng-1)*nx))], theta)
    if (hessian == TRUE){
      invH = newparam$invhessian
      SE = sqrt(diag(invH))
      if (ssigma){
        sdsigma = rep(exp(newparam$par[indsigma[1]])*SE[indsigma[1]], ng)
      }else{
        matsigma = diag(sigma)
        sdsigma = sqrt(diag(matsigma %*%  invH[indsigma, indsigma] %*% matsigma)) 
      }
      if (nx == 1){
        sdtmp = deltaTheta(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        sdbase = deltaThetaBase(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        if (nw!=0){
          SE = c(SE[indbeta], sdsigma, SE[inddelta], sdbase, sdtmp)  
        }else{
          SE = c(SE[indbeta], sdsigma, sqrt(sum(sdtmp**2)), sdtmp)
        }
      }else{
        if (nw!=0){
          SE = c(SE[indbeta], sdsigma, SE[inddelta], rep(NA, nx), SE[indtheta])  
        }else{
          SE = c(SE[indbeta], sdsigma, rep(NA, nx), SE[indtheta])
        }
      }
    }
  } else if (Method == "EM"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      sigma = rep(stats::sd(Y),ng)
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)], unlist(beta), sigma, unlist(delta))
      }else{
        paraminitEM = c(theta, unlist(beta), sigma, unlist(delta))
      }
    }else{
      if (nx == 1){paraminitEM = paraminit[-ng]}
    }
    if (max(Y)<ymax & min(Y)>ymin){
      if (ssigma == FALSE){
        param = EM_cpp(paraminitEM, ng, nx, nbeta, n, A, Y, X, ymin, ymax, TCOV, nw, itermax, EMIRLS, refgr)
      }else{
        param = EMSigmaunique_cpp(paraminitEM, ng, nx, nbeta, n, A, Y, X, ymin, ymax, TCOV, nw, itermax, EMIRLS, refgr)
      }
    }else{
      if (ssigma == TRUE){
        param = EMCensoredSigmaunique_cpp(paraminitEM, ng, nx, nbeta, n, A, Y, X, ymin, ymax, TCOV, nw, itermax, EMIRLS, refgr)
      }
      else{
        param = EMCensored_cpp(paraminitEM, ng, nx, nbeta, n, A, Y, X, ymin, ymax, TCOV, nw, itermax, EMIRLS, refgr)
      }
    }
    if (hessian == TRUE){
      SE = IEM_cpp(param, ng, nx, nbeta, n, A, Y, X, ymin, ymax, TCOV, nw, refgr)
      if (nx == 1){
        if (nw == 0){
          SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))  
        }else{
          SE = c(SE[ng:(ng + sum(nbeta) - 1)], SE[(ng + sum(nbeta)+sum(nw*ng)):(ng + sum(nbeta)+sum(nw*ng) + ng - 1)],
                 SE[(ng + sum(nbeta)):(ng + sum(nbeta)+sum(nw*ng) - 1)], 
                 sqrt(sum(SE[1:(ng-1)]**2)), SE[1:(ng-1)])
        }
      }else{
        if (nw == 0){
          SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])  
        }else{
          SE = c(SE[(ng*nx):(ng*nx + sum(nbeta) - 1)], SE[(ng*nx + sum(nbeta)+sum(nw*ng)):(ng*nx + sum(nbeta)+sum(nw*ng) + ng - 1)],
                 SE[(ng*nx + sum(nbeta)):(ng*nx + sum(nbeta)+sum(nw*ng) - 1)])
        }
      }
    }else{
      SE = NA
    }
    param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
  }
  if (hessian == TRUE){
    if (nx == 1 & Method == "L"){
      paramtmp = c(param[1:(length(param) - ng*nx)], exp(theta)/sum(exp(theta)))
      prob = 1-stats::pt(abs(paramtmp/SE), n*period-1)+stats::pt(-abs(paramtmp/SE), n*period-1) 
    }else{
      prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    }
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = prob
    )
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob = SE
    )
  }
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  namesigma = paste0("sigma", 1:ng)
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namesigma, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  if (nw == 0){
    delta = NA
    sigma = param[(sum(nbeta)+1):(sum(nbeta)+ng)]
    theta = param[-c(1:(sum(nbeta)+ng))]
  }else{
    delta = param[(sum(nbeta)+ng+1):(sum(nbeta)+ng+nw*ng)]
    sigma = param[(sum(nbeta)+1):(sum(nbeta)+ng)]
    theta = param[-c(1:(sum(nbeta)+ng+nw*ng))]
  }
  res = list(beta = beta,
             sigma = sigma,
             delta = delta,
             theta = theta,
             sd = SE, tab = d, Model = "CNORM",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(param = c(theta, beta, sigma, delta), model = "CNORM",
                                     method = Method,
                                     ng = ng, nx = nx, n = n,
                                     nbeta = nbeta, nw = nw, A= A, Y = Y, X = X,
                                     TCOV = TCOV, ymin = ymin, ymax = ymax),
             Time = A[1,], degre = degre - 1, min = ymin, max = ymax)
  class(res) = "Trajectory.CNORM"
  return(res)
}
#################################################################################
# find parameters for a LOGIT model
#################################################################################
#' Internal function to fit LOGIT Model
#'
#' @inheritParams trajeR
#' @inheritParams trajeR.CNORM
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @return  return a object of class Trajectory.LOGIT
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape.
#' }


trajeR.LOGIT <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nw, ntheta, period, degre, theta, beta, delta, pi, Method,
                         hessian, itermax, paraminit, EMIRLS, refgr){
  set_tour(1)
  set_storelik(10**100)
  theta = theta - theta[1:nx]
  theta = theta[-c(1:nx)]
  if (Method == "L"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      paraminit = c(theta, unlist(beta), unlist(delta))
    }
    if (!hessian){
      newparam = stats::optim(par = paraminit, fn = likelihoodLOGIT_cpp, gr = difLLOGIT_cpp,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                              ng = ng, nx = nx, n = n, A = A, Y = Y, X = X, nbeta = nbeta,
                              nw = nw, TCOV = TCOV)
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }else{
      newparam = ucminf::ucminf(par = paraminit, 
                                fn = LLOGIT, 
                                gr = difLLOGIT,
                                hessian = 2,
                                control = list(stepmax=10**(-5), maxeval = itermax),
                                ng = ng, nx = nx, n = n,  A = A, Y = Y, X = X, nbeta = nbeta,
                                nw = nw, TCOV = TCOV)
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }
    if (hessian == TRUE){
      invH = newparam$invhessian
      SE = sqrt(diag(invH))
      if (nx == 1){
        sdtmp = deltaTheta(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        sdbase = deltaThetaBase(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        SE = c(SE[-c(1:((ng-1)*nx))], sdbase, sdtmp)
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(NA, nx), SE[c(1:((ng-1)*nx))])
      }
    }
  } else if (Method == "EM"){
    # intial value for EM
    if (is.null(paraminit)){
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)],unlist(beta), unlist(delta))
      }else{
        paraminitEM = c(theta, unlist(beta), unlist(delta))
      }
    }else{
      if (nx == 1){
        paraminitEM = paraminit[-ng]
      }else{
        paraminitEM = paraminit
      }
    }
    param = EMLOGIT_cpp(paraminitEM, ng, nx, n, nbeta, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
    if (hessian == TRUE){
      SE = IEMLOGIT_cpp(param, ng, nx, nbeta, n, A, Y, X, TCOV, nw, refgr)
      if (nx == 1){
        SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)], sqrt(sum(SE[1:(ng-1)]**2)))
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
      }
    }else{
      SE = NA
    }
    if (nx == 1){
      param = c(param[-c(1:ng)], param[1:ng])
    }else{
      param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
    }
  } else if (Method == "EMIRLS"){
    # intial value for EM
    if (is.null(paraminit)){
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)],unlist(beta), unlist(delta))
      }else{
        paraminitEM = c(theta, unlist(beta), unlist(delta))
      }
    }else{
      if (nx == 1){
        paraminitEM = paraminit[-ng]
      }else{
        paraminitEM = paraminit
      }
    }
    param = EMLOGITIRLS_cpp(paraminitEM, ng, nx, n, nbeta, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
    if (hessian == TRUE){
      SE =  IEMLOGIT_cpp(param, ng, nx, nbeta, n, A, Y, X, TCOV, nw, refgr)
      if (nx == 1){
        SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
      }
    }else{
      SE = NA
    }
    if (nx == 1){
      param = param = c(param[-c(1:ng)], param[1:ng])
    }else{
      param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
    }
  }
  if (hessian == TRUE){
    if (nx == 1 & Method == "L"){
      paramtmp = c(param[1:(length(param) - ng*nx)], exp(theta)/sum(exp(theta)))
      prob = 1-stats::pt(abs(paramtmp/SE), n*period-1)+stats::pt(-abs(paramtmp/SE), n*period-1) 
    }else{
      prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    }
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = prob
    ) 
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob = SE
    )
  }
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  if (nw == 0){
    delta = NA
    theta = param[-c(1:(sum(nbeta)))]
  }else{
    delta = param[(sum(nbeta)+1):(sum(nbeta)+nw*ng)]
    theta = param[-c(1:(sum(nbeta)+nw*ng))]
  }
  res = list(beta = beta,
             delta = delta,
             theta = theta,
             sd = SE, tab = d, Model = "LOGIT",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(c(theta, beta, delta), model = "LOGIT", method = Method,
                                     ng = ng, nx = nx, n = n, nbeta = nbeta, nw = nw,
                                     A = A, Y = Y, X = X, TCOV = TCOV),
             Time = A[1,], degre = degre - 1)
  class(res) = "Trajectory.LOGIT"
  return(res)
}
#################################################################################
# find parameters for a ZIP model
#################################################################################
#' Internal function to fit ZIP Model
#' 
#' @inheritParams trajeR
#' @inheritParams trajeR.CNORM
#' @param nu Vector of real. The nu parameter.
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @return  return a object of class Trajectory.ZIP
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item nu - vector of the parameters nu.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape for the Poisson part.
#'   \item degre.nu - a vector with the degree of the polynomial shape for the exceeded zero state.
#' }

trajeR.ZIP <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nw, ntheta, period, degre, degre.nu, theta, beta, nu, delta, pi, Method,
                       hessian, itermax, paraminit, EMIRLS, refgr){
  set_tour(1)
  set_storelik(10**100)
  theta = theta - theta[1:nx]
  theta = theta[-c(1:nx)]
  degre.nu = degre.nu + 1
  nnu = degre.nu
  #  nu = unlist(c(sapply(1:ng, function(s){
  #    c(min(Y)+s*(max(Y)-min(Y))/(ng+1), rep(0, degre.nu[s]-1))
  #  })))
  if (Method == "L"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      paraminit = c(theta, unlist(beta), unlist(nu), delta)
    }
    if (!hessian){
      newparam = stats::optim(par = paraminit, fn = likelihoodZIP_cpp, gr = difLZIP_cpp,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                              ng = ng, nx = nx, n = n, nnu = nnu, A = A, Y = Y, X = X, nbeta = nbeta,
                              nw = nw, TCOV = TCOV)  
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }else{
      newparam = ucminf::ucminf(par = paraminit, 
                                fn = LZIP, 
                                gr = difLZIP,
                                hessian = 2,
                                control = list(stepmax=10**(-5), maxeval = itermax),
                                ng = ng, nx = nx, n = n, nnu = nnu, A = A, Y = Y, X = X, nbeta = nbeta,
                                nw = nw, TCOV = TCOV)
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }
    if (hessian == TRUE){
      invH = newparam$invhessian
      SE = sqrt(diag(invH))
      if (nx == 1){
        sdtmp = deltaTheta(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        sdbase = deltaThetaBase(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        SE = c(SE[-c(1:((ng-1)*nx))], sdbase, sdtmp)
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(NA, nx), SE[c(1:((ng-1)*nx))])
      }
    }
  } else if (Method == "EM"){
    # intial value for EM
    if (is.null(paraminit)){
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)], unlist(beta), unlist(nu), unlist(delta))
      }else{
        paraminitEM = c(rep(0,nx), theta, unlist(beta), unlist(nu), unlist(delta))
      }
    }else{
      if (nx == 1){
        paraminitEM = paraminit[-ng]
      }else{
        paraminitEM = paraminit
      }
    }
    param = EMZIP_cpp(paraminitEM, ng, nx, n, nbeta, nnu, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
    if (hessian == TRUE){
      SE = IEMZIP_cpp(param, ng, nx, nbeta, nnu, n, A, Y, X, TCOV, nw, refgr)
      if (nx == 1){
        SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
      }
    }else{
      SE = NA
    }
    if (nx == 1){
      param = c(param[-c(1:ng)], param[1:ng])
    }else{
      param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
    }
  }else if (Method == "EMIRLS"){
    # intial value for EM
    if (is.null(paraminit)){
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)], unlist(beta), unlist(nu), unlist(delta))
      }else{
        paraminitEM = c(rep(0,nx), theta, unlist(beta), unlist(nu), unlist(delta))
      }
    }else{
      if (nx == 1){
        paraminitEM = paraminit[-ng]
      }else{
        paraminitEM = paraminit
      }
    }
    param = EMZIPIRLS_cpp(paraminitEM, ng, nx, n, nbeta, nnu, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
    if (hessian == TRUE){
      SE = IEMZIP_cpp(param, ng, nx, nbeta, nnu, n, A, Y, X, TCOV, nw, refgr)
      if (nx == 1){
        SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
      }
    }else{
      SE = NA
    }
    if (nx == 1){
      param = c(param[-c(1:ng)], param[1:ng])
    }else{
      param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
    }
  }
  if (hessian == TRUE){
    if (nx == 1 & Method == "L"){
      paramtmp = c(param[1:(length(param) - ng*nx)], exp(theta)/sum(exp(theta)))
      prob = 1-stats::pt(abs(paramtmp/SE), n*period-1)+stats::pt(-abs(paramtmp/SE), n*period-1) 
    }else{
      prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    }
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = prob
    )
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob = SE)
  }
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  namenu =c()
  for (k in 1:ng){
    namenu = c(namenu, paste0("Nu", k, 1:nnu[k]))
  }
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namenu, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  if (nw == 0){
    delta = NULL
    nu = param[(sum(nbeta)+1):(sum(nbeta)+sum(nnu))]
    theta = param[-c(1:(sum(nbeta)+sum(nnu)))]
  }else{
    nu = param[(sum(nbeta)+1):(sum(nbeta)+sum(nnu))]
    delta = param[(sum(nbeta)+sum(nnu)+1):(sum(nbeta)+sum(nnu)+nw*ng)]
    theta = param[-c(1:(sum(nbeta)++sum(nnu)+nw*ng))]
  }
  method = ifelse(nx != 1, "L", Method)
  res = list(beta = beta,
             delta = delta,
             theta = theta,
             nu = nu,
             sd = SE, tab = d, Model = "ZIP",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(param = c(theta, beta, nu, delta), model = "ZIP", 
                                     method = method, ng =ng, 
                                     nx = nx, n = n, nbeta = nbeta, nw = nw, A = A, Y = Y, X = X, 
                                     TCOV = TCOV, nnu = nnu),
             Time = A[1,], degre = degre - 1, degre.nu = degre.nu - 1)
  class(res) = "Trajectory.ZIP"
  return(res)
}
#################################################################################
# find parameters for a Poisson model
#################################################################################
#' Internal function to fit poisson Model
#' 
#' @inheritParams trajeR
#' @inheritParams trajeR.CNORM
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @return  return a object of class Trajectory.Pois
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape for the Poisson part.
#' }

trajeR.POIS <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nw, ntheta, period, degre, theta, beta, delta, pi, Method,
                        hessian, itermax, paraminit, EMIRLS, refgr){
  set_tour(1)
  set_storelik(10**100)
  theta = theta - theta[1:nx]
  theta = theta[-c(1:nx)]
  #  nu = unlist(c(sapply(1:ng, function(s){
  #    c(min(Y)+s*(max(Y)-min(Y))/(ng+1), rep(0, degre.nu[s]-1))
  #  })))
  if (Method == "L"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      paraminit = c(theta, unlist(beta), delta)
    }
    if (!hessian){
      newparam = stats::optim(par = paraminit, fn = likelihoodPois_cpp, gr = difLPois_cpp,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                              ng = ng, nx = nx, n = n, A = A, Y = Y, X = X, nbeta = nbeta,
                              nw = nw, TCOV = TCOV)  
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }
    # }else{
    #   newparam = ucminf::ucminf(par = paraminit, 
    #                             fn = LZIP, 
    #                             gr = difLZIP,
    #                             hessian = 2,
    #                             control = list(stepmax=10**(-5), maxeval = itermax),
    #                             ng = ng, nx = nx, n = n, nnu = nnu, A = A, Y = Y, X = X, nbeta = nbeta,
    #                             nw = nw, TCOV = TCOV)
    #   param = newparam$par
    #   theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
    #   param = c(param[-c(1:((ng-1)*nx))], theta)
    # }
    #   if (hessian == TRUE){
    #     invH = newparam$invhessian
    #     SE = sqrt(diag(invH))
    #     if (nx == 1){
    #       sdtmp = deltaTheta(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
    #       sdbase = deltaThetaBase(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
    #       SE = c(SE[-c(1:((ng-1)*nx))], sdbase, sdtmp)
    #     }else{
    #       SE = c(SE[-c(1:((ng-1)*nx))], rep(NA, nx), SE[c(1:((ng-1)*nx))])
    #     }
    #   }
  }
  # } else if (Method == "EM"){
  #   # intial value for EM
  #   if (is.null(paraminit)){
  #     if (nx == 1){
  #       paraminitEM = c(pi[1:(ng-1)], unlist(beta), unlist(nu), unlist(delta))
  #     }else{
  #       paraminitEM = c(rep(0,nx), theta, unlist(beta), unlist(nu), unlist(delta))
  #     }
  #   }else{
  #     if (nx == 1){
  #       paraminitEM = paraminit[-ng]
  #     }else{
  #       paraminitEM = paraminit
  #     }
  #   }
  #   param = EMZIP_cpp(paraminitEM, ng, nx, n, nbeta, nnu, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
  #   if (hessian == TRUE){
  #     SE = IEMZIP_cpp(param, ng, nx, nbeta, nnu, n, A, Y, X, TCOV, nw, refgr)
  #     if (nx == 1){
  #       SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))
  #     }else{
  #       SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
  #     }
  #   }else{
  #     SE = NA
  #   }
  #   if (nx == 1){
  #     param = c(param[-c(1:ng)], param[1:ng])
  #   }else{
  #     param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
  #   }
  # }else if (Method == "EMIRLS"){
  #   # intial value for EM
  #   if (is.null(paraminit)){
  #     if (nx == 1){
  #       paraminitEM = c(pi[1:(ng-1)], unlist(beta), unlist(nu), unlist(delta))
  #     }else{
  #       paraminitEM = c(rep(0,nx), theta, unlist(beta), unlist(nu), unlist(delta))
  #     }
  #   }else{
  #     if (nx == 1){
  #       paraminitEM = paraminit[-ng]
  #     }else{
  #       paraminitEM = paraminit
  #     }
  #   }
  #   param = EMZIPIRLS_cpp(paraminitEM, ng, nx, n, nbeta, nnu, A, Y, X, TCOV, nw, itermax, EMIRLS, refgr)
  #   if (hessian == TRUE){
  #     SE = IEMZIP_cpp(param, ng, nx, nbeta, nnu, n, A, Y, X, TCOV, nw, refgr)
  #     if (nx == 1){
  #       SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)],  sqrt(sum(SE[1:(ng-1)]**2)))
  #     }else{
  #       SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
  #     }
  #   }else{
  #     SE = NA
  #   }
  #   if (nx == 1){
  #     param = c(param[-c(1:ng)], param[1:ng])
  #   }else{
  #     param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
  #   }
  # }
  if (hessian == TRUE){
    if (nx == 1 & Method == "L"){
      paramtmp = c(param[1:(length(param) - ng*nx)], exp(theta)/sum(exp(theta)))
      prob = 1-stats::pt(abs(paramtmp/SE), n*period-1)+stats::pt(-abs(paramtmp/SE), n*period-1) 
    }else{
      prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    }
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = prob
    )
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob = SE)
  }
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  if (nw == 0){
    delta = NULL
    theta = param[-c(1:(sum(nbeta)))]
  }else{
    delta = param[(sum(nbeta)+1):(sum(nbeta)+nw*ng)]
    theta = param[-c(1:(sum(nbeta)+nw*ng))]
  }
  method = ifelse(nx != 1, "L", Method)
  res = list(beta = beta,
             delta = delta,
             theta = theta,
             sd = SE, tab = d, Model = "POIS",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(param = c(theta, beta, delta), model = "POIS", 
                                     method = method, ng =ng, 
                                     nx = nx, n = n, nbeta = nbeta, nw = nw, A = A, Y = Y, X = X, 
                                     TCOV = TCOV),
             Time = A[1,], degre = degre - 1)
  class(res) = "Trajectory.POIS"
  return(res)
}
#################################################################################
# find parameters for Non Linear Model
#################################################################################
#' Internal function to fit Non Linear Model
#' 
#' @inheritParams trajeR
#' @inheritParams trajeR.CNORM
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @return  return a object of class Trajectory.NL
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item sigma - vector of the parameters sigma.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape.
#'   \item fct - the defintion of the function used int this model.
#' }

trajeR.NL <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nw, ntheta, period, degre, theta, beta, sigma, pi, Method, ssigma,
                      hessian, itermax, paraminit, EMIRLS, refgr, fct, diffct, nls.lmiter){
  nsigma = ng
  if (Method == "L"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      sigma = rep(stats::sd(Y),ng)
      paraminit = c(theta, unlist(beta), sigma)
    }
    paraminitL = c(paraminit[1:(ng*nx+sum(nbeta))], log(paraminit[(ng*nx+sum(nbeta)+1):(ng*nx+sum(nbeta)+ng)]))
    if (ssigma == FALSE){
      # different sigma
      newparam = stats::optim(par = paraminitL, fn = LikelihoodalphaNL, gr = difLalphaNL,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=6, maxit = itermax),
                              ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta, 
                              TCOV = TCOV, fct = fct, diffct = diffct)
    } else {
      # same sigma
      newparam = stats::optim(par = paraminitL, fn = LikelihoodalphaNL, gr = difLalphauniqueNL,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=6, maxit = itermax),
                              ng = ng, nx = nx, n =n, A = A, Y = Y, X = X, nbeta = nbeta, 
                              TCOV = TCOV)
    }
    param = newparam$par
    if (nw!=0){
      param = c(param[(ng*nx+1):(ng*nx+sum(nbeta))],exp(param[(ng*nx+sum(nbeta)+1):(ng*nx+sum(nbeta)+ng)]),
                param[(ng*nx+sum(nbeta)+ng+1):(ng*nx+sum(nbeta)+ng+ng*nw)], param[1:(ng*nx)])
    }else{
      param = c(param[(ng*nx+1):(ng*nx+sum(nbeta))],exp(param[(ng*nx+sum(nbeta)+1):(ng*nx+sum(nbeta)+ng)]),
                param[1:(ng*nx)])
    }
    if (hessian == TRUE){
      H = newparam$hessian
      Il = MASS::ginv(-H)
      SE = sqrt(diag(Il))
      SE = c(SE[-c(1:ng)], SE[1:ng])
    }
  } else if (Method == "EM"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      sigma = rep(stats::sd(Y),ng)
      if (nx == 1){
        paraminitEM = c(pi[1:(ng-1)], unlist(beta), sigma)
      }else{
        paraminitEM = c(theta, unlist(beta), sigma)
      }
    }else{
      if (nx == 1){paraminitEM = paraminit[-ng]}
    }
    if (ssigma == FALSE){
      param = EMNL(paraminitEM, ng, nx, nbeta, n, A, Y, X, TCOV, nw, itermax, EMIRLS, fct, diffct, nls.lmiter)
    }else{
      param = EMNLSigmaunique(paraminitEM, ng, nx, nbeta, n, A, Y, X, TCOV,  nw, itermax, EMIRLS, fct, diffct, nls.lmiter)
    }
    if (hessian == TRUE){
      SE = IEMNL(param, ng, nx, nbeta, n, A, Y, X, TCOV, nw, refgr, fct, diffct)
      if (nx == 1){
        SE = c(SE[-c(1:(ng-1))], SE[1:(ng-1)], sqrt(sum(SE[1:(ng-1)]**2)))
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(0, nx), SE[1:((ng-1)*nx)])
      }
    }else{
      SE = NA
    }
    if (nx == 1){
      param = c(param[-c(1:(ng-1))], param[1:(ng-1)], 1-sum(param[1:(ng-1)]))
    }else{
      param = c(param[-c(1:(ng*nx))], param[1:(ng*nx)])
    }
  }
  if (hessian == TRUE){
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    )
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob =SE
    )
  }
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  namesigma = paste0("sigma", 1:ng)
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namesigma, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  if (nw == 0){
    delta = NA
    sigma = param[(sum(nbeta)+1):(sum(nbeta)+ng)]
    theta = param[-c(1:(sum(nbeta)+ng))]
  }else{
    delta = param[(sum(nbeta)+ng+1):(sum(nbeta)+ng+nw*ng)]
    sigma = param[(sum(nbeta)+1):(sum(nbeta)+ng)]
    theta = param[-c(1:(sum(nbeta)+ng+nw*ng))]
  }
  method = ifelse(nx != 1, "L", Method)
  res = list(beta = beta,
             sigma = sigma,
             delta = delta,
             theta = theta,
             sd = SE, tab = d, Model = "CNORM",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(param = c(theta, beta, sigma), model = "NL",
                                     ng = ng, nx = nx, n = n,
                                     nbeta = nbeta, nw = nw, A= A, Y = Y, X = X,
                                     TCOV = TCOV, fct = fct),
             Time = A[1,], degre = degre - 1, fct = fct)
  class(res) = "Trajectory.NL"
  return(res)
}

#################################################################################
# find parameters for Beta Model
#################################################################################
#' Internal function to fit Beta regression
#' 
#' @inheritParams trajeR
#' @inheritParams trajeR.CNORM
#' @param phi Vector of real. The phi parameter.
#' @param nphi Vector of integers. Number of phi parameters for each group.
#' @param X Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix
#' with one column  with value 1.
#' @return  return a object of class Trajectory.NL
#' \itemize{
#'   \item beta -  vector of the parameter beta.
#'   \item sigma - vector of the parameters sigma.
#'   \item delta - vector of the parameter delta. Only if we use time covariate.
#'   \item theta - vector with the parameter theta if there exist a coavriate X that modify
#'   the probability or the probability of group membership.
#'   \item sd - vector of the standard deviation of the parameters.
#'   \item tab - a matrix with all the parameters and standard deviation.
#'   \item Model - a string with the model used.
#'   \item groups -  a integer with the number of group.
#'   \item Names - strings with the name of the parameters.
#'   \item Method  -  a string with the method used.
#'   \item Size - a integer with the number of individuals.
#'   \item Likelihood -  a real with the Likelihood obtained by the parameters.
#'   \item Time - a vector with the first row of time values.
#'   \item degre - a vector with the degree of the polynomial shape.
#' }
trajeR.BETA <- function(Y, A, X, TCOV, ng, nx, n, nbeta, nphi, nw, ntheta, period, degre, theta, beta, phi, delta, pi, Method,
                        hessian, itermax, paraminit, EMIRLS, refgr){
  set_tour(1)
  set_storelik(10**100)
  theta = theta - theta[1:nx]
  theta = theta[-c(1:nx)]
  if (Method == "L"){
    # initial value for Likelihood's method
    if (is.null(paraminit)){
      paraminit = c(theta, unlist(beta), unlist(phi), unlist(delta))
    }
    if (!hessian){
      newparam = stats::optim(par = paraminit, fn = LikelihoodBETA_cpp, gr = difLBETA_cpp,
                              method = "BFGS",
                              hessian = hessian,
                              control = list(fnscale=-1, trace=1, REPORT=1, maxit = itermax),
                              ng = ng, nx = nx, n = n, A = A, Y = Y, X = X, nbeta = nbeta, nphi = nphi,
                              nw = nw, TCOV = TCOV)
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }else{
      newparam = ucminf::ucminf(par = paraminit, 
                                fn = LBETA, 
                                gr = difLBETA,
                                hessian = 2,
                                control = list(stepmax=10**(-5), maxeval = itermax),
                                ng = ng, nx = nx, n = n,  A = A, Y = Y, X = X, nbeta = nbeta, nphi = nphi, 
                                nw = nw, TCOV = TCOV)
      param = newparam$par
      theta = c(rep(0, nx), param[c(1:((ng-1)*nx))])
      param = c(param[-c(1:((ng-1)*nx))], theta)
    }
    if (hessian == TRUE){
      invH = newparam$invhessian
      SE = sqrt(diag(invH))
      if (nx == 1){
        sdtmp = deltaTheta(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        sdbase = deltaThetaBase(theta, invH[1: (ng-1), 1:(ng - 1)], X, ng - 1)
        SE = c(SE[-c(1:((ng-1)*nx))], sdbase, sdtmp)
      }else{
        SE = c(SE[-c(1:((ng-1)*nx))], rep(NA, nx), SE[c(1:((ng-1)*nx))])
      }
    }
  }
  if (hessian == TRUE){
    if (nx == 1 & Method == "L"){
      paramtmp = c(param[1:(length(param) - ng*nx)], exp(theta)/sum(exp(theta)))
      prob = 1-stats::pt(abs(paramtmp/SE), n*period-1)+stats::pt(-abs(paramtmp/SE), n*period-1) 
    }else{
      prob = 1-stats::pt(abs(param/SE), n*period-1)+stats::pt(-abs(param/SE), n*period-1)
    }
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = param/SE,
                   Prob = prob
    ) 
  }else{
    SE = rep(NA,length(param))
    d = data.frame(Estimate =  param,
                   StandardError = SE,
                   TValue = SE,
                   Prob = SE
    )
  }
  
  
  namedegre = c("Intercept", "Linear", "Quadratic", "Cubic", "Quartic", "Quintic", "Sextic", "Septic", "Octic")
  namebeta = c()
  for (i in 1:length(nbeta)){
    namebeta = c(namebeta, namedegre[1:nbeta[i]])
  }
  namephi = c()
  for (i in 1:length(nphi)){
    namephi = c(namephi, namedegre[1:nphi[i]])
  }
  nametheta = rep(c("Intercept", colnames(X)[-1]), ng)
  if (nw == 0){
    namedelta = NULL
  }else{
    namedelta = rep(paste0("TCOV", 1:nw), ng)
  }
  d.names = c(namebeta, namephi, namedelta, nametheta)
  colnames(d) = c("Estimate", "Std. Error", "T for H0 : Parameter=0", "Prob>|T|")
  beta = param[1:(sum(nbeta))]
  phi = param[(sum(nbeta) + 1):(sum(nbeta) + sum(nphi))]
  if (nw == 0){
    delta = NA
    theta = param[-c(1:(sum(nbeta)+sum(nphi)))]
  }else{
    delta = param[(sum(nbeta)+1+sum(nphi)):(sum(nbeta)+nw*ng+sum(nphi))]
    theta = param[-c(1:(sum(nbeta)++sum(nphi)+nw*ng))]
  }
  res = list(beta = beta,
             phi = phi,
             delta = delta,
             theta = theta,
             sd = SE, tab = d, Model = "BETA",
             groups = ng, Names = d.names, Method = Method, Size = n,
             Likelihood = Likelihood(c(theta, beta, phi, delta), model = "BETA", method = Method,
                                     ng = ng, nx = nx, n = n, nbeta = nbeta, nphi = nphi, nw = nw,
                                     A = A, Y = Y, X = X, TCOV = TCOV),
             Time = A[1,], degre = degre - 1, degre.phi = nphi - 1, invH = invH)
  class(res) = "Trajectory.BETA"
  return(res)
}