R/sensitivity.R

In markovchain: Easy Handling Discrete Time Markov Chains

#' Sensitivity of the stationary distribution to a state's transition row
#'
#' Computes, for a finite, irreducible discrete-time Markov chain, the
#' coefficients needed to obtain the first-order change in the stationary
#' distribution caused by an infinitesimal perturbation of one row of the
#' transition matrix.
#'
#' A single entry \eqn{p_{kl}} of a stochastic matrix cannot be perturbed on
#' its own without leaving row \eqn{k}: some other entry (or entries) of that
#' row must move to compensate, so that the row still sums to one. Any
#' admissible perturbation of row \eqn{k} is therefore a direction vector
#' \eqn{d\in\mathbb{R}^n} with \eqn{\sum_l d_l = 0}, giving the perturbed
#' matrix \eqn{P(\varepsilon) = P + \varepsilon\, e_k d^{\mathsf T}} for small
#' \eqn{\varepsilon}. This function returns the \eqn{n\times n} matrix
#' \eqn{S} such that, for \emph{every} such \eqn{d} and every state \eqn{j},
#' \deqn{\left.\frac{d\pi_j}{d\varepsilon}\right|_{\varepsilon=0} =
#'   \sum_l d_l\, S_{lj} = \left(d^{\mathsf T} S\right)_j.}
#'
#' @param object A \code{markovchain} object representing a finite,
#'   irreducible discrete-time Markov chain.
#' @param state A single state name (character) or state index (single
#'   positive integer), identifying the row \eqn{k} of the transition matrix
#'   to be perturbed.
#'
#' @return An \eqn{n\times n} numeric matrix \eqn{S}, with both dimensions
#'   named after \code{states(object)}. Row \eqn{l} of \eqn{S} corresponds to
#'   the perturbation direction "increase \eqn{p_{kl}}" (paired with a
#'   compensating decrease elsewhere in row \eqn{k}); column \eqn{j}
#'   corresponds to the affected stationary probability \eqn{\pi_j}. See
#'   Details for how to read individual entries.
#'
#' @details
#' \strong{Closed form.} Let \eqn{Z=(I-P+\mathbf 1\pi^{\mathsf T})^{-1}} be
#' the fundamental matrix already used by \code{\link{kemenyConstant}}. Then
#' \deqn{S_{lj} = \pi_k\left(Z_{lj} - \pi_j\right).}
#' This particular centering (subtracting \eqn{\pi_j}, the same constant for
#' every row \eqn{l}) is what makes \eqn{S} usable directly with \emph{any}
#' zero-sum direction \eqn{d}, because \eqn{\sum_l d_l \pi_j = \pi_j\sum_l
#' d_l = 0} drops out of the sum above -- adding any other per-column
#' constant to \eqn{S} would give the same directional derivatives, but this
#' one has the convenient side effect that \code{sensitivity(object,
#' state)[state, ]} is the sensitivity of "leaving row \code{state}
#' unchanged", which is informative on its own (it need not be zero: the
#' *direction* \eqn{d=e_{\mathrm{state}}} is generally not itself a valid
#' zero-sum perturbation by itself, only differences of rows are).
#'
#' \strong{The common two-state case.} The usual textbook question --
#' "increase \eqn{p_{k,\mathrm{to}}} by \eqn{\varepsilon}, decrease
#' \eqn{p_{k,\mathrm{from}}} by \eqn{\varepsilon}, how does \eqn{\pi} move?"
#' -- is answered by taking the difference of two rows of \eqn{S}:
#' \deqn{\left.\frac{d\pi}{d\varepsilon}\right|_{\varepsilon=0} =
#'   S_{\mathrm{to}, \cdot} - S_{\mathrm{from}, \cdot}.}
#' See the second example below, which checks this against a direct
#' finite-difference recomputation of the stationary distribution.
#'
#' Only irreducibility is required, not aperiodicity: \eqn{Z} and \eqn{\pi}
#' are well defined for any irreducible chain regardless of periodicity.
#'
#' The implementation calls \code{\link{steadyStates}} once and then solves
#' one dense linear system for \eqn{Z}; both are \eqn{O(n^3)} time and
#' \eqn{O(n^2)} memory for a dense \eqn{n}-state transition matrix, the same
#' cost as \code{\link{kemenyConstant}}. It supports both row- and
#' column-stochastic storage; \eqn{S} is always returned with rows/columns
#' indexed by state name in the chain's own state order.
#'
#' @references
#' Schweitzer, P. J. (1968). Perturbation theory and finite Markov chains.
#' \emph{Journal of Applied Probability}, 5(2), 401-413.
#'
#' Meyer, C. D. (1980). The condition of a finite Markov chain and
#' perturbation bounds for the limiting probabilities. \emph{SIAM Journal on
#' Algebraic and Discrete Methods}, 1(3), 273-283.
#'
#' Cho, G. E. and Meyer, C. D. (2001). Comparison of perturbation bounds for
#' the stationary distribution of a Markov chain. \emph{Linear Algebra and
#' its Applications}, 335(1-3), 137-150.
#'
#' @seealso \code{\link{kemenyConstant}}, \code{\link{steadyStates}},
#'   \code{\link{is.irreducible}}
#'
#' @examples
#' statesNames <- c("a", "b", "c")
#' mc <- new("markovchain", states = statesNames,
#'   transitionMatrix = matrix(c(0.5, 0.3, 0.2,
#'                               0.2, 0.6, 0.2,
#'                               0.1, 0.1, 0.8), byrow = TRUE, nrow = 3,
#'                             dimnames = list(statesNames, statesNames)))
#' S <- sensitivity(mc, "a")
#' S
#'
#' # Check against a finite-difference recomputation of the stationary
#' # distribution: increase p("a"->"c") and decrease p("a"->"b") by eps.
#' eps <- 1e-6
#' P2 <- mc@transitionMatrix
#' P2["a", "c"] <- P2["a", "c"] + eps
#' P2["a", "b"] <- P2["a", "b"] - eps
#' mc2 <- new("markovchain", states = statesNames, transitionMatrix = P2)
#' (steadyStates(mc2) - steadyStates(mc)) / eps   # finite difference
#' S["c", ] - S["b", ]                            # closed-form prediction
#'
#' @exportMethod sensitivity
setGeneric("sensitivity", function(object, state) standardGeneric("sensitivity"))

