---
title: "Using Custom Distributions and Bootstrap Cross-Validation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using Custom Distributions and Bootstrap Cross-Validation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

This vignette demonstrates how to define a custom distribution (where user representation is supplied), fit it to data, compute generalized Process Capability Indices, and run a bootstrap cross-validation coverage diagnostic check.

## Setup

Load packages:

```{r setup}
library(gpci)
library(ggplot2)
```

## Defining a Custom Distribution

Suppose our process follows a skewed Weibull distribution. We define this distribution template using `define_distribution()`:

```{r define-custom}
# Define custom Weibull distribution template
custom_weibull <- define_distribution(
  name = "custom_weibull",
  cdf = function(x, shape, scale) pweibull(x, shape = shape, scale = scale),
  quantile = function(p, shape, scale) qweibull(p, shape = shape, scale = scale),
  params = list(shape = 2.0, scale = 10.0), # initial parameters
  support = c(0, Inf)
)
```

Let's verify that the PDF representation works by evaluating the PDF at $x = 5$:

```{r test-pdf}
# Theoretical PDF at x = 5 (using derived PDF)
do.call(custom_weibull$pdf, c(list(5), custom_weibull$params))

# Compare with the built-in dweibull:
dweibull(5, shape = 2, scale = 10)
```

They match exactly!

## Fitting the Custom Distribution to Skewed Data

We simulate process data from a Weibull distribution with shape = 2.5 and scale = 12:

```{r sim-skewed}
set.seed(42)
process_data <- rweibull(80, shape = 2.5, scale = 12.0)
```

We fit our custom Weibull template to this dataset using Maximum Likelihood Estimation:

```{r fit-custom}
fitted_weibull <- fit_distribution(
  data = process_data,
  dist = custom_weibull,
  method = "mle"
)

# Print fitted parameters
print(fitted_weibull$params)
```

## Robust Capability Analysis

For skewed processes, classical moment-based indices like $C_p$ can be misleading. We instead run capability analysis using the robust **quantile-based mode** (Clements / Pearn-Chen approach), which replaces the process mean with the median and the $6\sigma$ spread with $Q(0.99865) - Q(0.00135)$. We specify:
* LSL = 3.0
* USL = 20.0
* Target = 11.0

We compute robust indices: `Cp_q` (Clements Cp), `Cpk_q` (Clements Cpk), `CNpk` (Pearn-Chen symmetric Cpk), and the specialized literature indices: `CpTk` (Maiti et al., 2010), `Spmk` (Dey & Saha, 2019), and `CNpmc` (Alotaibi et al., 2022) with a tolerance cost function:

```{r capability-robust}
fit_robust <- capability(
  data = process_data,
  distribution = fitted_weibull,
  USL = 20,
  LSL = 3,
  target = 11,
  indices = c("Cp_q", "Cpk_q", "CNpk", "CpTk", "Spmk", "CNpmc"),
  mode = "quantile",
  fit = FALSE # Already fitted
)

print(fit_robust)
```

We can plot the process density and specifications:

```{r plot-robust-density}
plot(fit_robust, type = "density")
```

## Bootstrap Cross-Validation Diagnostics

To check how trustworthy the confidence intervals are for our custom distribution and sample size, we run the **bootstrap cross-validation** diagnostic module. 

Treating the point estimate on the original data as the reference value, this module generates $B_2$ synthetic samples, computes capability and bootstrap CIs on each, and tracks how often the CI covers the reference value.

```{r run-cv, eval = FALSE}
# Evaluate coverage calibration for percentile CIs
cv <- boot_cv(
  fit = fit_robust,
  B2 = 50,           # Number of synthetic samples
  B = 200,          # Bootstrap reps per synthetic sample
  alpha = c(0.10, 0.05),
  method = "percentile",
  type = "parametric",
  parallel = FALSE
)

# Print validation metrics (bias, RMSE, empirical coverage vs. nominal)
print(cv)
```