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