---
title: "Benchmark Examples for UniCensorEM"
author: "Shikhar Tyagi"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Benchmark Examples}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
library(UniCensorEM)
```

## Introduction

This vignette demonstrates the use of UniCensorEM with several benchmark distributions. The package is designed to work with any user-supplied probability density function (PDF), cumulative distribution function (CDF), and survival function.

## Exponential Distribution

The exponential distribution is a simple one-parameter distribution often used for lifetime data.

### Complete Data

```{r exponential-complete}
# Define exponential distribution functions
pdf_exp <- function(x, theta) theta[1] * exp(-theta[1] * x)
cdf_exp <- function(x, theta) 1 - exp(-theta[1] * x)
survival_exp <- function(x, theta) exp(-theta[1] * x)

# Generate complete data
set.seed(123)
x_exp <- rexp(100, rate = 2)

# Fit using EM algorithm
fit_exp <- em_fit(
  data = x_exp,
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

# View results
print(fit_exp)
summary(fit_exp)
```

### Right Censored Data

```{r exponential-right-censor}
# Create right-censored data
set.seed(123)
x_exp_rc <- rexp(100, rate = 2)
status_exp_rc <- rep(1L, 100)
status_exp_rc[1:20] <- 0L  # Right-censor first 20 observations

# Fit with right censoring
fit_exp_rc <- em_fit(
  data = data.frame(x = x_exp_rc, status = status_exp_rc),
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "right_censor",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_exp_rc)
```

## Weibull Distribution

The Weibull distribution is a flexible two-parameter distribution commonly used in reliability analysis.

```{r weibull}
# Define Weibull distribution functions
pdf_weibull <- function(x, theta) {
  theta[1] * theta[2] * (x^(theta[2] - 1)) * exp(-theta[1] * x^theta[2])
}
cdf_weibull <- function(x, theta) 1 - exp(-theta[1] * x^theta[2])
survival_weibull <- function(x, theta) exp(-theta[1] * x^theta[2])

# Generate Weibull data
set.seed(123)
x_weibull <- rweibull(100, shape = 2, scale = 1)

# Fit using EM algorithm
fit_weibull <- em_fit(
  data = x_weibull,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  survival = survival_weibull,
  theta0 = c(1.0, 2.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_weibull)
```

## Gamma Distribution

The gamma distribution is a two-parameter family of continuous probability distributions.

```{r gamma}
# Define gamma distribution functions
pdf_gamma <- function(x, theta) {
  (theta[1]^theta[2] / gamma(theta[2])) * x^(theta[2] - 1) * exp(-theta[1] * x)
}
cdf_gamma <- function(x, theta) pgamma(x, shape = theta[2], rate = theta[1])
survival_gamma <- function(x, theta) pgamma(x, shape = theta[2], rate = theta[1], lower.tail = FALSE)

# Generate gamma data
set.seed(123)
x_gamma <- rgamma(100, shape = 2, rate = 1)

# Fit using EM algorithm
fit_gamma <- em_fit(
  data = x_gamma,
  pdf = pdf_gamma,
  cdf = cdf_gamma,
  survival = survival_gamma,
  theta0 = c(1.0, 2.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_gamma)
```

## Log-Normal Distribution

The log-normal distribution is useful for modeling variables that are positively skewed.

```{r lognormal}
# Define log-normal distribution functions
pdf_lnorm <- function(x, theta) {
  (1 / (x * theta[2] * sqrt(2 * pi))) * exp(-(log(x) - theta[1])^2 / (2 * theta[2]^2))
}
cdf_lnorm <- function(x, theta) plnorm(x, meanlog = theta[1], sdlog = theta[2])
survival_lnorm <- function(x, theta) plnorm(x, meanlog = theta[1], sdlog = theta[2], lower.tail = FALSE)

# Generate log-normal data
set.seed(123)
x_lnorm <- rlnorm(100, meanlog = 0, sdlog = 1)

# Fit using EM algorithm
fit_lnorm <- em_fit(
  data = x_lnorm,
  pdf = pdf_lnorm,
  cdf = cdf_lnorm,
  survival = survival_lnorm,
  theta0 = c(0.0, 1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_lnorm)
```

## Normal Distribution

The normal (Gaussian) distribution is the most common continuous probability distribution.

```{r normal}
# Define normal distribution functions
pdf_normal <- function(x, theta) {
  (1 / (theta[2] * sqrt(2 * pi))) * exp(-(x - theta[1])^2 / (2 * theta[2]^2))
}
cdf_normal <- function(x, theta) pnorm(x, mean = theta[1], sd = theta[2])
survival_normal <- function(x, theta) pnorm(x, mean = theta[1], sd = theta[2], lower.tail = FALSE)

# Generate normal data
set.seed(123)
x_normal <- rnorm(100, mean = 0, sd = 1)

# Fit using EM algorithm
fit_normal <- em_fit(
  data = x_normal,
  pdf = pdf_normal,
  cdf = cdf_normal,
  survival = survival_normal,
  theta0 = c(0.0, 1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_normal)
```

## Model Selection

Compare different models using information criteria.

```{r model-selection}
# Fit exponential to Weibull data
fit_exp_on_weibull <- em_fit(
  data = x_weibull,
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

# Compare information criteria
cat("Exponential fit on Weibull data:\n")
print(model_selection(fit_exp_on_weibull))

cat("\nWeibull fit on Weibull data:\n")
print(model_selection(fit_weibull))
```

## Goodness-of-Fit Statistics

Compute goodness-of-fit statistics for the fitted model.

```{r gof-stats}
gof <- gof_stats(fit_weibull)
print(gof)
```

## Bootstrap Standard Errors

Compute bootstrap standard errors for parameter estimates.

```{r bootstrap}
# Bootstrap with fewer iterations for speed
boot_results <- bootstrap_se(
  fit_exp,
  B = 100,
  type = "parametric",
  ci_type = "percentile",
  level = 0.95
)

print(boot_results)
```

## Diagnostic Plots

Create diagnostic plots for the EM algorithm fit.

```{r diagnostic-plots, fig.width=8, fig.height=8}
plot(fit_weibull, which = c(1, 2, 3), ask = FALSE)
```

## Residual Analysis

Perform residual analysis to check model fit.

```{r residual-analysis, fig.width=8, fig.height=8}
plot_residuals(fit_weibull, type = "cox-snell", which = c(1, 2, 3))
```

## EM Diagnostics

Examine EM algorithm convergence diagnostics.

```{r em-diagnostics}
diag <- em_diagnostics(fit_weibull)
print(diag)
```

## Interval Censoring Example

Example with interval-censored data.

```{r interval-censor}
# Create interval-censored data
set.seed(123)
x_interval <- rexp(100, rate = 2)
status_interval <- rep(1L, 100)
time2_interval <- rep(NA, 100)

# Make some observations interval-censored
interval_idx <- sample(1:100, 20)
status_interval[interval_idx] <- 3L
time2_interval[interval_idx] <- x_interval[interval_idx] + runif(20, 0.5, 1.5)

# Fit with interval censoring
fit_interval <- em_fit(
  data = data.frame(x = x_interval, status = status_interval, time2 = time2_interval),
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "interval_censor",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_interval)
```

## Conclusion

This vignette demonstrates the flexibility of UniCensorEM in handling various distributions and censoring schemes. The package can be extended to any user-defined distribution by providing the appropriate PDF, CDF, and survival functions.