## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

## ----setup--------------------------------------------------------------------
library(BivLaplaceRL)

## ----gumbel-------------------------------------------------------------------
set.seed(42)
dat <- rgumbel_biv(200, k1 = 1, k2 = 1.5, theta = 0.3)
head(dat)
sgumbel_biv(1, 1, k1 = 1, k2 = 1.5, theta = 0.3)

## ----fgm----------------------------------------------------------------------
set.seed(1)
dat_fgm <- rfgm_biv(100, theta = 0.5)
pfgm_biv(0.5, 0.6, theta = 0.5)

## ----schur--------------------------------------------------------------------
set.seed(2)
dat_sc <- rschur_biv(100, lambda = 1)
summary(rowSums(dat_sc))

## ----closed_form--------------------------------------------------------------
blt_residual_gumbel(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5,
                    k1 = 1, k2 = 1, theta = 0.5)

## ----numerical----------------------------------------------------------------
my_surv <- function(x1, x2) exp(-x1 - x2 - 0.3 * x1 * x2)
blt_residual(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5,
             surv_fn = my_surv, star = TRUE)

## ----np_estim-----------------------------------------------------------------
set.seed(123)
dat <- rgumbel_biv(500, k1 = 1, k2 = 1, theta = 0.5)
np_blt_residual(dat, s1 = 1, s2 = 1, t1 = 0.3, t2 = 0.3)

# True value
blt_residual_gumbel(s1 = 1, s2 = 1, t1 = 0.3, t2 = 0.3, theta = 0.5)

## ----sim_study----------------------------------------------------------------
sim_blt_residual(n_obs = 200, n_sim = 50, s1 = 1, s2 = 1,
                 t1 = 0.3, t2 = 0.3, k1 = 1, k2 = 1, theta = 0.5)

## ----hazard_mrl---------------------------------------------------------------
biv_hazard_gradient(t1 = 1, t2 = 1, k1 = 1, k2 = 1, theta = 0.3)
biv_mean_residual(t1 = 0.5, t2 = 0.5, k1 = 1, k2 = 1, theta = 0)

## ----aging--------------------------------------------------------------------
res <- nbuhr_test(t2 = 0.5, k1 = 1, k2 = 1, theta = 0.3,
                  t1_grid = seq(0.1, 3, 0.3))
cat("Component 1:", res$class1, "\n")
cat("Component 2:", res$class2, "\n")

## ----reversed-----------------------------------------------------------------
blt_reversed(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5, theta = 0.3)
biv_rhazard_gradient(x1 = 0.5, x2 = 0.5, theta = 0.3)
biv_rmrl(x1 = 0.5, x2 = 0.5, theta = 0.3)

## ----blt_order----------------------------------------------------------------
sX <- function(x1, x2) sgumbel_biv(x1, x2, k1 = 1, k2 = 1, theta = 0.2)
sY <- function(x1, x2) sgumbel_biv(x1, x2, k1 = 2, k2 = 2, theta = 0.2)
res <- blt_order_residual(sX, sY, s1 = 1, s2 = 1,
                          t1_grid = c(0.5, 1), t2_grid = c(0.5, 1))
cat("X <=_BLt-rl Y:", res$order_holds, "\n")

res2 <- biv_whr_order(sX, sY, t1_grid = c(0.5, 1), t2_grid = c(0.5, 1))
cat("X <=_whr Y:", res2$order_holds, "\n")

## ----uni_lt-------------------------------------------------------------------
f  <- function(x) dexp(x, 1)
Fb <- function(x) pexp(x, 1, lower.tail = FALSE)

# At t=0: L_X(1, 0) = 1/(1+1) = 0.5
lt_residual(f, Fb, s = 1, t = 0)

# At t=0.5
lt_residual(f, Fb, s = 1, t = 0.5)

## ----uni_hr_mrl---------------------------------------------------------------
# Exp(1): constant hazard rate = 1
hazard_rate(f, Fb, t = c(0.5, 1, 2))

# Exp(1): constant MRL = 1 (memoryless)
mean_residual(Fb, t = c(0, 0.5, 1, 2))

# Gamma(2,1): increasing hazard, decreasing MRL
fG  <- function(x) dgamma(x, shape = 2, rate = 1)
FbG <- function(x) pgamma(x, shape = 2, rate = 1, lower.tail = FALSE)
hazard_rate(fG, FbG, t = c(0.5, 1, 2))
mean_residual(FbG, t = c(0, 0.5, 1))

## ----uni_np-------------------------------------------------------------------
set.seed(7)
x <- rexp(500, rate = 1)

# Nonparametric estimate at s=1, t=0: true value = 0.5
np_lt_residual(x, s = 1, t = 0)

# At t=0.5
np_lt_residual(x, s = 1, t = 0.5)

## ----uni_orders---------------------------------------------------------------
f1  <- function(x) dexp(x, 1)
Fb1 <- function(x) pexp(x, 1, lower.tail = FALSE)
f2  <- function(x) dexp(x, 2)
Fb2 <- function(x) pexp(x, 2, lower.tail = FALSE)

# LT-rl order: Exp(1) <=_Lt-rl Exp(2)?
lt_rl_order(f1, Fb1, f2, Fb2,
            s_grid = c(0.5, 1, 2), t_grid = c(0, 0.5, 1))$order_holds

# Hazard rate order: Exp(1) <=_hr Exp(2)?
hr_order(f1, Fb1, f2, Fb2, t_grid = c(0.5, 1, 2))$order_holds

# MRL order: Exp(2) <=_mrl Exp(1)?
mrl_order(Fb2, Fb1, t_grid = c(0, 0.5, 1, 2))$order_holds

## ----entropy------------------------------------------------------------------
# Shannon entropy of Exp(1) = 1
shannon_entropy(function(x) dexp(x, 1))

# Golomb IGF of Exp(1) at alpha=2: 0.5
info_gen_function(function(x) dexp(x, 1), alpha = 2)