R/bd_simu.R

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

#' Simulation of a linear birth-death-immigration-emigration process
#'
#' \code{bd_simu} simulates afinite statespace birth-death
#' process by returning a possible sequence
#' and a plot about states and time.
#'
#' @param BD is an object in class 'stat2003.bd'
#' @param time A numeric scalar and it is the total time that the
#' process will run in this simulation.
#' @param ini_state A string and represents initial state of the process.
#' @param label An integer used to specify the interval between two labels
#' on the y-axis.
#' @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{bd_simu} will returns a possible sequence and
#' a plot about states and time.
#'
#' @seealso \code{\link{stat2003.bd-class}} A class type of linear
#' birth-death in stat2003 package.
#' @seealso \code{\link{bd_equi}} will return the equilibrium
#' distribution of a birth-death process
#' @seealso \code{\link{bd_rate}} will return the generator matrix by
#' given birth-death rates for a linear birth death process.
#' @seealso \code{\link{bd_trans}} will return the transition matrix by
#' given birth-death rates for a linear birth death process.
#'
#' @examples
#' A <- new("stat2003.bd", lb = 2, ld = 3, im = 12, em = 10,
#'           nstate = 30)
#' bd_simu(A, time = 100,
#'         ini_state = 5, label = 3)
#'
#'@export
bd_simu <- function (BD,
                     time, ini_state, label, graph = TRUE) {
  if (BD@nstate != 0) {
    statespace <- 0 : (BD@nstate - 1)

    if (!(ini_state %in% statespace)) {
      stop("initial state is not in state space!")
    }
  }

  if (class(BD) != "stat2003.bd") {
    stop("BD has to be an object in class 'stat2003.bd' !")
  }

  if (time < 0) {
    stop("time must greater or equal to zero!")
  }

  if (label %% 1 != 0) {
    stop("label must be an integer!")
  }

  nstate <- BD@nstate

  if (nstate != 0) {
    nstate <- length(statespace)
    q <- bd_rate(BD, nstate)
    p <- bd_trans(BD, nstate)
    time_series <- 0
    time_sum <- 0
    random <- 0
    current_state <- which(statespace == ini_state)
    state_series <- current_state
    rate <- - q[current_state, current_state]

    while (time_sum < time ) {
      t <- (-(log(runif(1)) / rate))
      # if current_state is an absorbing state, i.e. rate is zero,
      # and t = (-(log(runif(1)) / rate)),
      # so t will be infinity, which means system will skip this
      # loop, and then absorbing situation will not be shown
      # observably in our simulation graph, so add one more point
      # manually.
      if (is.infinite(t)) {
        time_sum <- sum(time_sum, t)
        state_series <- append(state_series, current_state)
        time_series <- append(time_series, time)

      } else {

        time_sum <- sum(time_sum, t)
        time_series <- append(time_series, time_sum)

        random <- runif(1)
        if (current_state == 1) {
          current_state <- current_state + 1
        } else if (current_state == nstate) {
          current_state <- current_state - 1
        } else if (random < p[current_state, current_state - 1]) {
          current_state <- current_state - 1
        } else {
          current_state <- current_state + 1
        }

        state_series <- append(state_series, current_state)
        rate <- - q[current_state, current_state]
      }

    }

    if (graph == TRUE) {

      graphics::plot(time_series, state_series, type = "l", yaxt = "n",
           xlab = "Time", ylab = "States",
           ylim = c(1, nstate),
           xlim=c(0, time))

      if (label != 0) {
        label_series <- seq(from = 1, to = nstate,
                            by = label)
        graphics::axis(side = 2, at = c(label_series),
             labels = statespace[label_series], las = 2)
      } else {
        graphics::axis(side = 2, at = c(1 : nstate), labels = statespace, las = 2)
      }
      return(statespace[state_series])

    } else {
      return(statespace[state_series])
    }

  } else {
    # if nstate == 0
    time_series <- 0
    time_sum <- 0
    random <- 0
    current_state <- ini_state + 1

    state_series <- current_state
    rate <- ((current_state - 1) * BD@ld + BD@em +
               (current_state - 1) * BD@lb + BD@im)
    while (time_sum < time ) {
      t <- (-(log(runif(1)) / rate))
      # if current_state is an absorbing state, i.e. rate is zero,
      # and t = (-(log(runif(1)) / rate)),
      # so t will be infinity, which means system will skip this
      # loop, and then absorbing situation will not be shown
      # observably in our simulation graph, so add one more point
      # manually.
      if (is.infinite(t)) {
        time_sum <- sum(time_sum, t)
        state_series <- append(state_series, current_state)
        time_series <- append(time_series, time)

      } else {

        time_sum <- sum(time_sum, t)
        time_series <- append(time_series, time_sum)

        random <- runif(1)
        if (current_state == 1) {
          if (BD@im != 0) {
            current_state <- current_state + 1
          }
        } else if (random < ((current_state - 1) * BD@ld + BD@em) / rate ) {
          current_state <- current_state - 1
        } else {
          current_state <- current_state + 1
        }

        state_series <- append(state_series, current_state)

        if (current_state == 1) {
          rate <- 0 + BD@im
        } else {
          rate <- ((current_state - 1) * BD@ld + BD@em +
                     (current_state - 1) * BD@lb + BD@im)
        }
      }
    }

    if (graph == TRUE) {
      statespace <- 0 : (max(state_series) - 1)

      graphics::plot(time_series, state_series, type = "l", yaxt = "n",
           xlab = "Time", ylab = "States",
           ylim = c(1, length(statespace)),
           xlim=c(0, time))

      if (label != 0) {
        label_series <- seq(from = 1, to = length(statespace),
                            by = label)
        graphics::axis(side = 2, at = c(label_series),
             labels = statespace[label_series], las = 2)
      } else {
        graphics::axis(side = 2, at = c(1 : length(statespace)), labels = statespace, las = 2)
      }
      return(statespace[state_series])

    } else {
      return(statespace[state_series])
    }

  }

}