R/gamlssMX-10-11-07.R

#-------------------------------------------------------------------------------
## this  a general function for fitting mixures within gamlss
## created by Mikis Stasinopoulos and Bob Rigby 
## LAST CHECKED 5-8-14 
## latest change Tuesday, April 10, 2007 at 08:54
## Some of the models which can be fitted 
## with the function gamlssNP can be identical 
## In this version the prior probabilities pi can be modelled 
## TO DO 
## it needs a summary functions  but the EM-algorithm SE's are not correct
## do we need EBP here?
## do we need means? 
## the plotMX() need modification for BI type data Done but not checked
## 
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
gamlssMX  <- function ( formula = formula(data), 
                     pi.formula = ~1, 
                         family = "NO", # note is character
                         weights, 
                              K = 2, 
                           prob = NULL,
                           data = sys.parent(),
                        control = MX.control(...), 
                      g.control = gamlss.control(trace=FALSE,...),
                 zero.component = FALSE,     
                      ...)
 {
#-------------------------------------------------------------------------------
# this is to replicate rqres within gamlss enviroment DS Friday, March 31, 2006 at 10:30
rqres <- function (pfun = "pNO", 
                   type = c("Continuous", "Discrete", "Mixed"),
                   censored = NULL,  
                   ymin = NULL, 
                   mass.p = NULL, 
                   prob.mp = NULL,
                   y = y,
                   ... )
{ }
body(rqres) <-  eval(quote(body(rqres)), envir = getNamespace("gamlss"))
##------------------------------------------------------------------------------
##------------------------------------------------------------------------------
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
# extra functions within -------------------------------------------------------
# for getting the commulative function------------------------------------------
get.the.p.function <- function(object, ...)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(object))  stop(paste("This is not an gamlss object", "\n", ""))
   fname <- object$family[[1]]
 DistPar <- object$parameters
   nopar <- length(DistPar) 
    dfun <- paste("p",fname,sep="")
 # binomial denominators
 switch(nopar,  
    {pfun <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,q=object$y, bd=object$bd,  mu=fitted(object,"mu")))
             else  eval(call(dfun,q=object$y, mu=fitted(object,"mu")))},
    {pfun <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,q=object$y, bd=object$bd,  mu=fitted(object,"mu"),sigma=fitted(object,"sigma")))
             else  eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"))) },
    {pfun <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,q=object$y, bd=object$bd,  mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))
             else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))},
    {pfun <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,q=object$y, bd=object$bd,  mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))
             else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))})
