R/PARETO1.R

In gamlss.dist: Distributions for Generalized Additive Models for Location Scale and Shape

############################################
#  MIKIS STASINOPOULOS #####################
# the PARETO one parameter distribution#####
#######     y>0 AND mu>0       #############
#------------------------------------------------------------
# the fitting function
#-------------------------------------------------------------
PARETO1 <- function (mu.link = "log") 
{
  mstats <- checklink("mu.link", "PARETO1", substitute(mu.link),c("inverse", "log", "sqrt", "identity")) 
  structure(
    list(   family = c("PARETO1", "Pareto one Y>0"),
        parameters = list(mu = TRUE),
             nopar = 1, 
              type = "Continuous", 
           mu.link = as.character(substitute(mu.link)), 
        mu.linkfun = mstats$linkfun, 
        mu.linkinv = mstats$linkinv, 
             mu.dr = mstats$mu.eta, 
              dldm = function(y,mu)
             {
               dldm <- 1/mu-log(1+y) 
               dldm
             },
             d2ldm2 = function(y,mu)
             {
               d2ldm2 <-  -1/mu^2   
               d2ldm2 
             },
     G.dev.incr  = function(y,mu,...) -2*dPARETO1(x = y, mu = mu, log = TRUE),
           rqres = expression(rqres(pfun="pPARETO1", type="Continuous", ymin=0, y=y, mu=mu)), 
      mu.initial = expression({mu <- rep(1,length(y)) } ),
        mu.valid = function(mu) all(mu > 0), 
         y.valid = function(y)  all(y > 0)
    ),
    class = c("gamlss.family","family"))
}
#-----------------------------------------------------------------------------------------
dPARETO1<- function(x, mu = 1, log = FALSE)
{
  if (any(mu <= 0 ) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
#  if (any(x <= 0) )  stop(paste("x must be >0", "\n", ""))
  ly <- max(length(x),length(mu)) 
      x <- rep(x, length = ly)      
     mu <- rep(mu, length = ly)   
   logL <- log(mu) -(mu+1) * log(1+x) #
    lik <- if (log) logL else exp(logL)
    lik <- ifelse(x <= 0, 0, lik)
  as.numeric(lik)
}
#----------------------------------------------------------------------------------------
pPARETO1 <- function(q, mu = 1, lower.tail = TRUE, log.p = FALSE)
{
  if (any(mu <= 0) )  stop(paste("mu must be greater than 0 ", "\n", "")) 
 # if (any(q <= 0) )  stop(paste("q must be >0", "\n", ""))
      cdf <- 1-(1+q)^(-mu)
  if (lower.tail == TRUE) cdf <- cdf  
  else cdf <- 1 - cdf
  if (log.p == FALSE) cdf <- cdf
   else cdf < - log(cdf)
  cdf <-  ifelse(q <= 0, 0, cdf)   
  cdf
}
#----------------------------------------------------------------------------------------
qPARETO1  <- function(p, mu = 1, lower.tail = TRUE, log.p = FALSE)
{      
  if (any(mu <= 0) )  stop(paste("mu 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
     q <- (1/(1-p)^(1/mu))-1
  q
}
#----------------------------------------------------------------------------------------
rPARETO1 <- function(n, mu = 1)
{ 
  if (any(mu <= 0) )  stop(paste("mu 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 <- qPARETO1(p, mu = mu)
  r
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------