R/ExpectedValue.R

In feralaes/dampack: Package with useful decision-analytic modeling tools

#' A Markov function to compute expected value
#'
#' Computes the expected value from a Markov transition array
#' @param trans Markov transition array generated by `CalculateMarkovTrace` function
#' @param rwd Matrix or array with state rewards (diagonal) and transition rewards
#         (off-diagonal)
#' @param r Discount factor in the scale of a cycle [0, 1]; default = 0
#' @param half If TRUE apply Half cycle correction. Defaults to TRUE
#' @param verbose If TRUE messages are printed. Defaults to FALSE
#' @keywords Markov
#' @return exp.value Expected value
#'
ExpectedValue <- function(trans, rwd, r = 0, half = TRUE, verbose = FALSE){
  # Extract name of states
  state.names <- colnames(rwd)
  # Number of states
  n.states <- dim(trans)[1]
  # Number of cycles
  n.cycles <- dim(trans)[3]
  # Verify if rwd is 2D or 3D matrix
  n.dim <- length(dim(rwd))
  # Extract name of states
  #state.names <- colnames(trans)
  # If M is a 2D Matrix, convert to 3D by repeating it n.cycles times in a 3D array
  if (n.dim < 3){  
    rwd <- array(rep(rwd, n.cycles), 
                 dim = c(n.states, n.states, n.cycles), 
                 dimnames = list(state.names, state.names, 
                                 paste("Cycle", 0:(n.cycles-1), sep = "")))
    if(verbose) print("Reward matrix input by user is 2D, 3D stacking performed")
  } else {
    if(verbose) print("Reward matrix input by user is 3D, no stacking perfomed")
  }
  
  # Discount vector
  disc <- 1/(1 + r)^seq(0, (n.cycles-1))
  # Apply half cycle correction
  if (half){
    trans[, , 1]    <- trans[, , 1]*0.5
    trans[, , n.cycles] <- trans[, , n.cycles]*0.5
  }
  #   trans.rwd <- array(apply(trans, 3, function(x) x*rwd), 
  #                      dim = c(n.states, n.states, n.cycles))
  trans.rwd <- trans*rwd
  trace.rwd <- t(colSums(trans.rwd))*disc
  
  exp.value <- sum(trace.rwd)
  
  return(exp.value)
}