R/bd_matrix.R

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

#' Generator matrix for a birth-death-immigration-emigration process
#'
#' \code{bd_rate} is used to produce generator matrices (rate matrices)
#' of the input process.
#'
#' @param BD is an object in class 'stat2003.bd'
#' @param showstate If there are so many statespace, and users just want to
#' take first n elements of the generator matrix, setting showstate = n.
#' If showstate = 0, \code{bd_rate} will return whole
#' generator matrix. Default value of showstate is zero.
#'
#' @return \code{bd_rate} will return the generator matrix by given rates.
#'
#' @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_trans}} will return the transition matrix by
#' given birth-death rates for a linear birth death process.
#' @seealso \code{\link{bd_simu}} simulates a finite statespace Birth-death
#' process by returning a plot about states and time.
#'
#' @examples
#' A <- new("stat2003.bd", lb = 1, ld = 1, im = 3, em = 4,
#' nstate = 0)
#'
#' bd_rate(A, 4)
#'
#'@export
bd_rate <- function (BD, showstate = 0) {
  if (class(BD) != "stat2003.bd") {
    stop("BD has to be an object in class 'stat2003.bd' !")
  }

  if (BD@nstate == 0 & showstate == 0) {
    stop("nstate and showstate cannot be both zero!")
  }
  if (BD@nstate != 0 & BD@nstate < showstate) {
    stop("showstate cannot larger then nstate!")
  }

  if (BD@nstate != 0) {
    q <- matrix(0, nrow = BD@nstate,
                   ncol = BD@nstate)

    # death rate
    for (i in 2 : nrow(q)) {
      q[i, i - 1] <- (i - 1) * BD@ld + BD@em
    }
    # birth rate
    for (i in 2 : nrow(q)) {
      q[i - 1, i] <- (i - 2) * BD@lb + BD@im
    }
    # diagonal
    for (i in 1 : nrow(q)) {
      q[i, i] <- - sum(q[i,])
    }

  } else {
    q <- matrix(0, nrow = showstate + 1,
                   ncol = showstate + 1)

    # death rate
    for (i in 2 : nrow(q)) {
      q[i, i - 1] <- (i - 1) * BD@ld + BD@em
    }
    # birth rate
    for (i in 2 : nrow(q)) {
      q[i - 1, i] <- (i - 2) * BD@lb + BD@im
    }
    # diagonal
    for (i in 1 : nrow(q)) {
      q[i, i] <- - sum(q[i,])
    }
  }

  rownames(q) <- paste(0 : (nrow(q) - 1))
  colnames(q) <- paste(0 : (nrow(q) - 1))

  if (showstate == 0) {
      return(q)
  } else {
    return(q[1 : showstate, 1 : showstate])
  }
}


#' Transition matrix for the
#' jump chain of a birth-death-immigration-emigration process
#'
#' \code{bd_trans} is used to produce transition matrices of the
#' jump chain of the input process.
#'
#' @param BD is an object in class 'stat2003.bd'
#' @param showstate If there are so many statespace, and users just want to
#' take first n elements of the generator matrix, setting showstate = n.
#' If showstate = 0, \code{bd_trans} will return whole
#' transition matrix. Default value of showstate is zero.
#'
#' @return \code{bd_trans} will return the transition matrix by given rates.
#'
#' @seealso \code{\link{stat2003.bd-class}} A class type of linear
#' birth-death in stat2003 package.
#' @seealso \code{\link{bd_equi}} will returns 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_simu}} simulates a finite statespace Birth-death
#' process by returning a plot about states and time.
#'
#' @examples
#' A <- new("stat2003.bd", lb = 1, ld = 1, im = 3, em = 4,
#' nstate = 0)
#'
#' bd_trans(A, 4)
#'
#'
#'@export
bd_trans <- function(BD, showstate = 0) {
  if (class(BD) != "stat2003.bd") {
    stop("BD has to be an object in class 'stat2003.bd' !")
  }

  if (BD@nstate == 0 & showstate == 0) {
    stop("nstate and showstate cannot be both zero!")
  }
  if (BD@nstate != 0 & BD@nstate < showstate) {
    stop("showstate cannot larger then nstate!")
  }

  q <- bd_rate(BD, showstate)
  p <- q

  for (i in 1 : nrow(p)) {
    sum <- - sum(p[i, i])
    for (j in 1 : nrow(p)) {
      p[i, j] <- p[i, j] / sum
    }
  }

  # setting diagonal to zero
  for (i in 1 : nrow(p)) {
    p[i, i] = 0
  }


  # if there are any aborbing state
  for (i in 1 : nrow(q)) {
    if (sum(q[i, ] == 0) == nrow(q)) {
      p[i, ] <- 0
      p[i, i] <- 1
    }
  }

  return(p)
}