R/dmc_simu.R

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

#' Simulation of a discrete-time finite statespace Markov chain
#'
#' \code{dmc_simu} simulates a discrete-time finite statespace
#' Markov chain by returning a possible sequence
#' and a states against steps plot.
#'
#'
#' @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 states
#' against steps plot, otherwise the function will not return the plot
#'
#' @return \code{dmc_simu} returns a possible sequence for simulation of
#' a discrete-time Markov chain and a states against steps plot.
#'
#' @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_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(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 = c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6), p = m,
#'          statespace = letters[1 : 6] )
#'
#' dmc_simu(A, 100, graph = TRUE)
#'
#'
#'@export
dmc_simu <- function (MC, time, graph = TRUE) {
  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!")
  }

  stime <- 0
  random <- runif(1)  # a random number.
  state <- NULL

  while (random > 0) {
    stime <- stime + 1
    # for example, the process has two state space 0 and 1,
    # and Initial distribution is (0.5, 0.5), random number initially equal
    # to 0.65 > 0.5, which means process is in state 1 at step 0.
    # so we need to calculate after how many times, the random number
    # becomes negative, then the number of
    # times is the state that the process will go next.
    random <- random - MC@p_start[stime]
  }

  stime_pre <- stime  # stime_pre means the state at previous step,
  state <- stime_pre  # using state to record every state at each step.

  for (i in 1 : time) {
    stime <- 0
    random <- runif(1)
    while (random > 0) {
      stime <- stime + 1
      # the state st previous step is stime_pre, so coordinate is [stime_pre, ].
      random <- random - MC@p[stime_pre, stime]
    }
    # we moved to next step, so the state at previous step also changed.
    stime_pre <- stime
    state <- append(state, stime_pre)
  }

  state <- as.vector(state)

  if (graph == TRUE) {
    graphics::plot(0 : time, state, type = "o", yaxt = "n",
         xlab = "Steps", ylab = "States",
         ylim = c(1, nrow(MC@p)), yaxp=c(1, nrow(MC@p), nrow(MC@p) - 1))
    graphics::axis(side = 2, at = (1 : nrow(MC@p)), labels = MC@statespace, las = 2)
    return(MC@statespace[state])

  } else {
    return(MC@statespace[state])
  }
}