---
title: "Introduction to BivLaplaceRL"
author: "Mahesh Divakaran"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to BivLaplaceRL}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

**BivLaplaceRL** is an R package that implements statistical methods from three
research papers on bivariate Laplace transforms and entropy measures in
reliability theory:

1. Bivariate Laplace transform of residual lives (Jayalekshmi, Rajesh & Nair, 2022)
2. Bivariate Laplace transform order of reversed residual lives (Jayalekshmi & Rajesh)
3. Residual entropy generating function (Smitha, Rajesh & Jayalekshmi, 2024)

```{r setup}
library(BivLaplaceRL)
```

## Parametric Distributions

The package provides four bivariate distributions commonly used in reliability
modelling.

### Gumbel Bivariate Exponential

```{r gumbel}
# Simulate 200 observations
set.seed(42)
dat <- rgumbel_biv(200, k1 = 1, k2 = 1.5, theta = 0.3)
head(dat)

# Survival function at a point
sgumbel_biv(1, 1, k1 = 1, k2 = 1.5, theta = 0.3)
```

### Farlie-Gumbel-Morgenstern (FGM)

```{r fgm}
set.seed(1)
dat_fgm <- rfgm_biv(100, theta = 0.5)
head(dat_fgm)
pfgm_biv(0.5, 0.6, theta = 0.5)
```

### Schur-Constant

```{r schur}
set.seed(2)
dat_sc <- rschur_biv(100, lambda = 1)
# Marginals are standard exponential
summary(rowSums(dat_sc))  # T = X1 + X2 ~ Exp(lambda)
```

## Bivariate Laplace Transform of Residual Lives

### Closed-Form (Gumbel Distribution)

The key result from Jayalekshmi et al. (2022) is:

$$L^*_{X_{t_1|t_2}}(s_1) = \frac{1}{k_1 + s_1 + \theta t_2}, \quad
  L^*_{X_{t_2|t_1}}(s_2) = \frac{1}{k_2 + s_2 + \theta t_1}.$$

```{r closed_form}
blt_residual_gumbel(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5,
                    k1 = 1, k2 = 1, theta = 0.5)
```

### Numerical Computation

For a user-supplied survival function:

```{r 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)
```

### Nonparametric Estimation

```{r 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)
```

### Simulation Study

```{r 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)
```

## Bivariate Hazard Gradient and MRL

```{r hazard_mrl}
# Hazard gradient
biv_hazard_gradient(t1 = 1, t2 = 1, k1 = 1, k2 = 1, theta = 0.3)

# Mean residual life
biv_mean_residual(t1 = 0.5, t2 = 0.5, k1 = 1, k2 = 1, theta = 0)
```

## NBUHR / NWUHR Aging Classification

```{r 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")
```

## Bivariate Laplace Transform of Reversed Residual Lives

```{r reversed}
# FGM distribution
blt_reversed(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5, theta = 0.3)

# G form (used for characterisation)
blt_reversed(s1 = 1, s2 = 1, t1 = 0.5, t2 = 0.5, theta = 0.3, g_form = TRUE)

# Reversed hazard gradient
biv_rhazard_gradient(x1 = 0.5, x2 = 0.5, theta = 0.3)

# Reversed MRL
biv_rmrl(x1 = 0.5, x2 = 0.5, theta = 0.3)
```

## Stochastic Orders

### BLt-rl Order

```{r 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")
```

### Weak Bivariate Hazard Rate Order

```{r whr_order}
res <- biv_whr_order(sX, sY, t1_grid = c(0.5, 1), t2_grid = c(0.5, 1))
cat("X <=_whr Y:", res$order_holds, "\n")
```

## Residual Entropy Generating Function

### Shannon Entropy and Golomb IGF

```{r entropy_basics}
# Shannon entropy of Exp(1) = 1
shannon_entropy(function(x) dexp(x, 1))

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

### Dynamic REGF

```{r regf}
f  <- function(x) dexp(x, 1)
Fb <- function(x) pexp(x, 1, lower.tail = FALSE)

# REGF at various t values
t_vals <- seq(0, 2, 0.5)
sapply(t_vals, function(t) residual_info_gen(f, Fb, alpha = 2, t = t))
```

### Nonparametric Estimation

```{r np_regf}
set.seed(7)
x <- rexp(300, rate = 1)
np_residual_info_gen(x, alpha = 2, t = 0.5)

# True value
residual_info_gen(f, Fb, alpha = 2, t = 0.5)
```

### REGF Profile and Characterisation

```{r regf_profile, fig.cap="REGF profile across alpha"}
prof <- regf_profile(f, Fb, t = 0.5, plot = TRUE)
head(prof)
```

```{r regf_charac}
# Log-linearity in t characterises exponential distribution
ch <- regf_characterise(f, Fb, alpha = 2, t_grid = seq(0, 2, 0.4))
ch
```

### Simulation Study for REGF Estimator

```{r sim_regf}
sim_regf(lambda = 1, alpha = 2, t = 0, n_obs = 200, n_sim = 50)
```

## References

Jayalekshmi S., Rajesh G., Nair N.U. (2022). Bivariate Laplace transform of
residual lives and their properties. *Communications in Statistics—Theory and
Methods*. https://doi.org/10.1080/03610926.2022.2085874

Jayalekshmi S., Rajesh G. Bivariate Laplace transform order and ordering of
reversed residual lives. *International Journal of Reliability, Quality and
Safety Engineering*.

Smitha S., Rajesh G., Jayalekshmi S. (2024). On residual entropy generating
function. *Journal of the Indian Statistical Association*, 62(1), 81–93.
