R/bd_equi.R

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

#' Equilibrium distribution of the jump chain of a
#' linear birth-death-immigration-emigration process
#'
#' \code{bd_equi} is used to find the equilibrium distribution of the
#' jump chain of the input process.
#'
#' In continuous time M.C..
#' If {X (t), t >= 0} is an irreducible continuous-time Markov chain.
#' If there exists an invariant distribution pi then it is unique,
#' and the pi is also the equilibrium distribution of the chain.
#' If there is no invariant distribution,
#' then no equilibrium distribution exists.
#' And When a process is not irreducible, cannot guarantee
#' that an invariant distribution is an equilibrium distribution.
#'
#' @param BD is an object in class 'stat2003.bd'
#' @param showstate If there is so many state space, and you just want to
#' take first n elements of the equilibrium distribution,
#' then setting showstate = n. If showstate = 0, \code{bd_equi} will
#' return whole equilibrium distribution. Default value of showstate is zero.
#'
#' @return \code{bd_equi} will return a vector if the process
#' has equilibrium distribution.
#'
#' @seealso \code{\link{stat2003.bd-class}} A class type of linear
#' birth-death in stat2003 package.
#' @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.
#' @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_equi(A, 10)
#'
#'
#'@export
bd_equi <- 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 than nstate!")
  }

  lb <- BD@lb
  ld <- BD@ld
  im <- BD@im
  em <- BD@em
  nstate <- BD@nstate

  # There are so many if-else in this function

  # In continuous time M.C..
  # If {X (t), t >= 0} is an irreducible continuous-time Markov chain.
  # If there exists an invariant distribution pi then it is unique,
  # and the pi is also the equilibrium distribution of the chain.
  # If there is no invariant distribution,
  # then no equilibrium distribution exists.
  # And When a process is not irreducible, cannot guarantee
  # that an invariant distribution is an equilibrium distribution.

  # In continuous time M.C..
  # we can find invariant distribution by using piQ = 0,
  # where Q is the transition rates matrix of the jump chain of the M.C..
  # And after calculation, suppose pi = (pi0, pi1, pi2, ....)
  # and denotes A = (product of birth rate from state 0 to i-1)/
  #                 (product of death rate from state 1 to i)
  # where birth rate includes linear birth and immigration,
  # death rate includes linear death and emigration.
  # pi(i) = Ai * pi0 and pi0 = (1 + sum from i = 1 to infinity Ai )^(-1)
  #
  # If (sum from i = 1 to infinity Ai) is divergent,
  # an invariant distribution (and therefore an equilibrium
  # distribution) cannot exist.
  # Note: the sum in this function always means
  # (sum from i = 1 to infinity Ai)
  # Note: If the state space is finite, then the sum must be finite.
  #
  # Will use mathod above to find invariant distribution,
  # and then find equilibrium distribution.

  # The very first thing in this function is to find which birth-death
  # processes are irreducible.
  # B denotes Linear birth, D denotes Linear death,
  # I denotes Immigration, E denotes Emigration.
  #
  # Irreducible: BDIE,
  #              BDI, BIE, DIE
  #              DI, IE
  #
  # Not irreducible: BDE
  #                  BE, BD, BI, DE
  #                  B, I, D, E
  #
  # i.e. if a process has no Immigration, it is not irreducible,
  # if a process has Immigration but has no both of Emigration
  # and Linear death it is not irreducible.
  irre <- NULL

  if ((im == 0) | (im != 0 & (ld == 0 & em == 0))) {
    irre <- FALSE
  } else {
    irre <- TRUE
  }

  if (irre) {
    # we still have to consider whether the state space is finite
    # if it is not (nstate = 0 means infinity state space),
    if (nstate == 0) {
      PI <- matrix(0, nrow = 1, ncol = showstate)
      # when th M.C. has infinity state space, we can ignore
      # immigration and emigration.
      # By ratio test,
      # if linear birth rate < linear birth rate,
      # the sum converges,
      # so it do has the invariant and equilibrium distribution.
      if (lb < ld) {

        sum <- 0
        product <- 1

        for (i in 1 : (1000 * nstate)) {
          product <- product * ((lb * (i - 1) + im) / (ld * i + em))
          sum <- sum + product
        }

        PI[1, 1] <- (1 + sum)^(- 1)

        product <- 1
        for (i in 1 : (showstate - 1)) {
          product <- product * ((lb * (i - 1) + im) / (ld * i + em))
          PI[1, 1 + i] <- PI[1, 1] * product
        }

        } else if (lb == ld) {
        # if linear birth rate = linear death rate,
        # we have to consider im and em rate
          if (im < em) {
          # if im rate < em + ld rate,
          # the sum converges,
          sum <- 0
          product <- 1

          for (i in 1 : (1000 * nstate)) {
            product <- product * ((lb * (i - 1) + im) / (ld * i + em))
            sum <- sum + product
          }

          PI[1, 1] <- (1 + sum)^(- 1)

          product <- 1
          for (i in 1 : (showstate - 1)) {
            product <- product * ((lb * (i - 1) + im) / (ld * i + em))
            PI[1, 1 + i] <- PI[1, 1] * product
          }

        } else {
          # if im rate >= em rate, the sum diverges,
          # so no invariant or equilibrium distribution.
          return("No equilibrium distribution of the continuous-time M.C.")
        }

      } else {
        # if linear birth rate > linear death rate,
        # the sum diverges.

        return("No equilibrium distribution of the continuous-time M.C.")
      }
    } else {
      # if nstate != 0, i.e. it has finity state space.
      # then the sum must be finite.
      # then use pi(i) = Ai * pi0 and
      # pi0 = (1 + sum from i = 1 to infinity Ai )^(-1) directly.

      sum <- 0
      product <- 1
      PI <- matrix(0, nrow = 1, ncol = nstate)

      for (i in 1 : (nstate - 1)) {
        product <- product * ((lb * (i - 1) + im) / (ld * i + em))
        sum <- sum + product
      }

      PI[1, 1] <- (1 + sum)^(- 1)

      product <- 1
      for (i in 1 : (nstate - 1)) {
        product <- product * ((lb * (i - 1) + im) / (ld * i + em))
        PI[1, 1 + i] <- PI[1, 1] * product
      }

    }

  } else {
    # if not irreducible
    if (nstate == 0){
      # if the M.C. has infinity state space.
      PI <- matrix(0, nrow = 1, ncol = showstate)

      if (lb == 0 & im == 0) {
        # if the process only has death rate (including ld and em)
        # the state 0 wil be an absorbing state,
        # and therefore equilibrium distribution = (1, 0, 0, ...)
        PI[1, 1] <- 1

      } else if (ld == 0 & em == 0) {
        # if the process only has birth rate (including lb and im)
        # the M.C. diverges. (infinity state space.)
        return("No equilibrium distribution of the continuous-time M.C.")

      } else {
        if (lb <= ld) {
          # if lb rate <= ld rate,
          # the state 0 wil be an absorbing state,
          # and therefore equilibrium distribution = (1, 0, 0, ...)
          PI[1, 1] <- 1
        } else {
          # if lb rate > ld rate, the M.C. diverges.(infinity state space)
          return("No equilibrium distribution of the continuous-time M.C.")
        }

      }
    } else {
      #if the M.C. has finite state space,
      PI <- matrix(0, nrow = 1, ncol = nstate)

      if (lb == 0 & im == 0) {
        # if the process only has death rate (including ld and em)
        # the state 0 wil be an absorbing state,
        # and therefore equilibrium distribution = (1, 0, 0, ...)
        PI[1, 1] <- 1

      } else if ((lb != 0 & im != 0 & ld == 0 & em == 0) |
                 (lb == 0 & im != 0 & ld == 0 & em == 0)) {
        # if the process only has im rate or im and lb rate
        # the last state wil be an absorbing state,
        # and therefore equilibrium distribution = (0, 0, ..., 1)
        PI[1, nstate] <- 1

      } else if (lb != 0 & im == 0 & ld == 0 & em == 0) {
        # if the process only has lb rate
        # state 0 and the last state wil be two absorbing states,
        # so no equilibrium distribution
        return("No equilibrium distribution of the continuous-time M.C.")

      } else {
        if (lb <= ld) {
          # if lb rate <= ld rate,
          # the state 0 wil be an absorbing state,
          # and therefore equilibrium distribution = (1, 0, 0, ...)
          PI[1, 1] <- 1

        } else {
          # if lb rate > ld rate, the M.C. diverges.
          return("No equilibrium distribution of the continuous-time M.C.")
        }

      }
    }
  }

  rownames(PI) <- " "
  colnames(PI) <- 0 : (ncol(PI) - 1)

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

}