# R/4c.biviph.R In matrixdist: Statistics for Matrix Distributions

#### Documented in biviph

#' Bivariate inhomogeneous phase-type distributions
#'
#' Class of objects for bivariate inhomogeneous phase-type distributions.
#'
#' @slot name Name of the phase type distribution.
#' @slot gfun A list comprising of the parameters.
#'
#' @return Class object.
#' @export
#'
setClass("biviph",
contains = c("bivph"),
slots = list(
gfun = "list"
)
)

#' Constructor function for bivariate inhomogeneous phase-type distributions
#'
#' @param bivph An object of class \linkS4class{bivph}.
#' @param alpha A probability vector.
#' @param S11 A sub-intensity matrix.
#' @param S12 A matrix.
#' @param S22 A sub-intensity matrix.
#' @param dimensions The dimensions of the bivariate phase-type (if no parameters are provided).
#' @param gfun Vector of inhomogeneity transforms.
#' @param gfun_pars List of parameters for the inhomogeneity functions.
#'
#' @return An object of class \linkS4class{biviph}.
#' @export
#'
#' @examples
#' under_bivph <- bivph(dimensions = c(3, 3))
#' biviph(under_bivph, gfun = c("weibull", "pareto"), gfun_pars = list(c(2), c(3)))
biviph <- function(bivph = NULL,
gfun = NULL,
gfun_pars = NULL,
alpha = NULL,
S11 = NULL,
S12 = NULL,
S22 = NULL,
dimensions = c(3, 3)) {
if (all(is.null(c(gfun, gfun_pars)))) {
stop("input inhomogeneity function and parameters")
}
d <- length(gfun)

if (is.null(bivph)) {
bivph <- bivph(alpha = alpha, S11 = S11, S12 = S12, S22 = S22, structure = structure, dimensions = dimensions)
}

if (!all(gfun %in% c("pareto", "weibull", "lognormal", "loglogistic", "gompertz", "gev", "identity"))) {
stop("invalid gfun for at least one marginal")
}

if (all(gfun %in% c("pareto", "weibull", "lognormal", "gompertz"))) {
for (i in 1:2) {
if (is.null(gfun_pars[[i]])) gfun_pars[[i]] <- 1
if (length(gfun_pars[[i]]) != 1 | sum(gfun_pars[[i]] <= 0) > 0) {
stop(paste("gfun parameter for marginal", i, "should be positive and of length one"))
} else {
names(gfun_pars[[i]]) <- paste("beta", i, sep = "")
}
}
}

if (any(gfun %in% c("gev"))) {
index_gev <- which(gfun == "gev")
for (i in index_gev) {
index_gev2 <- which(gfun[[i]] == "gev")
if (is.null(gfun_pars[[i]])) gfun_pars[[i]] <- c(0, 1, 1)
if (length(gfun_pars[[i]]) != 3 | (gfun_pars[[i]][2] > 0) == FALSE) {
stop("gfun parameter should be of length three: mu, sigma, xi, and sigma > 0")
} else {
names(gfun_pars) <- c("mu", "sigma", "xi")
}
}
}

if (any(gfun == "loglogistic")) {
index_log <- which(gfun == "loglogistic")
for (i in index_log) {
if (is.null(gfun_pars[[i]])) gfun_pars[[i]] <- c(1, 1)
if (length(gfun_pars[[i]]) != 2 | (gfun_pars[[i]][1] <= 0) | (gfun_pars[[i]][2] <= 0)) {
stop("gfun parameter should be positive and of length two: alpha, theta > 0")
} else {
names(gfun_pars[[i]]) <- c(paste("alpha", i, sep = ""), paste("theta", i, sep = ""))
}
}
}

ginv <- list()
ginv_prime <- list()
lambda <- list()
lambda_prime <- list()

for (i in 1:2) {
f1 <- function(beta, t) t^(beta)
f2 <- function(beta, t) log(t / beta + 1)
f3 <- function(beta, t) log(t + 1)^(beta)
f4 <- function(beta, t) log((t / beta[1])^(beta[2]) + 1)
f5 <- function(beta, t) (exp(t * beta) - 1) / beta
f6 <- function(beta, t, w) revers_data_trans(t, w, beta)
nb <- which(gfun[i] == c("weibull", "pareto", "lognormal", "loglogistic", "gompertz", "gev"))
ginv[[i]] <- base::eval(parse(text = paste("f", nb, sep = "")))

f1 <- function(beta, t) t^(beta) * log(t)
f2 <- function(beta, t) -t / (beta * t + beta^2)
f3 <- function(beta, t) log(t + 1)^(beta) * log(log(t + 1))
f4 <- NA
f5 <- function(beta, t) exp(t * beta) * (t * beta - 1) / beta^2
f6 <- NA
nb <- which(gfun[i] == c("weibull", "pareto", "lognormal", "loglogistic", "gompertz", "gev"))
ginv_prime[[i]] <- base::eval(parse(text = paste("f", nb, sep = "")))

f1 <- function(beta, t) beta * t^(beta - 1)
f2 <- function(beta, t) (t + beta)^(-1)
f3 <- function(beta, t) beta * log(t + 1)^(beta - 1) / (t + 1)
f4 <- NA
f5 <- function(beta, t) exp(t * beta)
f6 <- NA
nb <- which(gfun[i] == c("weibull", "pareto", "lognormal", "loglogistic", "gompertz", "gev"))
lambda[[i]] <- base::eval(parse(text = paste("f", nb, sep = "")))

f1 <- function(beta, t) t^(beta - 1) + beta * t^(beta - 1) * log(t)
f2 <- function(beta, t) -(t + beta)^(-2)
f3 <- function(beta, t) log(t + 1)^(beta - 1) / (t + 1) + beta * log(t + 1)^(beta - 1) * log(log(t + 1)) / (t + 1)
f4 <- NA
f5 <- function(beta, t) t * exp(t * beta)
f6 <- NA
nb <- which(gfun[i] == c("weibull", "pareto", "lognormal", "loglogistic", "gompertz", "gev"))
lambda_prime[[i]] <- base::eval(parse(text = paste("f", nb, sep = "")))
}
name <- if (is(bivph, "biviph")) bivph@name else paste("inhomogeneous ", bivph@name, sep = "")

methods::new("biviph",
name = name,
pars = bivph@pars,
gfun = list(
name = gfun,
pars = gfun_pars,
inverse = ginv,
inverse_prime = ginv_prime,
intensity = lambda,
intensity_prime = lambda_prime
)
)
}