#' @rdname sensitivity
setMethod("sensitivity", "markovchain", function(object, state) {
  if (!is.irreducible(object)) {
    stop("sensitivity is defined here only for irreducible Markov chains.")
  }

  stateNames <- states(object)
  n <- length(stateNames)

  if (length(state) != 1L || anyNA(state)) {
    stop("state must be a single state name or a single state index.")
  }
  if (is.numeric(state)) {
    if (state != as.integer(state) || state < 1L || state > n) {
      stop("state, if numeric, must be a single integer between 1 and the number of states.")
    }
    k <- as.integer(state)
  } else {
    k <- match(state, stateNames)
    if (is.na(k)) {
      stop("Unknown state: ", state)
    }
  }

  P <- as.matrix(object@transitionMatrix)
  if (!object@byrow) {
    P <- t(P)
  }
  if (n != ncol(P) || any(!is.finite(P))) {
    stop("The transition matrix must be square and finite.")
  }

  pi <- as.numeric(steadyStates(object))
  if (length(pi) != n || any(!is.finite(pi)) || sum(pi) <= 0) {
    stop("Unable to obtain a valid stationary distribution.")
  }
  pi <- pi / sum(pi)

  # Same fundamental-matrix construction used by kemenyConstant(): avoids
  # explicitly forming the n-by-n outer product 1 %*% t(pi).
  A <- diag(n) - P
  A <- sweep(A, 2L, pi, FUN = "+")
  Z <- solve(A)

  # S[l, j] = pi[k] * (Z[l, j] - pi[j]); sweep() subtracts pi from every row
  # of Z (recycled across columns), then the whole matrix is scaled by
  # pi[k].
  S <- pi[k] * sweep(Z, 2L, pi, FUN = "-")

  if (!all(is.finite(S))) {
    stop("Unable to compute a finite sensitivity matrix.")
  }
  dimnames(S) <- list(stateNames, stateNames)
  S
})