pfun 
}
#-------------------------------------------------------------------------------
get.likelihood <- function(obj)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(obj))  stop(paste("This is not an gamlss object", "\n", ""))
   fname <- obj$family[[1]]
 DistPar <- obj$parameters
   nopar <- length(DistPar) 
    dfun <- paste("d",fname,sep="")
 switch(nopar,  
    {lik <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,x=obj$y, bd=obj$bd,  mu=fitted(obj,"mu")))
             else   eval(call(dfun, x=obj$y, mu=fitted(obj)))},
    {lik <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,x=obj$y, bd=obj$bd,  mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma") ))
             else   eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"))) },
    {lik <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,x=obj$y, bd=obj$bd,  mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu") ))
            else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu")))},
    {lik <- if (fname%in%.gamlss.bi.list)  eval(call(dfun,x=obj$y, bd=obj$bd,  mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau") ))
            else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau")))})
lik
}           
#-------------------------------------------------------------------------------
expand.vc <- function(x,ni)
{
  if (length(ni)==1)
  {
      if (ni<1) stop("ni should be greater than 1")
      xx <- x
      if (ni==1) return(xx)
      for ( i in 2:ni) xx <- rbind(x,xx)
      xx
   } 
   else 
   {
      n <- dim(x)[[1]]
      c <- dim(x)[[2]]
     xx <- NULL
      for ( i in seq(1,n))
      {
          xx <- rbind(xx,matrix(rep(x[i,],ni[i]),ncol=c,byrow=TRUE))
      }
      xx
  }
}
#-------------------------------------------------------------------------------
# the proper function starts here
#-------------------------------------------------------------------------------
#library(gamlss)
#library(nnet)
 gamlssMXcall <- match.call()  #   the function call
## checking for NA in the data 
 if(!missing(data) & any(is.na(data)))   
      stop("The data contains NA's, use data = na.omit(mydata)") 
## prob are the starting values of the pi's   
 if(zero.component == TRUE)
 { KK <- K+1 
  if (is.null(prob)) 
   { 
     prob <- rep(1/KK, KK)
   }
  else
   {
    if (length(prob)!=KK) stop("the length of prob must be  K+1(for zero component) ")
   } 
 }
 else
 {
   if (is.null(prob)) 
     { 
     prob <- rep(1/K, K)
     }
 else
   {
    if (length(prob)!=K) stop("the length of prob must be K")
   } 
 }
## expand the formula if not list
 if (!is.list(formula))
  {
   allFormula <- vector("list",K)
   for (i in 1:K) allFormula[[i]] <- formula
  }
 else allFormula <- formula
## expand the family if not list
 if (!is.list(family))
  {
   allFamily <- vector("list",K)
   for (i in 1:K) allFamily[[i]] <- family
  }
  else allFamily <- family
## check that the length for both the family and formula is K
if (length(allFamily)!=K) stop("the length of the list for the family is not equal the number of componets in the mixture")
if (length(allFormula)!=K) stop("the length of the list for the formula is not equal the number of componets in the mixture")
## get response for his length
        Y <- model.extract(model.frame(allFormula[[1]],data=data),"response")
        N <-  if(is.null(dim(Y))) length(Y) else dim(Y)[1]# calculate the dimension for y  
pweights. <- pweights. <<- if (missing(weights)) rep(1,N) else weights  
## get things from control
  set.seed(control$seed) # set the seed 
prob.sample <- if (is.null(control$sample)) 10/N else control$sample  
## creating the matrix of posterior probabilities (weights)
        W <- matrix(1, nrow=N, ncol=K ) 
allModels <- vector("list",K)
## get starting values for the weights (posterior probabilities)
 for (i in 1:K)
   {         wSam <- sample(c(0,1), N, replace = TRUE, prob=c(1-min(.5,prob.sample),min(.5,prob.sample)))  
              ww. <- ww. <<- wSam # W[,i] # put in 
             #data1 <<- data.frame(data, ww, pweights)
    allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control, ...)
             W[,i] <- get.likelihood(allModels[[i]]) 
   } 
 if(zero.component == TRUE) W[,KK] <- prob[KK]*ifelse(Y==0,1,0)
      SumLik <- rowSums(W)
           W <- W/SumLik
           W <- ifelse(W<1e-10,0,W)
        prob <- colSums(W)/N
 # -2 * logLikelihood
       newdv <- -2*sum(log(SumLik))    # the global deviance
       olddv <- newdv + 1  
    iter.num <- 1 
 trace.print <- 0 
    dev.fits <- rep(0,control$n.cyc)
#-------------------------------------------------------------------------------
# if require model for prior propabilities pi's 
     modelPi <- if (pi.formula[[2]]==1) FALSE else TRUE
  if (modelPi)
   {
    if (missing(data)) stop("the data argument is needed if pi's are modelled")
     if(zero.component == TRUE) 
      {
      dataKK <- expand.vc(data,KK)# expand data.
     fac.fit <- gl(KK,N)
      }
      else
      { 
      dataKK. <<- expand.vc(data,K)# expand data
      fac.fit <- gl(K,N) # get the factor
      }
     res.var <- model.matrix(~fac.fit-1) # get the matrix of zeroes ans ones
   form.prob <-  update(pi.formula, res.var~.)# formula(paste("res.var~", pi.formula[[2]])) # get the formula 
        PROB <-  matrix(prob, ncol=K, nrow=N, byrow=TRUE) # expand prior prob
   }
## EM starts here--------------------------------------------------------------- 
 while ( abs(olddv-newdv) > control$cc && iter.num < control$n.cyc ) # MS Wednesday, June 26, 2002 
  {
  for (i in 1:K)
    {          ww. <- ww. <<- W[,i]
    allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control,...)
             W[,i] <- if (modelPi) PROB[,i]*get.likelihood(allModels[[i]]) else prob[i]*get.likelihood(allModels[[i]])
    } 
     if(zero.component == TRUE)
      {
      W[,KK] <-  if (modelPi) PROB[,KK]*ifelse(Y==0,1,0) else prob[KK]*ifelse(Y==0,1,0)  
      }
            SumLik <-rowSums(W)
                 W <- W/SumLik 
                 W <- ifelse(W<1e-10,0,W) 
                di <- -2*log(SumLik) 
             olddv <- newdv
             newdv <- sum(di)
       if (modelPi)
          {
          dataKK.$wWw <- wWw <- as.vector(W)
          dataKK.$res.var <- res.var
          mult.mod <- multinom(form.prob, weights=wWw, data=dataKK., trace=FALSE)
           fitted(mult.mod)[1:11,]
              PROB <- fitted(mult.mod)[1:N,] # take only N of them 
          }
       else { prob <- colSums(W)/N}
   if (control$trace == TRUE) cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n"  )
   else if (is.numeric(control$trace))
        { 
          trace.print <- trace.print+1
         if (trace.print==control$trace)
          {cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n"  )
          trace.print <- 0 
          }
        }
dev.fits[iter.num] <- newdv 
          iter.num <- iter.num+1    
  }
if (control$plot==TRUE) plot(dev.fits[dev.fits!=0], type="l", xlab="EM iterations", ylab="-2 logLik")
## EM finish here--------------------------------------------------------------- 
## the output starts here 
## residuals -------------------------------------------------------------------
                WF <- WF2 <- matrix(0, ncol=K, nrow=N)
   if (modelPi)
    {
     for (i in 1:K)     WF[,i] <- PROB[,i]*get.the.p.function(allModels[[i]]) 
    }
    else 
    {
    for (i in 1:K)     WF[,i] <- prob[i]*get.the.p.function(allModels[[i]]) 
    }
               res <-qnorm(rowSums(WF))
    for (i in 1:K) WF2[,i] <- W[,i]*get.the.p.function(allModels[[i]])
              res2 <-qnorm(rowSums(WF2))   
#------------------------------------------------------------------------------- 
## degrees of freedom
df.fit <- 0 
for (i in 1:K) df.fit <- df.fit + allModels[[i]]$df.fit
df.fit <- if(modelPi) df.fit+mult.mod$edf  else df.fit+K-1 # 
if (modelPi)  colnames(PROB) <- paste("pi",seq(1,K), sep="")
familyAll <- rep("NULL",K)
for (i in 1:K) familyAll[i] <- allModels[[i]]$family[1]
## all output
out <- list()
out <- list(models = allModels,
          model.pi = if (modelPi) mult.mod else NULL,
        G.deviance = newdv, 
            df.fit = df.fit, 
       df.residual = N-df.fit,
         post.prob = W, 
              prob = if (modelPi) PROB else prob,
            family = familyAll,
              call = gamlssMXcall,
               aic = newdv+2*df.fit,
               sbc = newdv+log(N)*df.fit,
                 K = K,
                 N = N,
           weights = pweights., 
         residuals = res,
        Presiduals = res2,
              seed = control$seed
                 )
class(out) <- list("gamlssMX", "gamlss")
if (modelPi) on.exit(rm(ww.,pweights.,dataKK., envir=sys.frame(0)))
else on.exit(rm(ww.,pweights., envir=sys.frame(0)))
out
 }  
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
MX.control <- function(cc=0.0001, n.cyc=200, trace=FALSE, seed=NULL, plot=TRUE, sample=NULL, ...)
{
  ##list(cc=0.0001, n.cyc=100, trace=TRUE),
  ##  Control iteration for GAMLSS
  ##  Mikis Stasinopoulos Monday, March 25, 2002 at 16:17
  ##
  if(cc <= 0) {
    warning("the value of cc supplied is zero or negative the default value of 0.0001 was used instead")
    c.crit <- 0.0001}
  if(n.cyc < 1) {
    warning("the value of no cycles supplied is zero or negative the default value of 20 was used instead")
    n.cyc <- 100}
  if(is.logical(trace)) trace <- trace 
  else if (is.numeric(trace) & trace <= 0) 
  {warning("the value of trace supplied is less or equal t zero the default of 1 was used instead")
   trace <- 1
  }
  seed <- if (is.null(seed) ) floor(runif(min=0, max=100000, n=1)) else seed 
  if (!is.null(sample)) 
  { if (sample>1||sample<0)
  {
    sample <-0.1 
    warning("the value of sample supplied is not probability the default of 0.10 was used instead")
  }
  else sample <- sample  
  } 
  plot <- if(is.logical(plot)) plot else TRUE     
  list(cc = cc, n.cyc = n.cyc, trace = trace, seed=seed, plot=plot, sample=sample)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# METHODS for "gamlssMX"
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
print.gamlssMX <- function (x, digits = max(3, getOption("digits") - 3), ...) 
{
    K <- x$K 
    cat("\nMixing Family: ", deparse(x$family), "\n", fill=TRUE)
    cat("Fitting method: EM algorithm \n")
    cat("\nCall: ", deparse(x$call,width.cutoff=50),  "\n", fill=TRUE)
     for (i in 1:K)
       {
        if ("mu" %in% x$models[[i]]$parameters) 
         {      
          cat("Mu Coefficients for model:", i,  "\n")   
          print.default(format(coef(x$models[[i]], "mu"), digits = digits), print.gap = 2, quote = FALSE)        
         }
        if ("sigma" %in% x$models[[i]]$parameters) 
         { 
          cat("Sigma Coefficients for model:", i,  "\n")           
          print.default(format(coef(x$models[[i]], "sigma"), digits = digits), print.gap = 2, quote = FALSE)        
         }
        if ("nu" %in% x$models[[i]]$parameters) 
         {       
         cat("Nu Coefficients for model:", i,  "\n")    
         print.default(format(coef(x$models[[i]], "nu"), digits = digits), print.gap = 2, quote = FALSE)        
         }
        if ("tau" %in% x$models[[i]]$parameters) 
         {
          cat("Tau Coefficients for model:", i,  "\n")          
          print.default(format(coef(x$models[[i]], "tau"), digits = digits), print.gap = 2, quote = FALSE)        
         }
       }
       if(is.matrix(x$prob)) 
         {
          cat("model for pi: \n")
          print(coef(x$model.pi))
          cat("\nEstimated probabilities: \n")
          print(x$prob[1:3,])
          cat("... \n" )
         }
       else   cat("\nEstimated probabilities:", x$prob, "\n" )
    cat("\nDegrees of Freedom for the fit:", x$df.fit, "Residual Deg. of Freedom  ", 
        x$df.residual, "\n")
   cat("Global Deviance:    ", format(signif(x$G.deviance)), 
        "\n            AIC:    ", format(signif(x$aic)), "\n            SBC:    ", 
        format(signif(x$sbc)), "\n")
    invisible(x)
}

################################################################################
#                         fitted.gamlssMX
################################################################################
# this gives the compoment fv if K is set otherwise average the componets 
# using the fitted probabilities
# MS is this justified if mu not the mean? 
#  also does it makes sense for other parameters??
fitted.gamlssMX<-function (object, K=1, ... ) 
{
if (K%in%seq(1:object$K))
 {
  x <- fitted(object$models[[K]], ...)
 }
else 
  {
  WF <- matrix(0, ncol=object$K, nrow=object$N)
  for (i in 1:object$K)     WF[,i] <- object$prob[i]*fitted(object$models[[i]]) 
  x <-rowSums(WF)
  }
 x
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# coefficients of the componets
coef.gamlssMX <- function(object, K=1, ...)
{
coef(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
formula.gamlssMX <- function(x, K=1, ...)
{
formula(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
model.matrix.gamlssMX <- function(object, K=1, ...)
{
model.matrix(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
terms.gamlssMX <- function(x, K=1, ...)
{
terms(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
predict.gamlssMX <- function(object, K=1, ...)
{
predict(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
residuals.gamlssMX <- function(object,type=c("prior","post"), ...)
{
  type <- match.arg(type)
  switch(type,"prior"=object$residuals,"post"= object$Presiduals)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values 
# and select the one with smallest deviance 
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values 
# and select the one with smallest deviance 
# last modification Thursday, April 30, 2009 
gamlssMXfits <- function(    n = 5,
                       formula = formula(data), 
                    pi.formula = ~1, 
                        family = "NO", # note is character
                        weights, 
                             K = 2, 
                          prob = NULL,
                          data = sys.parent(),
                       control = MX.control(), 
                     g.control = gamlss.control(trace=FALSE),
                     zero.component = FALSE,   
                     ...)
{
 gamlssMXfitscall <- match.call()  
        modellist <- list()
              dev <- rep(0,length=n) 
             seed <- floor(runif(min=0,max=10000, n=n))
for (i in (1:n))
 {
  m1 <- try(gamlssMX(formula=formula,  pi.formula = pi.formula, family = family ,K=K,  prob=prob, data = 
                    data, control=MX.control(seed=seed[i]), g.control = gamlss.control(trace=FALSE), 
                    zero.component=zero.component)) #
             if (any(class(m1)%in%"try-error")) 
               { 
                cat("model=", i, "failed", "\n") 
                dev[i] <- NA    
        modellist[[i]] <- NA
                 next
               }
 modellist[[i]] <- m1 
         dev[i] <- deviance(modellist[[i]])
 cat("model=", i,"\n")                    
 }
                   II <- which.min(dev)
                model <- modellist[[II]]
   gamlssMXfitscall$n <- NULL
gamlssMXfitscall[[1]] <- as.name("gamlssMX")
           model$call <- gamlssMXfitscall
          model$extra <- list(dev=dev, seed=seed, which=II)
 model
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
update.gamlssMX <- function (object, 
                          formula., 
                          ..., 
                          what = c("mu", "sigma", "nu", "tau"),
                          parameter= NULL,
                          evaluate = TRUE) 
{
    call <- object$call
    if (is.null(call)) 
        stop("need an object with call component")
    extras <- match.call(expand.dots = FALSE)$...
    if (!missing(formula.)) 
        {
      what <- if (!is.null(parameter))  {
        match.arg(parameter, choices=c("mu", "sigma", "nu", "tau"))} else  match.arg(what)
      if (!what%in%object$par) stop(paste(what,"is not a parameter in the object","\n")) 
        if (what=="mu") 
         { call$formula <- update.formula(formula(object,what), formula.) }
        else  
         {
          call[[paste(what,"formula",sep=".")]] <- 
         if (length(update.formula(formula(object,what), formula.))==2)
          update.formula(formula(object,what), formula.)
         else
          update.formula(formula(object,what), formula.)[-2]
         }
        }
    if (length(extras) > 0) 
        {
        existing <- !is.na(match(names(extras), names(call)))
        for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
        if (any(!existing)) 
           {
            call <- c(as.list(call), extras[!existing])
            call <- as.call(call)
           }
        }
    if (evaluate) 
        eval(call, parent.frame())
    else call
}
################################################################################
#
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the pdf of y fY(y)=p1*f1(y)+p2*f(y), of a MX model
# MS + BR 
dMX <- function(y, 
                 mu = list(mu1=1,mu2=5), 
              sigma = list(sigma1=1,sigma2=1), 
                 nu = list(nu1=1,nu2=1), 
                tau = list(tau1=1,tau2=1), 
                 pi = list(pi1=.2,pi2=.8),
                # K = 2, 
             family = list(fam1="NO", fam2="NO"),  
                log = FALSE,
                 ...)
  {
#  gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
#  .gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
 ## checking the probabilities 
        K <- length(pi) 
     sump <- sum(unlist(pi))
       ## set the length using the length of propabilities
#  for (i in 1:K) sump <-pi[[i]]+ sump 
  if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
 ## get the length of y and the length of the parameters
 if(is.null(dim(y))) N <- length(y) else N <- dim(y)[1]  
       Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K)) 
     for (i in 1:K) # go throught the components
   { 
        fam <-  as.gamlss.family(family[[i]])
      fname <- fam$family[[1]]
    DistPar <- fam$parameters
      nopar <- length(DistPar)  
       dfun <- paste("d",fname,sep="") 
       if (fname%in%.gamlss.bi.list) # for binomial data
    {
     if (NCOL(y) == 1)
        {
        y <- y 
       bd <- rep(1, length(y))
        }
     else if (NCOL(y) == 2)
        {
        y <- y[,1]
       bd <- y[,1]+y[,2]
        }
     else
        {
      stop("wrong response variable for binomial data")
        }
     }           
   switch(nopar,  
      {lik <- if (fname%in%.gamlss.bi.list){eval(call(dfun, x=y, bd=bd,  mu=mu[[i]]))}
            else  eval(call(dfun, x=y, mu=mu[[i]]))
      },
      {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd,  mu=mu[[i]], sigma=sigma[[i]]))
            else  eval(call(dfun, x=y, mu=mu[[i]], sigma=sigma[[i]])) 
      },
      {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd,  mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]])) 
            else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
      {lik <-  if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd,  mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]])) 
            else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
    Prob[,i]<-pi[[i]]*lik
    }  # finish go throught the components
       fy <- rowSums(Prob)
       fy <- if(log == FALSE) fy else log(fy)
       fy 
  }
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the cdf y from a MX model
pMX <- function(q, 
                 mu = list(mu1=1,mu2=5), 
              sigma = list(sigma1=1,sigma2=1), 
                 nu = list(nu1=1,nu2=1), 
                tau = list(tau1=1,tau2=1), 
                  pi = list(pi1=.2,pi2=.8),
                #  K = 2, 
             family = list(fam1="NO", fam2="NO"),  
                 log = FALSE,
                 ...)
  {
 # gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
## checking the probabilities 
          K <- length(pi) 
       sump <- sum(unlist(pi))
   if (sump!=1) stop(paste("the vector pi should sum to 1"))          
#  for (i in 1:K) sump <-pi[[i]]+ sump 
  if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
 ## get the length of q and the length of the parameters
 if(is.null(dim(q))) N <- length(q) else N <- dim(q)[1]  
       Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K))  
     for (i in 1:K)
   { 
        fam <- as.gamlss.family(family[[i]])
      fname <- fam$family[[1]]
    DistPar <- fam$parameters
      nopar <- length(DistPar)  
       dfun <- paste("p",fname,sep="")
        if (fname%in%.gamlss.bi.list) # for binomial data
    {
     if (NCOL(q) == 1)
        {
        q <- q 
       bd <- rep(1, length(q))
        }
     else if (NCOL(q) == 2)
        {
        q <- q[,1]
       bd <- q[,1]+q[,2]
        }
     else
        {
      stop("wrong response variable for binomial data")
        }
     }            
   switch(nopar,  
     {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]],  mu=mu[[i]]))
            else   eval(call(dfun, q=q, mu=mu[[i]]))
     },
     {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]],  mu=mu[[i]], sigma=sigma[[i]]))
            else   eval(call(dfun, q=q, mu=mu[[i]], sigma=sigma[[i]])) 
     },
     {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]],  mu=mu[[i]], sigma=sigma[[i]],  nu=nu[[i]]))
            else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
     {lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]],  mu=mu[[i]], sigma=sigma[[i]],  nu=nu[[i]], tau=tau[[i]]))
            else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
    Prob[,i]<-pi[[i]]*lik
    } 
            fy <-rowSums(Prob)
       fy <- if(log == FALSE) fy else log(fy)
       fy 
  }
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# A function creating a function which can be use to plot a 
# fitted mixture distrubution
# it needs a gamlssMX object and the observation number 
getpdfMX <- function(object=NULL, observation=1)
{
  if (class(object)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
          K <- object$K
     family <- object$family
  ParamList <- list()
         MU <- SIGMA <- NU<- TAU <- PI <- npar<- list()
modelForPi  <- class(object$model.pi)[1]=="multinom"
  for (i in 1:K)
  {
    ParamList[[i]] <- object$models[[i]]$parameters
    if ("mu"%in%ParamList[[i]])
          MU[[i]] <- fitted(object, K=i)[observation]
    if ("sigma"%in%ParamList[[i]])
       SIGMA[[i]] <- fitted(object, K=i, parameter="sigma")[observation]
    if ("nu"%in%ParamList[[i]])
          NU[[i]] <- fitted(object, K=i, parameter="nu")[observation]
    if ("tau"%in%ParamList[[i]])
         TAU[[i]] <- fitted(object, K=i, parameter="tau")[observation]
          PI[[i]] <- if (modelForPi)  object$prob[observation,i]  else object$prob[i]
        npar[[i]] <- length(ParamList[[i]])
  }  
        npar <- unlist(npar)
  if (any(npar==npar[1])) 
    nopar <- npar[1] else 
      stop("the different distributions should have same number of parameters for the function to work")
  plotfun <- function(y)
  {  
    if (object$family[1]%in%.gamlss.bi.list) 
    {
      yy <- object$y
      bd <- object$bd 
      switch(nopar,
             dMX(y, bd=bd, mu=MU, pi=PI, family=family),      
             dMX(y, bd=bd, mu=MU, sigma=SIGMA, pi=PI, family=family),
             dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
             dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))  
    } else 
    {
      switch(nopar,
             dMX(y, mu=MU, pi=PI, family=family),      
             dMX(y, mu=MU, sigma=SIGMA, pi=PI, family=family),
             dMX(y, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
             dMX(y, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))
    }  
  }
  plotfun
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# at the moment only for -Inf +Inf continuous 
# " discrete? 
meanMX <- function(obj, observations=NULL, limit=c(-Inf, Inf))
{
  # checking whether a gamlss object
  if (class(obj)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
  #   
   obs <- if (is.null(observations)) 1:obj$N else observations
  mean <- rep(0, length(obs))
     j <- 1
  for (i in obs)
  {
    fnLO <- getpdfMX(obj, observation=i)
    f <-function(x) x*fnLO(x)
    mean[j] <- integrate(f, limit[1], limit[2])$value
    j <- j+1
  }
  mean
}
# 
# 

Try the gamlss.mx package in your browser

Any scripts or data that you put into this service are public.

gamlss.mx documentation built on May 29, 2024, 10:33 a.m.