R/dmc_equi.R

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

#' Equilibrium distribution of a discrete-time finite statespace Markov chain
#'
#' \code{dmc_equi} is used to find the equilibrium distribution of
#' a discrete-time finite statespace Markov chain.
#'
#' @param MC is an object in class 'stat2003.d'
#'
#' @return \code{dmc_equi} will return a vector if the Marchov chain
#' has equilibrium distribution.
#'
#' @seealso \code{\link{stat2003.d-class}}  A class type of discrete-time finite
#' state space Markov chain in stat2003 package.
#' @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_irreclass}} focuses on irreducible classes for
#' a given discrete-time Markov chain
#' @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), nrow = 4, ncol=4, byrow = TRUE)
#'
#' A <- new("stat2003.d", p_start = c(1/2, 0, 0, 1/2), p = m,
#'           statespace = letters[1 : 4] )
#'
#' dmc_equi(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), nrow = 5, ncol=5, byrow = TRUE)
#'
#' A <- new("stat2003.d", p_start = c(0, 0, 0, 1, 0), p = m,
#'          statespace = c("1", "2", "C", "4", "E") )
#'
#' dmc_equi(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), nrow = 7, ncol=7, byrow = TRUE)
#'
#' A <- new("stat2003.d", p_start = c(0, 0, 0, 1, 0, 0, 0), p = m,
#'         statespace = LETTERS[1 : 7] )
#'
#' dmc_equi(A)
#'
#'
#' m <- matrix(c(1, 0, 0, 0, 0, 0,
#'               1/4, 0, 3/4, 0, 0, 0,
#'               0, 1/2, 0, 1/2, 0, 0,
#'               0, 0, 1/2, 0, 1/3, 1/6,
#'               0, 0, 0, 0, 1/4, 3/4,
#'               0, 0, 0, 0, 1/3, 2/3), nrow = 6, ncol=6, byrow = TRUE)
#' A <- new("stat2003.d", p_start = c(0, 0, 0, 1, 0, 0), p = m,
#'         statespace = c("1", "2", "3", "4", "5", "6") )
#'
#' dmc_equi(A)
#'
#'
#'@export
dmc_equi <- function (MC) {
  if (class(MC) != "stat2003.d") {
    stop("MC has to be an object in class 'stat2003.d' !")
  }

  # An equilibrium distribution exists if and only if
  # there is an ergodic (aperiodic + positive recurrent) class C
  # and the Markov chain is certain to be absorbed into C eventually
  # (only C is colsed), wherever it starts.
  #
  # we are currently only dealing with finite discrete time M.C..
  # A finite closed class must be recurrent
  # and a finite Markov chain cannot contain any null recurrent states.
  #
  # so, in our case, closed class must be positive recurrent.
  # i.e. we can decide whether a class is positive recurrent
  # from whether it is closed. (Am I right? still thinking......)
  # (This is just a temporary method,
  # and I am still thinking about how to create a function to distinguish
  # positive recurrent, null recurrent and transient)


  decide <- dmc_irreclass(MC)

  # if total number of closed class is not exactly 1,
  # no equalibrium distribution.
  if ((sum(decide[,1] == "A closed class")) != 1) {
    return("There is no equilibrium distribution of this discrete-time M.C.")
  }

  # if the unique closed class is not aperiodic, no equalibrium distribution.
  for (i in 1 : length(decide[,1])) {
    if (decide[,1][i] == "A closed class" & decide[,2][i] != "aperiodic") {
      return("There is no equilibrium distribution of this discrete-time M.C.")
    }
  }

  # If there is only one closed class and the closed class is ergodic,
  # an equilibrium distribution exists.

  state_space <- nrow(MC@p)
  random <- runif(1, min = 0, max = 1)  # a random number in [0, 1]
  # Create a random initial distribution
  vector <- matrix(c(random, 1 - random,
                     runif(state_space - 2, min = 0, max = 0)),
                   nrow = 1, ncol = state_space)
  equ_distri1 <- (vector %*% MC@p)

  for (i in 1:(6000 * state_space)) {
    equ_distri1 <- equ_distri1 %*% MC@p
  }

  # Actually we don't have to calculate equ_distri2,
  # beacause we have already checked this M.C. in previous steps.
  # Keeping this just for double check.
  random <- runif(1, min = 0, max = 1)
  vector <- matrix(c(random,
                     runif(state_space - 2, min = 0, max = 0), 1 - random),
                   nrow = 1, ncol = state_space)

  equ_distri2 <- (vector %*% MC@p)
  for (i in 1:(6000 * state_space + 1)) {
    equ_distri2 <- equ_distri2 %*% MC@p
  }

  if (all.equal(equ_distri1, equ_distri2) == T) {
    # Very small numbers are rounded down to zero by using the function round
    cat("Equilibrium distribution is :", round(equ_distri2, digits = 10000000))
  } else {
    print("There is no equilibrium distribution of this discrete-time M.C.")
  }
}