Using Custom Distributions and Bootstrap Cross-Validation

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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:

library(gpci)
library(ggplot2)

Defining a Custom Distribution

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

# 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\):

# Theoretical PDF at x = 5 (using derived PDF)
do.call(custom_weibull$pdf, c(list(5), custom_weibull$params))
#> [1] 0.07788008

# Compare with the built-in dweibull:
dweibull(5, shape = 2, scale = 10)
#> [1] 0.07788008

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:

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:

fitted_weibull <- fit_distribution(
  data = process_data,
  dist = custom_weibull,
  method = "mle"
)

# Print fitted parameters
print(fitted_weibull$params)
#> $shape
#> [1] 2.100232
#> 
#> $scale
#> [1] 11.61896

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:

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)
#> --- Process Capability Analysis (Class: gpcifit) ---
#> Distribution:  custom_weibull 
#> Parameters:    shape = 2.1002, scale = 11.619 
#> Spec Limits:  LSL = 3 , USL = 20 , Target = 11 
#> Mode:          quantile 
#> Expected Nonconforming (p_hat):  10.0318 %
#> 
#> Point Estimates of Capability Indices:
#>   Cp_q  Cpk_q   CNpk   CpTk   Spmk  CNpmc 
#> 0.6060 0.5450 0.4819 0.8029 0.5426 0.5736

We can plot the process density and specifications:

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.

# 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)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.