R/cont_phase_type.R

In aumath-advancedr2019/ticphasetype: Phase-type distributions for bioinformatics

#' \code{cont_phase_type} class
#'
#' Description of the class \code{cont_phase_type}, which represents continuous phase-type
#' distributions.
#'
#' \code{cont_phase_type} is the generator function for the continuous phase-type distribution class
#' of the same name, which inherits from \code{list}. The class is generated by supplying a sub-intensity
#' matrix and an optional initial probability vector. If the initial probabilities are not specified,
#' then they are automatically generated as the first state having a probability of 1 and the rest
#' a probability of 0.
#'
#' @param subint_mat matrix
#' @param init_probs vector, a one-row matrix or \code{NULL}.
#'
#' @usage cont_phase_type(subint_mat = NULL, init_probs = NULL)
#'
#' @examples
#' subintensity_matrix = matrix(c(-1.5, 0, 0, 1.5, -1, 0, 0, 1, -0.5), ncol = 3)
#' cont_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)
#' cont_phase_type(subintensity_matrix, initial_probabilities)
#'
#' @export



cont_phase_type <- function(subint_mat = NULL, init_probs = NULL) {
  if (is.null(subint_mat)) {
    stop('Unable to construct the phase-type distribution. Please provide either the type or the subintensity matrix.')
  } else if (is.matrix(subint_mat)) {
    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.')
    } else if (sum(rowSums(subint_mat) > 0) != 0) {
      stop('The rowsums in subintensity matrix have to be non-positive')
    } else if ((is.vector(init_probs) & is.atomic(init_probs)) | is.matrix(init_probs)) {
      if (nrow(subint_mat) == length(init_probs)) {
        init_probs <- matrix(init_probs, nrow = 1)
      } else {
        stop('The length of the initial probabilities does not match the size of the subintensity matrix.')
      }
    } else {
      stop('The initial probabilities must be a a matrix with one row or a vector.')
    }
  } else {
    stop('The subintensity matrix must be a matrix.')
  }
  value <- list(subint_mat = subint_mat,
                init_probs = init_probs,
                defect = 1-sum(init_probs))
  attr(value, "class") <- "cont_phase_type"
  value
}


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



#' @export

mean.cont_phase_type <- function(x, ...) {
  moment_ph(x, 1)
}

#' Variance of cont_phase_type distributions
#'
#' Calculates the variance of continuous and discrete phase-type distributions,
#' represented by the \code{cont_phase_type} and the \code{disc_phase_type} classes
#' respectively
#'
#' @param obj a cont_phase_type or disc_phase_type object.
#'
#' @export

var <- function(obj) {
  UseMethod('var', obj)
}

#' @export

var.cont_phase_type <- function(obj) {
  moment_ph(obj, 2)-moment_ph(obj, 1)**2
}


#' @export

summary.cont_phase_type <- function(object, ...) {
  cat('\nSubintensity matrix:\n')
  print(object$subint_mat)
  cat('\nInitial probabilities:\n')
  print(object$init_probs)
  cat('\nDefect:\n')
  print(object$defect)
  cat('\nMean: ', mean(object), '\n', sep = '')
  cat('\nVariance: ', var(object), '\n\n', sep = '')
}