R/upclassifymodel.R

In upclass: Updated Classification Methods using Unlabeled Data

.packageName <- 'upclass'

upclassifymodel <- function (Xtrain, cltrain, Xtest, cltest = NULL, modelName = "EEE", 
                             tol = 10^-5, iterlim = 1000, Aitken = TRUE,  
                             ...) 
{
  functioncall <- match.call(expand.dots = TRUE) 
  
  # get no. of observations in training data
  if (is.matrix(Xtrain)) {
    Ntrain <- nrow(Xtrain)
  }
  else {
    Ntrain <- length(Xtrain)
  }
  # get no of observations in test data
  if (is.matrix(Xtest)) {
    Ntest <- nrow(Xtest)
  }
  else {
    Ntest <- length(Xtest)
  }
  
  grps<-names(table(cltrain)) # need this for the unmap command
  ltrain <- unmap(cltrain)
  G <- ncol(ltrain)
  if (!is.null(cltest)) {
    ltest <-unmap(cltest,groups=grps) 
  }
  if (is.matrix(Xtrain) & (!is.matrix(Xtest))) {
    Xtest <- matrix(Xtest, 1, length(Xtest))
    if (!is.null(cltest)) {
      ltest <- matrix(ltest, 1, length(ltest))
    }
  }
  if (is.matrix(Xtrain)) {
    d <- ncol(Xtrain)
  }
  else {
    d <- 1
  }
  G <- dim(ltrain)[2] 
  fitm <- mstep(modelName = modelName, data = Xtrain, z = ltrain)
  llold <- -Inf
  llstore <- rep(-Inf, 3) # for Aitken
  criterion <- TRUE
  iter <- 0
  while (criterion) {
    iter <- iter + 1
    fite <- do.call("estep", c(list(data = Xtest), fitm))
    emptyz<-TRUE
    if(all(!is.na(fite$z)))
    {
      emptyz<-FALSE
      z <- fite$z
      zall <- rbind(ltrain, z)
      if (d > 1) {
        Xall <- rbind(Xtrain, Xtest)
      }
      else {
        Xall <- c(Xtrain, Xtest)
      }
      fitm <- mstep(modelName = modelName, data = Xall, z = zall)
      mTau.train <- matrix(log(fitm$parameters$pro), Ntrain, 
                           fitm$G, byrow = TRUE)
      lDensity.train <- do.call("cdens", c(list(data = Xtrain, 
                                                logarithm = TRUE), fitm))
      sum.train <- mTau.train + lDensity.train
      mat.train <- ltrain * sum.train
      ll.train <- sum(mat.train)
      ll.test <- sum(do.call("dens", c(list(data = Xtest, logarithm = TRUE), 
                                       fitm)))
      ll <- ll.train + ll.test
      llstore <- c(llstore[-1], ll)
      if (Aitken) {
        criterion <- (Aitken(llstore)$linf - ll) > tol
      }
      else {
        criterion <- (ll - llold) > tol
      }
      criterion <- (criterion) & (iter < iterlim)
      llold <- ll
    }
    else {
      criterion<-FALSE
    }
  }
  res <- list() # is this an empty list
  if(!emptyz)
  {
    res$call <- functioncall
    res$Ntrain <- Ntrain
    res$Ntest <- Ntest
    res$d <- d
    res$G <- G
    res$iter <- iter
    res$converged <- (iter < iterlim)
    res$modelName <- modelName
    res$parameters <- fitm$parameters
    
      fitetrain <- do.call("estep", c(list(data = Xtrain), 
                                      fitm))
      ztrain <- fitetrain$z
      cltrain <- map(ztrain)
      tab0 <- table(map(ltrain), cltrain)
      tab <- matrix(0, G, G)
      rownames(tab) <- colnames(tab) <- 1:G
      tab[rownames(tab0), colnames(tab0)] <- tab0
      misclasstrain <- sum(tab) - sum(diag(tab))
      Brier <- sum((ltrain - ztrain)^2)*100/(2*Ntrain)
      res$train <- list()
      res$train$z <- ztrain
      res$train$cl <- cltrain
      res$train$misclass <- misclasstrain
      res$train$rate <- misclasstrain*100/Ntrain
      res$train$Brierscore <- Brier
      res$train$tab <- tab
      cl <- map(z)
      res$test <- list()
      res$test$z <- z
      res$test$cl <- cl
      if (!is.null(cltest)) {
        cltest <- map(ltest)
        tab0 <- table(cltest, cl)
        tab <- matrix(0, G, G)
        rownames(tab) <- colnames(tab) <- 1:G
        tab[rownames(tab0), colnames(tab0)] <- tab0
        ratetest <- sum(tab) - sum(diag(tab))
        Brier <- sum((ltest - z)^2)*100/(2*Ntest)
        res$test$misclass <- ratetest
        res$test$rate <-ratetest*100/length(cltest)
        res$test$Brierscore <- Brier
        res$test$tab <- tab
      }
  #  }
    bic.all <- bic(modelName, ll, fitm$n, fitm$d, fitm$G)
    res$ll <- ll
    res$bic <- bic.all
  }
  else {
    res$error<-"Training groups too small for this model type"
    res$ll<-NA
    res$bic <- NA
  }
  res
}