library(gamlss)

ZITNO = function (mu.link = "identity", sigma.link = "log", nu.link = "identity")
{
    mstats = checklink("mu.link", "ZITNO", substitute(mu.link),
        c("inverse", "log", "identity", "own"))
    dstats = checklink("sigma.link", "ZITNO", substitute(sigma.link),
        c("inverse", "log", "identity", "own"))
    vstats = checklink("nu.link", "ZITNO", substitute(nu.link),
        c("logit", "probit", "cloglog", "log", "own"))

    structure(list(family = c("ZITNO", "Zero-Inflated Truncated Normal"),
               parameters = list(mu = TRUE, sigma = TRUE, nu = TRUE),
                    nopar = 3,
                     type = "Mixed",
                  mu.link = as.character(substitute(mu.link)),
               sigma.link = as.character(substitute(sigma.link)),
                  nu.link = as.character(substitute(nu.link)),
               mu.linkfun = mstats$linkfun,
            sigma.linkfun = dstats$linkfun,
               nu.linkfun = vstats$linkfun,
               mu.linkinv = mstats$linkinv,
            sigma.linkinv = dstats$linkinv,
               nu.linkinv = vstats$linkinv,
                    mu.dr = mstats$mu.eta,
                 sigma.dr = dstats$mu.eta,
                    nu.dr = vstats$mu.eta,
        


        G.dev.incr = function(y, mu, sigma, nu, ...){  # Global deviance
            -2 * dZITNO(y, mu, sigma, nu, log = TRUE)
        },

        rqres = expression({     # (Normalize quantile) residuals
            uval <- ifelse(y == 0, nu * runif(length(y), 0, 1),
                (1 - nu) * pZINO(y, mu, sigma, nu))
            rqres <- qnorm(uval)
        }),
           mu.initial = expression(mu <- mean(y[y>0])),
        sigma.initial = expression(sigma <- rep(1, length(y))),
           nu.initial = expression(nu <- rep(0.3, length(y))),
             mu.valid = function(mu) TRUE,
          sigma.valid = function(sigma) all(sigma > 0),
             nu.valid = function(nu) all(nu > 0) && all(nu < 1),
              y.valid = function(y) all(y >= 0)),
                class = c("gamlss.family", "family"))

}

#----------------------------------------------------------------------------------------
########## Density  function of ZINO ##########
dZITNO<-function(x, mu=1, sigma=1, nu=.1, log=FALSE)
 {        
          if (any(mu < 0) || any(is.na(mu)))  stop(paste("mu must be positive", "\n", "")) 
          if (any(sigma < 0) || any(is.na(sigma)))  stop(paste("sigma must be positive", "\n", "")) 
          if (any(x < 0))  stop(paste("x must be positive", "\n", ""))  
 log.lik <- ifelse(x==0, log(nu), log(1-nu) + log(dnorm(x,mean=mu,sd=sigma)))
     if(log==FALSE) fy  <- exp(log.lik) else fy <- log.lik
      fy 
  } 

#----------------------------------------------------------------------------------------
########## Acumulate  function distribution of ZINO ##########
pZITNO <- function(q, mu=1, sigma=1, nu=0.1, lower.tail = TRUE, log.p = FALSE)
  {       
    if (any(mu < 0))  stop(paste("mu must be positive", "\n", "")) 
    if (any(sigma < 0))  stop(paste("sigma must be positive", "\n", "")) 
    if (any(q < 0))  stop(paste("y must be positive", "\n", ""))  
    a = sqrt(2/sigma)
    b = (mu*sigma)/(sigma+1)    
    cdf1 <- pnorm((1/a)*(sqrt(q/b) - sqrt(b/q)))
    cdf <- cdf1
   
    if (any(is.na(cdf)))
      {
       index <- seq(along=q)[is.na(cdf)]
       for (i in index)
          {
        cdf[i] <- integrate(function(x) 
                 dbs(x, alpha = a[i], beta = b[i], log=FALSE), 0.001, q[i] )$value
          }    
      }
     cdf <- ifelse((q==0), nu, nu+(1-nu)*cdf)    
    if(lower.tail==TRUE) cdf  <- cdf else  cdf <- 1-cdf 
    if(log.p==FALSE) cdf  <- cdf else  cdf <- log(cdf) 
    cdf
   }   

#----------------------------------------------------------------------------------------
########## Quantile function of ZINO ##########

qZITNO = function (p, mu = 0.5, sigma = 1, nu = 0.1, lower.tail = TRUE,
        log.p = FALSE)
{
    if (any(mu <= 0))
        stop(paste("mu must be positive ", "\n", ""))
    if (any(sigma < 0))  #In this parametrization  sigma = delta
        stop(paste("sigma must be positive", "\n", ""))
    if (any(nu <= 0)|any(nu >= 1))  #In this parametrization  nu = p
        stop(paste("nu must be beetwen 0 and 1 ", "\n", ""))

    if (log.p == TRUE)
        p <- exp(p)
    else p <- p
    if (lower.tail == TRUE)
        p <- p
    else p <- 1 - p
    if (any(p < 0) | any(p > 1))
        stop(paste("p must be between 0 and 1", "\n", ""))

     #a = sqrt(2/sigma)
     #b = (mu*sigma)/(sigma+1)
      
    suppressWarnings(q <- ifelse((nu >= p),0, qnorm((p - nu)/(1-nu),mean = mu, sd = sigma, lower.tail = TRUE, log.p = FALSE)))
    q
} 


#----------------------------------------------------------------------------------------
######## Random generation function of ZINO########

rZITNO = function (n, mu = 0.5, sigma = 1, nu = 0.1)
{
    if (any(mu <= 0))
        stop(paste("mu must be positive", "\n", ""))
    if (any(sigma < 0))  #In this parametrization  sigma = delta
        stop(paste("sigma must be positive", "\n", ""))
    if (any(nu <= 0)|any(nu >= 1))  #In this parametrization  nu = p
        stop(paste("nu must be beetwen 0 and 1 ", "\n", ""))
    if (any(n <= 0))
        stop(paste("n must be a positive integer", "\n", ""))
    n <- ceiling(n)
    p <- runif(n)
    r <- qZITNO(p, mu = mu, sigma = sigma, nu = nu)
    r
}

 # plot the function density of ZINO ##########
plotZITNO = function (mu = .5, sigma = 1, nu = 0.1, from = 0, to = 0.999, n = 101, title="title",
    ...)
{
    y = seq(from = 0.001, to = to, length.out = n)
    pdf = dZITNO(y, mu = mu, sigma = sigma, nu = nu)
    pr0 = c(dZITNO(0, mu = mu, sigma = sigma, nu = nu))
    po = c(0)
    plot(pdf ~ y, main=title, ylim = c(0, max(pdf,pr0)), type = "l",lwd=3)
    points(po, pr0, type = "h",lwd=3)
    points(po, pr0, type = "p", col = "red",lwd=3)
}