R/dmc_irreclass.R

In paulnorthrop/stat2003: stat2003: Introduction to Applied Probability at UCL

#' Finding irreducible classes for discrete-time finite statespace Markov chain
#'
#' \code{dmc_irreclass} focuses on irreducible classes for
#' a given finite discrete-time Markov chain. This function tells you
#' how many irreducible classes in this M.C. and included states
#' of each class and it can also identify whether classes are
#' closed and period of each class.
#'
#' @param MC is an object in class 'stat2003.d'
#' @param period A logical scalar.  If TRUE (the default)
#' \code{dmc_irreclass} will also returns the period of each class.
#' @param whetherc A logical scalar. If TRUE (the default)
#' \code{dmc_irreclass} will also identifies whether classes are closed.
#' @param r_or_t A logical scalar. If TRUE (the default)
#' \code{dmc_irreclass} will also identifies whether classes are
#' (positive) recurrent or transient
#'
#' @return \code{dmc_irreclass} returns a table about all irreducible classes.
#'
#' @seealso \code{\link{stat2003.d-class}}  A class type of discrete-time finite
#' state space Markov chain in stat2003 package.
#' @seealso \code{\link{dmc_equi}} returns the equilibrium
#' distribution for a discrete-time Markov chain.
#' @seealso \code{\link{dmc_inv}} returns the invariant
#' distribution for a discrete-time Markov chain.
#' @seealso \code{\link{dmc_simu}} simulates a discrete-time Markov chain by
#' returning one possible sequence and a states against steps plot.
#' @seealso \code{\link{dmc_tp}} can calculate transient probabilities at a
#' specific step, and also can give a graph about transient probabilities
#' from step zero to that specific step.
#' @seealso \code{\link{dmc_period}} returns the period of each state for
#' a given discrete-time Markov chain
#'
#'
#'
#' @examples
#'
#' m <- matrix(c(0, 1, 0, 0,
#'               1, 0, 0, 0,
#'               1/4, 1/4, 1/4,1/4,
#'               0, 0, 1, 0), nr = 4, nc=4, byrow = TRUE)
#' A <- new("stat2003.d", p_start = c(1/2, 0, 0, 1/2), p = m,
#'          statespace = c("0", "1", "C", "D") )
#'
#' dmc_irreclass(A)
#'
#'
#' m <- matrix(c(1, 0, 0, 0, 0,
#'               0, 1/2, 0, 1/2, 0,
#'               1/2, 0, 0, 1/2, 0,
#'               0, 1, 0, 0, 0,
#'               1/3, 0, 1/3, 0, 1/3), nr = 5, nc=5, byrow = TRUE)
#'
#' A <- new("stat2003.d", p_start = c(1/2, 0, 0, 1/2, 0), p = m,
#'          statespace = c("A!", "B!", "C!", "D!", "E!") )
#'
#' dmc_irreclass(A)
#'
#'
#'
#' m <- matrix(c(0, 1, 0, 0, 0, 0, 0,
#'               0, 0, 1, 0, 0, 0, 0,
#'               0, 1/3, 0, 2/3, 0, 0, 0,
#'               0, 0, 0, 0, 2/3, 1/3, 0,
#'               0, 2/3, 0, 1/3, 0, 0, 0,
#'               0, 0, 0, 0, 0, 1/3, 2/3,
#'               0, 0, 0, 0, 0, 1/4, 3/4), nr = 7, nc=7, byrow = TRUE)
#'
#' A <- new("stat2003.d", p_start = c(0, 0, 0, 1, 0, 0, 0), p = m,
#'          statespace = letters[1:7] )
#'
#' dmc_irreclass(A)
#'
#'
#'@export
dmc_irreclass <- function(MC, period = TRUE, whetherc = TRUE, r_or_t = TRUE) {
  if (class(MC) != "stat2003.d") {
    stop("MC has to be an obejct in class 'stat2003.d' !")
  }

  communi_matrix <- NULL

  # the following loop is used to check whether i communicates with j,
  # if it does, communi_matrix will record 1,
  # if it does not, communi_matrix will record zero.
  for (i in 1 : nrow(MC@p)) {
    for (j in 1 : nrow(MC@p)) {
      p_n <- MC@p
      k <- 1
      # if after 600 * nrow(MC@p) times' moving, p_n[i, j] still equals to zero,
      # it jumps out of the loop.
      while (p_n[i, j] == 0 & k < 600 * nrow(MC@p)) {
        k <- k + 1
        p_n <- p_n %*% MC@p
      }
      if (p_n[i, j] != 0) communi_matrix <- append (communi_matrix, 1)
      else communi_matrix <- append (communi_matrix, 0)
    }
  }

  communi_matrix <- matrix(communi_matrix, nrow = nrow(MC@p),
                           ncol = nrow(MC@p), byrow = T)
  class_vec <- NULL
  class_vec <- matrix(0, nrow = 1, ncol = nrow(MC@p),
                      dimnames = list(NULL, MC@statespace))

  # the following loop is used to check whether i intercommunicates with j,
  # if communi_matrix[i, j] and communi_matrix[j, i] both equal to 1,
  # i and j are intercommunicating.

  # class_vec is used to record class that the state belongs to.

  for (i in 1 : nrow(MC@p)) {
    for (j in 1 : nrow(MC@p)) {
      if (communi_matrix[i, j] == 1 & communi_matrix[j, i] == 1) {
        # if i and j intercommunicate and they have't had a class yet,
        # give them a new class called i
        if (class_vec[1, i] == 0 & class_vec[1, j] == 0) {
          class_vec[1, i] = class_vec[1, j] = i
        }
        # if i and j intercommunicate and one of them has a class,
        # give the other the same class
        if (class_vec[1, i] != 0 & class_vec[1, j] == 0) {
          class_vec[1, j] = class_vec[1, i]
        }
        if (class_vec[1, i] == 0 & class_vec[1, j] != 0) {
          class_vec[1, i] = class_vec[1, j]
          # Note that, it is not possible that i and j intercommunicate
          # while they have different classes,
          # because the loop is in the order of 1, 2, 3 .....
        }
      }
    }
    # if state i does not intercommunicate with others then leave itself in a class.
    if (class_vec[1, i] == 0) class_vec[1, i] = i
  }

  # delete repeated classes and count the total number of classes.
  # number_class <- NULL
  # number_class <- unique(as.vector(class_vec))
  # cat("There are/is", length(number_class), "irresducible class(es)")

  # we start to check whether the classes are closed.
  #
  # in the previous steps,
  # if i communicates with j, we set communi_matrix[i, j] = 1
  # if i does not communicate with j, we set communi_matrix[i, j] = 0.
  # now, we need a intercommuni_matrix that when i and j intercommunicate,
  # intercommuni_matrix[i, j] = intercommuni_matrix[j, i] = 1
  # otherwise intercommuni_matrix[i, j] = intercommuni_matrix[j, i] = 0.
  intercommuni_matrix <- communi_matrix

  for (i in 1 : nrow(MC@p)) {
    for (j in 1 : nrow(MC@p)) {
      if (communi_matrix[i, j] == 1 & communi_matrix[j, i] == 0) {
        # if i and j intercommunicate, we just leave them, because
        # the initial value of intercommuni_matrix = communi_matrix
        # so intercommuni_matrix[i, j] and [j, i] have been set to 1,
        # on the other hand,
        # if i communicates with j but j doesn't communicate with i,
        # we set intercommuni_matrix[i, j] = 0
        intercommuni_matrix[i, j] = 0
      }
    }
  }

  whetherclosed <- NULL
  whetherclosed <- matrix("A closed class", nrow = 1, ncol = nrow(MC@p))

  # the following loop is used to check whether a class is closed.
  #
  # Firstly, we check whether i and j are in the same class
  # if they are, just leave them.
  # if they are not, check whether there is a probaility from i to j
  # if there is, we define that the state itself is not closed,
  # if there isn't, move to the next 'j'.
  for (i in 1 : nrow(MC@p)) {
    for (j in 1 : nrow(MC@p)) {
      if (intercommuni_matrix[i, j] == 0 & MC@p[i, j] != 0) {
        whetherclosed[, i] = "Not a closed class"
      }
    }
  }

  # If there is an unclosed state in a class,
  # this class cannot be a closed class.
  # i.e., if a class is a closed class, all states in this class are closed
  # so let's check it !
  for (i in 1 : nrow(MC@p)) {
    for (j in 1 : nrow(MC@p)) {
      # if i and j are in the same class and one of them is not a closed state
      # we set both of them as "Not a closed class" in order to
      # aggregate them in the future.
      if (class_vec[1, i] == class_vec[1, j] &
          (whetherclosed[1, i] == "Not a closed class" |
           whetherclosed[1, j] == "Not a closed class" )) {
        whetherclosed[1, i] = "Not a closed class"
      }
    }
  }

  # This function will also return the period of each class.
  # As we have already created a 'period' function, use it straight away.
  period2 <- dmc_period(MC)[[1]]


  # The following code is used to format the output
  class_vec <- rbind (class_vec[1, ], whetherclosed, period2, MC@statespace)
  class_vec <- stats::aggregate(class_vec[4, ] ~
                           class_vec[1, ] + class_vec[2, ] + class_vec[3, ],
                         FUN = function(a) {y <- a; return(y)})

  if (class(class_vec[, 4]) == "matrix") {
    colnames(class_vec[, 4]) <- (1 : ncol(class_vec[, 4]))
  }

  class_vec[, 5] <- class_vec[, 2]

  colnames(class_vec) <- c(" ", "whether closed", "period", "states", "recurrent or transient")
  rownames(class_vec) <- paste("Class", 1:nrow(class_vec))

  for (i in 1:nrow(class_vec)) {
    if (class_vec[i, 2] == "Not a closed class") {
      class_vec[i, 5] <- "transient"
    } else {
      class_vec[i, 5] <- "positive recurrent"
    }
  }

  class_vec[, 1] <- NULL

  if (r_or_t == FALSE) {
    class_vec[, 4] <- NULL
  }

  if (period == FALSE) {
    class_vec[, 2] <- NULL
  }

  if (whetherc == FALSE) {
    # have used     class_vec[, 1] <- NULL  before,
    # however this first col is not same as before.
    class_vec[, 1] <- NULL
  }

  return((class_vec))
}