Использование пакета gamlss для получения дистрибутива ex gaussian для использования в пакете depmixs4 - PullRequest
Новая
       52

Использование пакета gamlss для получения дистрибутива ex gaussian для использования в пакете depmixs4

0 голосов
/ 20 февраля 2019

Я моделирую данные о времени реакции, используя скрытое моделирование Маркова с пакетом depmixs4.Настройкой по умолчанию для этого пакета является распределение по Гауссу, однако, поскольку в популяции, которую я использую, часто бывают долгие времена реакции, мы подумали, что было бы лучше использовать распределение по экс-Гауссу.Экс-гауссовский дистрибутив не является одним из параметров семейства по умолчанию в пакете, поэтому я могу использовать вашу помощь в интерпретации кода, чтобы он работал для моих данных: 

setClass("exgaus", contains="response")

#define a generic for the method defining the response class

setGeneric("exgaus", function(y, pstart = NULL, fixed = NULL) standardGeneric("exgaus"))

#define the method that creates the response class

setMethod("exgaus",
          signature(y="ANY"),
          function(y, pstart=NULL, fixed = NULL) {
            y <- matrix(y, length(y))
       x <- matrix(1)
       parameters <- list()
       npar <- 3
       if(is.null(fixed)) fixed <- as.logical(rep(0, npar))
       if(!is.null(pstart)){
         if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
         parameters$mu <- pstart[1]
         parameters$sigma <- log(pstart[2])
         parameters$nu <- log(pstart[3])
       }
       mod <- new("exgaus",parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
       mod
          }
        )

       setMethod("show","exgaus",
                 function(object) {
                   cat("Model of type exgaus (see ?gamlss for details) \n")
                   cat("Parameters: \n")
                   cat("mu: ", object@parameters$mu, "\n")
                   cat("sigma: ", object@parameters$sigma, "\n")
                   cat("nu: ", object@parameters$nu, "\n")
                 }
       )

       setMethod("dens","exgaus",
                 function(object,log=FALSE) {
                   dexGAUS(object@y, mu = predict(object),
                           sigma = exp(object@parameters$sigma), nu = exp(object@parameters$nu), log = log)
                 }
       )

       setMethod("getpars","response",
                 function(object,which="pars"){
                   switch(which,
                          "pars" = {
                            parameters <- numeric()
                            parameters <- unlist(object@parameters)
                            pars <- parameters
                          },
                          "fixed" = {
                            pars <- object@fixed
                          }
                   )
                   return(pars)
                 }
       )

setMethod("setpars","exgaus",
          function(object, values, which="pars", ...) {
            npar <- npar(object)
            if(length(values)!=npar) stop("length of 'values' must be",npar)
            # determine whether parameters or fixed constraints are being set
            nms <- names(object@parameters)
            switch(which,
                   "pars"= {
                     object@parameters$mu <- values[1]
                     object@parameters$sigma <- values[2]
                     object@parameters$nu <- values[3]
                   },
                   "fixed" = {
                     object@fixed <- as.logical(values)
                   }
            )
            names(object@parameters) <- nms
            return(object)
          }
)
setMethod("fit","exgaus",
          function(object,w) {
            if(missing(w)) w <- NULL
            y <- object@y
            fit <- gamlss(y~1,weights=w,family=exGAUS(),
                          control=gamlss.control(n.cyc=100,trace=FALSE),
                          mu.start=object@parameters$mu,
                          sigma.start=exp(object@parameters$sigma),
                          nu.start=exp(object@parameters$nu))
            pars <- c(fit$mu.coefficients,fit$sigma.coefficients,fit$nu.coefficients)
            object <- setpars(object,pars)
            object
          }
)
setMethod("predict","exgaus",
          function(object) {
            ret <- object@parameters$mu
            return(ret)
          }
)

#changed formula=corr~1 to formula=rt~1
#transInit - creates the initial transition probabilities but I don't them?

rModels <- list(
  list(
  exgaus(rt,pstart=NULL),
  GLMresponse(formula=rt~1, data=datanow5,
              family=multinomial("identity"))
),
list(
  exgaus(rt,pstart=NULL),
  GLMresponse(formula=rt~1, data=datanow5,
              family=multinomial("identity"), pstart=NULL
)
),
trstart=NULL,
transition <- list(),
transition[[1]] <- transInit(nstates=2,data=datanow5,pstart=c(trstart[1:2],0,0)),
transition[[2]] <- transInit(nstates=2,data=datanow5,pstart=c(trstart[3:4],0,0)),
instart=c(0.5,0.5),
inMod <- transInit(~1,ns=2,ps=instart,family=multinomial("identity"), data=data.frame(rep(1,3))),
mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,
                  homogeneous=FALSE),
fm3 <- fit(mod,emc=em.control(rand=FALSE)),
summary(fm3))

Выше приведено по умолчаниюКод дан в пакете с некоторыми изменениями на моем конце.Я запустил приведенный выше код и получил следующую ошибку:

Ошибка в (функция (классы, fdef, mtable): не удалось найти унаследованный метод для функции 'transInit' для подписи '' отсутствует ''

У меня нет начальных вероятностей модели (так как она никогда ранее не использовалась), а также я не знаю начальных вероятностей перехода. Могу ли я установить их в NULL или есть некоторые значения по умолчанию, чтобыиспользовать с моделью?

...