R/dmc_tp.R

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

#' Transient probabilities of a discrete-time Markov chain
#'
#' \code{dmc_tp} can calculate transient probabilities at a specific step,
#' and also can give a graph about transient probabilities from
#' initial step to that specific step. The graph describes transient probabilities
#' in gray levels; zero indicates "white" and one indicates "black".
#'
#'
#' @param MC is an object in class 'stat2003.d'
#' @param time an interger, total number of steps
#' that the process will run in this simulation.
#' @param graph Logical. If True the function will also return a graph about
#' transient probabilities, otherwise the function will not return the graph.
#'
#' @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_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_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.
#'
#' @return \code{dmc_tp} returns transient probabilities at a specific step
#' for each state and can also give a graph about transient probabilities
#' from step zero to that specific step.
#'
#'
#' @examples
#'
#'
#'
#' 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(1, 0, 0, 0, 0, 0, 0), p = m,
#'          statespace = c("A", "B", 2 : 6) )
#'
#' dmc_tp(A, time = 9, graph = TRUE)
#'
#'
#'
#'
#' 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), nrow = 4, ncol=4, byrow = TRUE)
#' p_start <- c(0,0,0,1)
#'
#' A <- new("stat2003.d", p_start = c(0, 0, 0, 1), p = m,
#'          statespace = c("A", "B", "C", "D") )
#'
#' dmc_tp(A, time = 10, graph = TRUE)
#'
#'
#'
#'@export
dmc_tp <- function(MC, time, graph = F) {
  if (class(MC) != "stat2003.d") {
    stop("MC has to be an object in class 'stat2003.d' !")
  }
  if (time %% 1 != 0) {
    stop("time must be an integer!")
  }

  tp <- MC@p_start %*% MC@p
  # we have to use rev(), beacause order will be different
  # when we change tp_ap to a matrix.
  tp_ap <- append (rev(MC@p_start), rev(tp))

  for (i in 1 : (time - 1)) {
    tp <- tp %*% MC@p
    tp_ap <- append (tp_ap, rev(tp))
  }

  if (graph == T) {
    # want to make zero indicates "white" and
    # one indicates "black", so 1 - tp_ap.
    tp_ap <- 1 - tp_ap
    tp_ap <- matrix((grDevices::grey(tp_ap)), nrow = nrow(MC@p), ncol = time + 1)
    image <- grDevices::as.raster(tp_ap)

    graphics::plot(c(0, time + 1), c(0, nrow(MC@p)), type = "n", xaxt="n",
         yaxt = "n", xlab = "Steps", ylab = "State Space")
    graphics::rasterImage(image, 0, 0, time + 1, nrow(MC@p), interpolate = FALSE)
    graphics::axis(side = 2, at = ((0 : nrow(MC@p)) + 0.5),
         labels = c(MC@statespace, 1), las = 2)
    graphics::axis(side = 1, at = ((0 : (time)) + 0.5), labels=(0 : (time)))

    print(matrix(tp, nrow = 1, ncol = nrow(MC@p),
                 dimnames = list("Probability", MC@statespace)))
  } else {
    print(matrix(tp, nrow = 1, ncol = nrow(MC@p),
                 dimnames = list("Probability", MC@statespace)))
  }
}