R/phase_type.R

In rivasiker/phasty: General-purpose phase-type functions

#' Generator Function for Phase-Type Distributions
#'
#' Description of the S3 classes \code{cont_phase_type}, \code{disc_phase_type},
#' \code{mult_cont_phase_type}, \code{mult_disc_phase_type} which represents the
#' different phase-type distributions.
#'
#' @details
#'
#' \code{phase_type} is the generator function for the four types of phase-type
#' classes, respectively univariate continuous or discrete and multivariate
#' continuous or discrete which inherits from \code{list}.
#' The class is generated by supplying a sub-intensity matrix and an optional
#' initial probability vector plus a reward matrix in the case of multivariate
#' phase-type.
#' If the initial probabilities are not specified, then the initial probability
#' will be \code{init_probs = c(1, 0, 0, ...)} with the same length as the
#' number of transient states.
#'
#' @param subint_mat is the square matrix containing the transition rates or
#'   probabilities between transient states for continous or discrete
#'   phase-type respectively.
#'   If the phase-type is continuous, the subintensity matrix diagonal should
#'   only contains negative values and the row sums should be non-positive.
#'   If the phase-type is discrete, the subintensity matrix should only contains
#'   values between 0 and 1.
#' @param init_probs a vector, a one-row matrix or \code{NULL} which gives the
#'   probabilities to start in each state. If init_probs is \code{NULL},
#'   the probability to start on the first state will be 1 and to start on any
#'   other state 0.
#' @param reward_mat is a matrix code{NULL}(default) where each row is a reward
#'   vector, and each column corresponds to a state. It should have the same
#'   number of columns as the length of the initial probabilities.
#' @param round_zero is an integer or \code{NULL}(default), which gives the
#'   decimal from which we should truncate the positive values of the
#'   subintensity matrix.
#'   It could be useful in the scenarios where there is a reward transformation
#'   leading to values with many decimals and potentially computational
#'   approximation and potentially to positive row sums in continuous phase-type
#'   .
#'
#' @usage phase_type(subint_mat = NULL, init_probs = NULL, reward_mat = NULL,
#'                   round_zero = NULL)
#'
#' @examples
#'
#' ##===========================##
#' ## For continuous univariate ##
#' ##===========================##
#'
#' subintensity_matrix <- matrix(c(-1.5, 0, 0,
#'                                1.5, -1, 0,
#'                                0, 1, -0.5), ncol = 3)
#' phase_type(subintensity_matrix)
#'
#' #---
#'
#' subintensity_matrix <- matrix(c(-1.5, 0, 0,
#'                                1.5, -1, 0,
#'                                0, 1, -0.5), ncol = 3)
#' initial_probabilities <- c(0.9, 0.1, 0)
#' phase_type(subintensity_matrix, initial_probabilities)
#'
#' ##=========================##
#' ## For discrete univariate ##
#' ##=========================##
#'
#' subintensity_matrix <- matrix(c(0.4, 0, 0,
#'                                0.24, 0.4, 0,
#'                                0.12, 0.2, 0.5), ncol = 3)
#' phase_type(subintensity_matrix)
#'
#' #---
#'
#' subintensity_matrix <- matrix(c(0.4, 0, 0,
#'                                0.24, 0.4, 0,
#'                                0.12, 0.2, 0.5), ncol = 3)
#' initial_probabilities <- c(0.9, 0.1, 0)
#' phase_type(subintensity_matrix, initial_probabilities)
#'
#' ##=============================##
#' ## For continuous multivariate ##
#' ##=============================##
#'
#' subintensity_matrix <- matrix(c(-3, 0, 0,
#'                                2, -2, 0,
#'                                0, 1, -1), nrow = 3, ncol = 3)
#' reward_matrix = matrix(sample(seq(0, 10, 0.1), 6), nrow = 3, ncol = 2)
#' phase_type(subintensity_matrix, reward_mat = reward_matrix)
#'
#' ##===========================##
#' ## For discrete multivariate ##
#' ##===========================##
#'
#' subintensity_matrix <- matrix(c(0.4, 0, 0,
#'                                0.24, 0.4, 0,
#'                                0.12, 0.2, 0.5), ncol = 3)
#' reward_matrix <- matrix(sample(seq(0, 10), 6), nrow = 3, ncol = 2)
#' phase_type(subintensity_matrix, reward_mat = reward_matrix)
#'
#' @export