#' Show method for bivariate inhomogeneous phase-type distributions
#'
#' @param object An object of class \linkS4class{biviph}.
#' @importFrom methods show
#' @export
#'
setMethod("show", "biviph", function(object) {
cat("object class: ", methods::is(object)[[1]], "\n", sep = "")
cat("name: ", object@name, "\n", sep = "")
cat("parameters: ", "\n", sep = "")
print(object@pars)
cat("g-function name:", object@gfun\$name, "\n")
cat("parameters: ", "\n", sep = "")
methods::show(object@gfun\$pars)
})

#' Simulation method for bivariate inhomogeneous phase-type distributions
#'
#' @param x An object of class \linkS4class{biviph}.
#' @param n An integer of length of realization.
#'
#' @return A realization of independent and identically distributed bivariate
#'  inhomogeneous phase-type vector.
#' @export
#'
#' @examples
#' under_bivph <- bivph(dimensions = c(3, 3))
#' obj <- biviph(under_bivph, gfun = c("weibull", "pareto"), gfun_pars = list(c(2), c(3)))
#' sim(obj, n = 100)
setMethod("sim", c(x = "biviph"), function(x, n = 1000) {
p1_aux <- dim(x@pars\$S11)[1]
p2_aux <- dim(x@pars\$S22)[1]
alpha_aux <- c(x@pars\$alpha, rep(0, p2_aux))
S_aux <- merge_matrices(x@pars\$S11, x@pars\$S12, x@pars\$S22)
R_aux <- matrix(c(c(rep(1, p1_aux), rep(0, p2_aux)), c(rep(0, p1_aux), rep(1, p2_aux))), ncol = 2)
rMIPHstar(n, alpha_aux, S_aux, R_aux, x@gfun\$name, x@gfun\$pars)
})

#' Marginal method for biviph class
#'
#' @param x An object of class \linkS4class{biviph}.
#' @param mar Indicator of which marginal.
#' @return An object of the of class \linkS4class{iph}.
#' @export
#'
#' @examples
#' under_bivph <- bivph(dimensions = c(3, 3))
#' obj <- biviph(under_bivph, gfun = c("weibull", "pareto"), gfun_pars = list(c(2), c(3)))
#' marginal(obj, 1)
setMethod("marginal", c(x = "biviph"), function(x, mar = 1) {
if (!(mar %in% 1:2)) {
stop("maringal provided not available")
}
if (mar == 1) {
x0 <- iph(ph(alpha = x@pars\$alpha, S = x@pars\$S11), gfun = x@gfun\$name[mar], gfun_pars = x@gfun\$pars[[mar]])
} else {
alpha0 <- x@pars\$alpha %*% base::solve(-x@pars\$S11) %*% x@pars\$S12
x0 <- iph(ph(alpha = alpha0, S = x@pars\$S22), gfun = x@gfun\$name[mar], gfun_pars = x@gfun\$pars[[mar]])
}
x0
})

#' Density method for bivariate inhomogeneous phase-type distributions
#'
#' @param x An object of class \linkS4class{biviph}.
#' @param y A matrix of locations.
#'
#' @return A vector containing the joint density evaluations at the given locations.
#' @export
#'
#' @examples
#' under_bivph <- bivph(dimensions = c(3, 3))
#' obj <- biviph(under_bivph, gfun = c("weibull", "pareto"), gfun_pars = list(c(2), c(3)))
#' dens(obj, matrix(c(0.5, 1), ncol = 2))
setMethod("dens", c(x = "biviph"), function(x, y) {
if (is.vector(y)) {
y <- t(y)
}
gfun.pars <- x@gfun\$pars
y1_inv <- x@gfun\$inverse[[1]](gfun.pars[[1]], y[, 1])
y1_int <- x@gfun\$intensity[[1]](gfun.pars[[1]], y[, 1])
y2_inv <- x@gfun\$inverse[[2]](gfun.pars[[2]], y[, 2])
y2_int <- x@gfun\$intensity[[2]](gfun.pars[[2]], y[, 2])
dens <- bivph_density(cbind(y1_inv, y2_inv), x@pars\$alpha, x@pars\$S11, x@pars\$S12, x@pars\$S22) * y1_int * y2_int
unname(dens)
})

#' Coef method for biviph class
#'
#' @param object An object of class \linkS4class{biviph}.
#'
#' @return Parameters of bivariate inhomogeneous phase-type model.
#' @export
#'
#' @examples
#' under_bivph <- bivph(dimensions = c(3, 3))
#' obj <- biviph(under_bivph, gfun = c("weibull", "pareto"), gfun_pars = list(c(2), c(3)))
#' coef(obj)
setMethod("coef", c(object = "biviph"), function(object) {
L <- object@pars
L\$gfun <- object@gfun\$name
L\$gpars <- object@gfun\$pars
L
})

## Try the matrixdist package in your browser

Any scripts or data that you put into this service are public.

matrixdist documentation built on Aug. 8, 2023, 5:06 p.m.