stat2003: stat2003: Introduction to Applied Probability at UCL

Documented in dmc_inv

#' Invariant distributions of a discrete-time finite statespace Markov chain
#'
#' \code{dmc_inv} can calculate invariant distributions of the input chain.
#' Moreover, for user's convenience, \code{dmc_inv} is also empowered to
#'calculate the period for each closed class,
#'which is similar to a part of function \code{dmc_irreclass}.
#'
#' @param MC is an object in class 'stat2003.d'
#' @param show_mat A logical scalar.  If TRUE (the default) the transition
#'   matrix (rounded to 1 decimal places) is printed to screen so that the
#'   user can check it.  Otherwise (show_mat = FALSE) the matrix is not
#'   printed.
#' @param show_closed A logical scalar.  If TRUE the number of colsed class
#'   and the details of each closed class is printed to screen so that the
#'   user can check it.  Otherwise (show_closed = FALSE) the classes is not
#'   printed.
#' @param lambda_tol A numeric scalar.  Tolerance for finding eigenvalues of
#'   probability transition matrix whose real part is equal to 1.  Any eignvalue
#'   within a distance lambda_tol of 1 is treated as being equal to 1.
#' @param imag_tol A numeric scalar.  Tolerance for checking that the imaginary
#'   part of an eigenvalue is close enough to zero.  This is only used when the
#'   real part of the eigenvalue is close enough to 1.  Also used to check that
#'   the imaging parts of the corresponding eigenvectors are close enough to
#'   zero.
#' @param row_tol A numeric scalar.  Tolerance for checking that the rows of
#'   the input probability transition matrix sum to 1.  Any matrix that has any
#'   row sums more than row_tol away from 1 is rejected.
#' @param too_big_to_print A numeric scalar.  If the statespace is bigger than this,
#'   the transition matrix will also not be printed.
#'
#' @return An n_inv by n_s numeric matrix, where n_stat is the number of
#'   invariant distributions and n_s is the size of the state space, i.e.
#'   the dimension of the symmetric input transition matrix.
#'   Therefore, each ROW of the matrix contains an invariant distribution of
#'   the Markov chain.
#'
#' @seealso \code{\link{stat2003.d-class}}  A class type of discrete-time finite
#' state space Markov chain in stat2003 package.
#' @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_equi}} returns the equilibrium
#' distribution for a discrete-time Markov chain.
#' @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(1/3, 2/3, 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_inv(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), nr = 6, nc=6, byrow = TRUE)
#' A <- new("stat2003.d", p_start = rep(1/6, 6), p = m,
#'          statespace = letters[1 : 6] )
#' dmc_inv(A)
#'
#'
#'@export
dmc_inv <- function(MC, show_mat = TRUE, show_closed = TRUE, lambda_tol = 1e-4,
                      imag_tol = 1e-4, too_big_to_print = 15, row_tol = 1e-4 ) {
  if (class(MC) != "stat2003.d") {
    stop("MC has to be an object in class 'stat2003.d' !")
  }
  P_dim <- dim(MC@p)

  if (show_mat & P_dim[1] < too_big_to_print) {
    cat("Finding invariant distributions of the discrete-time Markov chain
        with transition matrix", "\n")
    print(round(MC@p,2))
  }
  # just save p of MC to a new variable.
  MCp <- MC@p
  prob_tol <- 1e-4
  # Find the (left) eigenvalues and eigenvectors of MC@p, i.e.
  # any pairs of vector vec and scalar lambda such that
  # vec %** MC@p = lambda * vec.
  # If there are any such combinations for which lambda = 1 then
  # vec (after normalisation to sum to 1 so that it is a probability
  # vector) satisifies vec = vec P and therefore is an invariant
  # distribution of the chain.
  #
  # Find the (left) eigenvalues and eigenvectors.
  # We do this by finding the eigen decomposition of the transpose
  # of the transition matrix.  We would like to find the solutions
  # of pi = pi P that correspond to the closed classes of the Markov
  # Chain, that is, solutions for which only states in a particular
  # closed class have non-zero probability.  To do this we proceed
  # iteratively: if P has more than one eigenvalue that is equal to 1
  # then we
  # (a) take the first such eigenvalue,
  # (b) identify the states for which the probability in the corresponding
  #     eigenvector is non-zero, and
  # (c) remove these states from the Markov chain.
  # Then we repeat the eigen decomposition and steps (a) - (c) until the
  # are no more eigenvalues that are equal to 1.
  #
  # Eigen decomposition for the whole chain.
  eigen_list <- eigen(x = t(MC@p))
  # Function to find all the eigenvalue-eigenvector combinations such that
  # 1. The real part of eigenvalue is close enough to 1.
  # 2. The imaginary part of the eigenvalue is close enough to 0.
  # 3. The imaginary parts of the eigenvector are all close enough to 0.
  e_value_eq_one_check <- function(x_list) {
    x <- x_list$values
    real_cond <- abs(Re(x) - 1) < lambda_tol
    imag_cond <- abs(Im(x) - 0) < imag_tol
    y <- x_list$vectors
    vec_cond <- apply(y, 2, function(x) all(abs(Im(x) - 0) < imag_tol))
    which(real_cond & imag_cond)
  }
  # Which eigenvalues are close enough to 1, i.e. within lambda_tol?
  eigen_eq_one <- e_value_eq_one_check(eigen_list)
  n_ones <- length(eigen_eq_one)
  # If there are no unit eigenvalues then stop.
  if (n_ones == 0) {
    stop("No invariant distributions have been found")
  }
  inv_distns <- NULL
  # For each unit eigenvalue calculate the number of (essentially)
  # non-zero values in the corresponding eignvector.
  temp <- Re(eigen_list$vectors[, eigen_eq_one, drop = FALSE])
  n_non_zero <- apply(temp, 2, function(x) sum(abs(x) > prob_tol))
  which_e_value <- which.min(n_non_zero)
  # Find the first eigenvector with unit eigenvalue.
  first_vec <- Re(eigen_list$vectors[, eigen_eq_one[which_e_value]])
  # Set values close enough to zero to zero.
  first_vec[abs(first_vec) < prob_tol] <- 0
  # Normalise
  first_vec <- first_vec / sum(first_vec)
  # Add the vector to an inv_distns matrix.
  inv_distns <- rbind(inv_distns, first_vec)
  # If there is not exactly one unit eigenvalue.
  if (n_ones != 1) {
    # If we get to here then n_ones > 1.
    removed_states <- NULL
    remaining_states <- 1:P_dim[1]
    while (n_ones > 1 & length(remaining_states) > 0) {
      # Which states does we remove?
      rm_states <- remaining_states[which(first_vec != 0)]
      # Keep track of which states have been removed.
      removed_states <- c(removed_states, rm_states)
      # ... and which states remain.
      remaining_states <- remaining_states[!(remaining_states %in% rm_states)]
      # Remove states.
      MC@p <- MCp[remaining_states, remaining_states, drop = FALSE]
      # Rescale the transition matrix of the remaining states so that
      # the rows sum to 1.
      MC@p <- MC@p / rowSums(MC@p)
      # New eigen decomposition.
      eigen_list <- eigen(x = t(MC@p))
      eigen_eq_one <- e_value_eq_one_check(eigen_list)
      temp <- Re(eigen_list$vectors[, eigen_eq_one, drop = FALSE])
      n_non_zero <- apply(temp, 2, function(x) sum(abs(x) > prob_tol))
      which_e_value <- which.min(n_non_zero)
      first_vec <- Re(eigen_list$vectors[, eigen_eq_one[which_e_value]])
      first_vec[abs(first_vec) < prob_tol] <- 0
      first_vec <- first_vec / sum(first_vec)
      # Create a vector with zeros for the removed states and
      # the values in first_vec for the remaining states.
      inv_vec <- numeric(P_dim[1])
      inv_vec[remaining_states] <- first_vec

      # If this invariant distribution only contains 1, like (0, 0, 0, 1)
      # We have to go back to initial transition matrix to check whether
      # the state containing 1 is an absorbing state.
      # Because, sometimes, the previous steps (deleting and normalising)
      # will change the probability from others to 1.
      if (1 %in% inv_vec) {
        location <- match(1, inv_vec)
        if (MCp[location, location] == 1) {
          inv_distns <- rbind(inv_distns, inv_vec)
          n_ones <- length(eigen_eq_one)
        } else {
          # if it is not a real abosrbing state
          # just don't record it.
          n_ones <- length(eigen_eq_one)
        }
      } else {
        # checking whether the closed class is a fake closed class
        sum = 0
        for (k in which(inv_vec != 0)) {
          for (z in which(inv_vec == 0)) {
            if (MCp[k, z] == 0) {
              # come back to original transition matrix to check
              # whether the proabailities of other states which are
              # in the different class are equal to zero.
              sum = sum + 1
            }
          }
        }

        number = length(which(inv_vec != 0)) * length(which(inv_vec == 0))
        if (sum == number) {
          inv_distns <- rbind(inv_distns, inv_vec)
          n_ones <- length(eigen_eq_one)
        } else {
          # if it is not a real closed class
          # just don't record it.
          n_ones <- length(eigen_eq_one)
        }

      }
    }
  }
  # Return a matrix with an invariant distribution in each row.
  rownames(inv_distns) <- paste("inv_dist", 1:nrow(inv_distns))
  colnames(inv_distns) <- MC@statespace


  # if show_closed = TRUE, we have to find closed classes
  # and period of each class.
  number_closed <- nrow(inv_distns)
  if (show_closed == TRUE) {
    cat("There is/are ", number_closed, "closed class in this Markov Chain.\n")
    for (i in 1 : number_closed) {
      state <- which(inv_distns[i, ] != 0)
      # we changed MC@p in previous steps, so we have to use MCp now.
      closed_p <- MCp[state, state]
      # if this is not an absorbing state, normalise it.
      # if it is and have invariant dis.,
      # there is no need to normalise and calculate
      # because the probability must be 1 and period must be 1.
      if (length(state) != 1){
        for (j in 1 : nrow(closed_p)) {
          closed_p[j, ] <- closed_p[j, ] / sum(closed_p[j, ])
        }
        # Point (ii) at the bottom of page 123 of book
        # Theory of Stochastic Processes says if a chain is
        # irreducible and periodic with period t, the transition
        # matrix P has exactly t eigenvalues of unit modulus.
        # We can apply this to out function,
        # cause finite closed class is irreducible.
        closed_eigen <- eigen(closed_p)
        # finding eigenvalues of unit modulus
        modulus <- Mod(closed_eigen$values)
        period <- sum(mapply(function(x, y) {isTRUE(all.equal(x, y))},
                             modulus, 1))
        cat("Closed class", i, "with states: ", MC@statespace[state],
            " with period ", period, ".\n")
      } else {
        cat("Closed class", i, "with states: ", MC@statespace[state],
            " with period 1.\n")
      }
    }
  }
  return(inv_distns)
}