phase_type <- function(subint_mat = NULL, init_probs = NULL,
                       reward_mat = NULL, round_zero = NULL) {
  #############
  # Check the conditions necessary for every phase-type distribution
  #############

  if (is.null(subint_mat)) {
    stop('Unable to construct the phase-type distribution.
         Please provide either the type or the subintensity matrix.')
  }

  if (!is.matrix(subint_mat)) {
    stop('The subintensity matrix should be a matrix.')
  } else {
    subint_mat <- matrix(as.numeric(subint_mat), ncol = ncol(subint_mat))
  }

  if (nrow(subint_mat) != ncol(subint_mat)){
    stop('Subintensity matrix should be a square numerical matrix.')
  }

  if (is.null(init_probs)) {
    init_probs <- matrix(c(1, rep(0, nrow(subint_mat) - 1)),
                         1, nrow(subint_mat))
    warning('The initial probability vector is automatically generated.\n')
  }

  if ((is.vector(init_probs) & is.atomic(init_probs)) | is.matrix(init_probs)) {
    init_probs <- as.numeric(init_probs)
    init_probs <- matrix(init_probs, nrow = 1)
    if (!is.null(round_zero)){
      if (round(round_zero) == round_zero) { # avoid positive value due to
        #approximation
        init_probs[init_probs > 0] <- trunc(init_probs[init_probs > 0] *
                                              10^round_zero) / 10^round_zero
        subint_mat[subint_mat > 0] <- trunc(subint_mat[subint_mat > 0] *
                                              10^round_zero) / 10^round_zero
        print(trunc(subint_mat[subint_mat > 0],
                    round_zero))
      }
    }


    if (nrow(subint_mat) != length(init_probs)) {
      stop('The length of the initial probabilities does not match the size of
           the subintensity matrix.')
    }

    if (sum(init_probs) == 0){
      warning('The sum of the inital probability is equal to 0 with a defect
              of 1.\n')
    }

    if (sum(init_probs) < 0 || sum(init_probs) > 1){
      stop('The total initial probability should be between 0 and 1')
    } else if (sum(init_probs < 0) != 0 || sum(init_probs > 1) != 0) {
      stop('Each initial probability should be between 0 and 1.')
    }

  } else {
    stop('The initial probabilities must be a a matrix with one row or
         a vector.')
  }

  if (is.matrix(reward_mat)){
    if (sum(reward_mat < 0) < 0){
      stop('The reward matrix should only contains non-negative values.')
    }

    if (nrow(reward_mat) != length(init_probs)){
      stop('The reward matrix does not have the same number of columns as the
           number of states.')
    }
  }

  #############
  # Check the conditions necessary for continuous phase-type distribution
  #############

  if (length(which(diag(subint_mat) < 0)) == length(init_probs)) {
    # check if any rowsums are postive
    if (any(rowSums(subint_mat)>1e-14)){
      stop('The row sums of the subintensity matrix should be non-positive.')
    }
    if (is.matrix(reward_mat)){
      value <- list(subint_mat = subint_mat,
                    init_probs = init_probs,
                    reward_mat = reward_mat,
                    defect = 1 - sum(init_probs))
      attr(value, "class") <- "mult_cont_phase_type"
      return(value)
    } else {
      value <- list(subint_mat = subint_mat,
                    init_probs = init_probs,
                    defect = 1 - sum(init_probs))
      attr(value, "class") <- "cont_phase_type"
      return(value)
    }

    #############
    # Check the conditions necessary for discrete phase-type distribution
    #############

  } else if (sum(subint_mat < 0) == 0){
    if (sum(rowSums(subint_mat > 1)) > 0){
      stop('The rowsums should be between 0 and 1.')
    }

    if (sum(subint_mat > 1) > 0){
      stop('The subintensity matrix should only contains values between
           0 and 1.')
    }

    if (is.matrix(reward_mat)){
      if (sum(trunc(reward_mat)) != sum(reward_mat)){
        stop('The reward matrix should only contains integers.')
      }
      value <- list(subint_mat = subint_mat,
                    init_probs = init_probs,
                    reward_mat = reward_mat,
                    defect = 1 - sum(init_probs))
      attr(value, "class") <- "mult_disc_phase_type"
      return(value)
    } else {
      value <- list(subint_mat = subint_mat,
                    init_probs = init_probs,
                    defect = 1 - sum(init_probs))
      attr(value, "class") <- "disc_phase_type"
      return(value)
    }

  } else {
    stop('This matrix is not valid for either discrete or continuous phase-type distribution.')
  }
}


moment_ph <- function(obj, m) {
  e <- matrix(rep(1,nrow(obj$subint_mat)), nrow(obj$subint_mat), 1)
  inv <- solve(obj$subint_mat %^% m)
  moment <- as.numeric((-1) ** m * factorial(m) * obj$init_probs %*% inv %*% e)
  return(moment)
}


moment_individual <- function(element, obj){
  solve(-obj$subint_mat) %*% diag(obj$reward_mat[,element])
}

moment_row <- function(row, obj) {
  total <- diag(1, nrow(obj$subint_mat), ncol(obj$subint_mat))
  for (individual in lapply(row, moment_individual, obj = obj)) {
    total <- total %*% individual
  }
  sum(obj$init_probs %*% total)
}

perm <- function(v) {
  n <- length(v)
  if (n == 1) v
  else {
    X <- NULL
    for (i in 1:n) X <- rbind(X, cbind(v[i], perm(v[-i])))
    return(X)
  }
}

#' Moments of the multivariate phase-type distribution
#'
#' This function calculates the moments for the multivariate phase-type
#' distributions.
#'
#' The variables for which the moments are calculated can be specified in
#' the \code{v} vector as indices in the reward matrix of the
#' \code{mult_cont_phase_type} object.
#'
#' @param obj a mult_cont_phase_type.
#' @param v a vector or an integer.
#'
#' @usage moment_mph(obj, v)
#'
#' @export

moment_mph <- function(obj, v) {
  v <- matrix(v)
  X <- perm(v)
  sum(apply(X, 1, moment_row, obj = obj))
}

# require the partitions::parts() function (or write it again)
moment_mdph <- function(obj, v){
  Pv <- partitions::parts(v)
  udrj <- list()
  for (s in 1:v){
  }
}