R/PIG.R

In gamlss.dist: Distributions to be used for GAMLSS modelling.

## the Poison Inverse Gaussian distribution
## KA RAR and MS
## will not work with CG()
##-----------------------------------------------------------------------------------------
PIG <- function (mu.link = "log", sigma.link = "log") 
{
    mstats <- checklink("mu.link", "Beta Binomial", substitute(mu.link),
                         c("inverse", "log", "identity", "sqrt"))   
    dstats <- checklink("sigma.link", "Beta Binomial", substitute(sigma.link), 
                        c("inverse", "log", "identity", "sqrt"))   
    structure(
          list(family = c("PIG","Poisson.Inverse.Gaussian"),
           parameters = list(mu = TRUE, sigma = TRUE),
                nopar = 2,  
                 type = "Discrete", 
              mu.link = as.character(substitute(mu.link)),
           sigma.link = as.character(substitute(sigma.link)), 
           mu.linkfun = mstats$linkfun, 
        sigma.linkfun = dstats$linkfun, 
           mu.linkinv = mstats$linkinv, 
        sigma.linkinv = dstats$linkinv,
                mu.dr = mstats$mu.eta, 
             sigma.dr = dstats$mu.eta, 
                 dldm = function(y,mu,sigma) {
               ty <- as.double(.C("tofy_", as.single(y), as.single(mu), as.single(sigma),
                                  as.integer(length(y)),as.integer(max(y)+1),  
                                  PACKAGE="gamlss.dist")[[1]])
              dldm <- (y-ty)/mu
              dldm
                                    }, 
               d2ldm2 = function(y,mu,sigma) { 
                                    #d2ldm2 <- eval.parent(quote(-dldp*dldp))
                         ty <- as.double(.C("tofy_", as.single(y), as.single(mu), as.single(sigma),
                                   as.integer(length(y)),as.integer(max(y)+1), PACKAGE="gamlss.dist")[[1]])
                       dldm <- (y-ty)/mu
                     d2ldm2 <- -dldm*dldm
                     d2ldm2 <- ifelse(d2ldm2 < -1e-15, d2ldm2,-1e-15) 
                     d2ldm2
                                    },
                 dldd = function(y,mu,sigma) {
              ty <- as.double(.C("tofy_", as.single(y), as.single(mu), as.single(sigma),
                              as.integer(length(y)),as.integer(max(y)+1),
                              PACKAGE="gamlss.dist")[[1]])
             dldd <- ((ty*(1+sigma*mu)/mu)-(1+sigma*y))/(sigma^2)
              dldd
                                    },
               d2ldd2 = function(y,mu,sigma){
                                    #d2ldd2 <- eval.parent(quote(-dldp*dldp))
                       ty <- as.double(.C("tofy_", as.single(y), as.single(mu), as.single(sigma),
                              as.integer(length(y)),as.integer(max(y)+1),PACKAGE="gamlss.dist")[[1]])
                     dldd <- ((ty*(1+sigma*mu)/mu)-(1+sigma*y))/(sigma^2)
                   d2ldd2 <- -dldd*dldd
                   d2ldd2 <- ifelse(d2ldd2 < -1e-15, d2ldd2,-1e-15)  
                   d2ldd2
                                   },
              d2ldmdd = function(y,mu,sigma) {
             ty <- as.double(.C("tofy_", as.single(y), as.single(mu), as.single(sigma),
                                  as.integer(length(y)),as.integer(max(y)+1),  
                                  PACKAGE="gamlss.dist")[[1]])
                      dldm <- (y-ty)/mu
                      dldd <- ((ty*(1+sigma*mu)/mu)-(1+sigma*y))/(sigma^2)
                   d2ldmdd <- -dldm*dldd
                    d2ldmdd
                                    },
              
           G.dev.incr  = function(y,mu,sigma,pw=1,..) -2*dPIG(y, mu, sigma, log=TRUE),
                rqres = expression(
                  rqres(pfun="pPIG", type="Discrete", ymin=0, y=y, mu=mu, sigma=sigma)
                                   ), 
           mu.initial = expression(mu <- (y+mean(y))/2),
        sigma.initial = expression( sigma <- rep( max( ((var(y)-mean(y))/(mean(y)^2)),0.1),
                                        length(y)) ),
             mu.valid = function(mu) all(mu > 0) , 
          sigma.valid = function(sigma)  all(sigma > 0), 
              y.valid = function(y)  all(y >= 0)
          ),
                class = c("gamlss.family","family"))
}
##-----------------------------------------------------------------------------------------
# 12-3-06 it seems that this function do work with y=0
#t.of.y <- function (y, mu, sigma, bsum = FALSE, ...)
#{
#  tofy <- sumlty <- rep(0,length(y))
#  lofyp1 <- 0
# for (i in seq(along=y))
#    {
#         lofyp1 <- y[i]+1 
#        tofynew <- rep(0,lofyp1)
#       tofynew[1] <- mu[i]*((1+2*sigma[i]*mu[i])^(-0.5))       
#        if (lofyp1==1)
#         {
#          tofy[i] <- tofynew[lofyp1]
#          if (bsum)  sumlty[i] <-  0           
#         }
#        else
#         { 
#         for (j in 2:lofyp1)
#           {
#       tofynew[j] <- ((sigma[i]*(2*(j-1)-1)/mu[i])+(1/tofynew[j-1]))*(tofynew[1])^2
#           }
#          tofy[i] <- tofynew[lofyp1]
#        if (bsum)             sumlty[i] <- sum(log(tofynew))-log(tofynew[j])
#         }
#    }
# result <- cbind(tofy, sumlty)
# result
#}
#-----------------------------------------------------------------------------------------
dPIG<-function(x, mu = 0.5, sigma = 0.02 , log = FALSE)
 { 
          if (any(mu <= 0) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
          if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", "")) 
          if (any(x < 0) )  stop(paste("x must be >=0", "\n", ""))  
          ly <- length(x)                                                                    
      nsigma <- rep(sigma, length = ly)
         nmu <- rep(mu, length = ly)
      sumlty <- as.double(.C("tofy_", as.single(x), as.single(nmu),as.single(nsigma),as.integer(ly),
                          as.integer(max(x)+1),PACKAGE="gamlss.dist")[[2]])
       logfy <- -lgamma(x+1)+(1-sqrt(1+2*sigma*mu))/sigma +sumlty
          if(log==FALSE) fy <- exp(logfy) else fy <- logfy
          fy
  }
##-----------------------------------------------------------------------------------------  
pPIG <- function(q, mu=0.5, sigma=0.02, lower.tail = TRUE, log.p = FALSE)
  {     
  ## function to calculate the cdf
  ## we should write this function in FORTRAN
        tocdf <- function (y, mu, sigma, bsum = TRUE, ...)
        {
         ty <- cdf <- rep(0,length(y))
        for (z1 in length(y):1)
          {
        lyp1 <- dum1 <- dum2 <- 0
        lyp1 <- y[z1]+1 
        tynew <- rep(0,lyp1)
        lpnew <- rep(0,lyp1)
        tynew[1] <- dum1 <- mu[z1]*((1+2*sigma[z1]*mu[z1])^(-0.5)) 
        lpnew[1] <- dum2 <- (1-((1+2*sigma[z1]*mu[z1])^(0.5)))/sigma[z1]
        dum <- ifelse(lyp1==1, 1,2)
         j <- dum:lyp1 
            for (i in j)
            {
            if (i !=1)
            tynew[i] <- ((sigma[z1]*(2*(i-1)-1)/mu[z1])+(1/tynew[i-1]))*(tynew[1]*tynew[1])
            if (i !=1)
                lpnew[i] <- lpnew[i-1] + log(tynew[i-1]) - log(i-1)           
            tynew[1] <- dum1 
            lpnew[1] <- dum2
            ty[z1] <- tynew[lyp1] 
            }
        if (bsum ) 
            cdf[z1] <- sum(exp(lpnew))
          }
        cdf
        }
    ## function tocdf finish here    
          if (any(mu <= 0) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
          if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", "")) 
          if (any(q < 0) )  stop(paste("y must be >=0", "\n", ""))  
           lq <- length(q)                                                                    
       nsigma <- rep(sigma, length = lq)
          nmu <- rep(mu, length = lq)              
          cdf <- tocdf(y = q, mu = nmu, sigma = nsigma, bsum = TRUE)
          if(lower.tail==TRUE) cdf <- cdf else cdf=1-cdf
          if(log.p==FALSE) cdf <- cdf else cdf <- log(cdf)                                                                    
          cdf
   }
##-----------------------------------------------------------------------------------------
qPIG <- function(p, mu = 0.5, sigma = 0.02,  lower.tail = TRUE, log.p = FALSE, 
                 max.value = 10000)
  {      
          if (any(mu <= 0) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
          if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", "")) 
          if (any(p < 0) | any(p > 1.0001))  stop(paste("p must be between 0 and 1", "\n", "")) 
          if (log.p==TRUE) p <- exp(p) else p <- p
          if (lower.tail==TRUE) p <- p else p <- 1-p    
           ly <- length(p)                                                       
          QQQ <- rep(0,ly)                         
       nsigma <- rep(sigma, length = ly)
          nmu <- rep(mu, length = ly)       
        for (i in seq(along=p))                                                          
      { 
       cumpro <- 0                                                                        
     if (p[i]+0.000000001 >= 1) QQQ[i] <- Inf
     else  
        {  
            for (j in seq(from = 0, to = max.value))
            {
           cumpro <-  pPIG(j, mu = nmu[i], sigma = nsigma[i] , log.p = FALSE)  
           QQQ[i] <- j 
       if  (p[i] <= cumpro ) break 
            } 
        }
      }
          QQQ   
   }
##-----------------------------------------------------------------------------------------
rPIG <- function(n, mu = 0.5, sigma = 0.02)
  { 
          if (any(mu <= 0) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
          if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", "")) 
          if (any(n <= 0))  stop(paste("n must be a positive integer", "\n", ""))  
          n <- ceiling(n)
          p <- runif(n)
          r <- qPIG(p, mu=mu, sigma=sigma)
          r
  }
##-----------------------------------------------------------------------------------------