Process Capability Indices (PCIs) are fundamental tools in Statistical Quality Control. They quantify the ability of a manufacturing or service process to produce output within engineering specification limits. Classical indices such as \(C_p\) and \(C_{pk}\) assume the quality characteristic follows a normal distribution.
ProcessCapabilityR extends this framework by implementing the Generalized Process Capability Index \(C_{py}\) (Maiti, Saha & Nanda, 2010), which works for any continuous or discrete distribution the user plugs in. The classical indices are recovered as special cases under normality.
This vignette demonstrates the full API through worked examples.
Consider a centered process with \(USL = 63\), \(LSL = 57\), \(\mu = 60\), \(\sigma = 1\).
# All classical capability indices
cat("Cp =", cp(LSL = 57, USL = 63, sigma = 1), "\n")
#> Cp = 1
cat("Cpk =", cpk(LSL = 57, USL = 63, mu = 60, sigma = 1), "\n")
#> Cpk = 1
cat("Cpu =", cpu(USL = 63, mu = 60, sigma = 1), "\n")
#> Cpu = 1
cat("Cpl =", cpl(LSL = 57, mu = 60, sigma = 1), "\n")
#> Cpl = 1
cat("Z =", z_level(LSL = 57, USL = 63, mu = 60, sigma = 1), "\n")
#> Z = 3For a perfectly centered process the specification width equals \(6\sigma\), so \(C_p = C_{pk} = 1.0\) and \(Z = 3\).
# Taguchi indices with target = process mean
cat("Cpm =", cpm(LSL = 57, USL = 63, mu = 60, sigma = 1, target = 60), "\n")
#> Cpm = 1
cat("Cpmk =", cpmk(LSL = 57, USL = 63, mu = 60, sigma = 1, target = 60), "\n")
#> Cpmk = 1Now shift the mean to \(\mu = 61\) (off-center but within limits):
cat("Cp =", cp(LSL = 57, USL = 63, sigma = 1), "\n")
#> Cp = 1
cat("Cpk =", cpk(LSL = 57, USL = 63, mu = 61, sigma = 1), "\n")
#> Cpk = 0.6666667
cat("Cpu =", cpu(USL = 63, mu = 61, sigma = 1), "\n")
#> Cpu = 0.6666667
cat("Cpl =", cpl(LSL = 57, mu = 61, sigma = 1), "\n")
#> Cpl = 1.333333Note that \(C_p\) is unchanged (it ignores centering), while \(C_{pk}\) drops to 0.667 because the process is closer to the USL.
Set \(\mu = 60\), \(\sigma = 1\), and a target \(T = 59\) (target is different from the process mean):
cat("Cpm =", cpm(LSL = 57, USL = 63, mu = 60, sigma = 1, target = 59), "\n")
#> Cpm = 0.7071068
cat("Cpmk =", cpmk(LSL = 57, USL = 63, mu = 60, sigma = 1, target = 59), "\n")
#> Cpmk = 0.7071068\(C_{pm} = 1/\sqrt{2} \approx 0.707\) — the Taguchi denominator \(\sqrt{\sigma^2 + (\mu - T)^2}\) inflates because \(\mu \neq T\).
Using \(\bar{x} = 60\), \(s = 1.5\) (long-term SD, larger than \(\sigma\)):
cat("Pp =", pp(LSL = 57, USL = 63, s = 1.5), "\n")
#> Pp = 0.6666667
cat("Ppk =", ppk(LSL = 57, USL = 63, xbar = 60, s = 1.5), "\n")
#> Ppk = 0.6666667
cat("Ppu =", ppu(USL = 63, xbar = 60, s = 1.5), "\n")
#> Ppu = 0.6666667
cat("Ppl =", ppl(LSL = 57, xbar = 60, s = 1.5), "\n")
#> Ppl = 0.6666667With \(s > \sigma\), the performance indices (\(P_p = 0.667\)) are lower than the capability indices (\(C_p = 1.0\)), indicating extra variation from long-term sources.
\(C_{py}\) is defined as the ratio of actual to desired yield: \[C_{py} = \frac{F(USL) - F(LSL)}{F(UDL) - F(LDL)} = \frac{p}{p_0}\]
Under normality with \(LDL = \mu - 3\sigma\) and \(UDL = \mu + 3\sigma\), the actual and desired yields are identical, so \(C_{py} = 1\):
d_norm <- pci_dist_normal(mean = 60, sd = 1)
cat("Cpy (spec = desirable) =",
cpy(d_norm, LSL = 57, USL = 63, LDL = 57, UDL = 63), "\n")
#> Cpy (spec = desirable) = 1With wider spec limits (\(LSL = 56\), \(USL = 64\)), \(C_{py} > 1\):
cat("Cpy (wider spec) =",
cpy(d_norm, LSL = 56, USL = 64, LDL = 57, UDL = 63), "\n")
#> Cpy (wider spec) = 1.002644d_weibull <- pci_dist(
pdf = function(x, shape, scale) dweibull(x, shape = shape, scale = scale),
cdf = function(x, shape, scale) pweibull(x, shape = shape, scale = scale),
params = list(shape = 2, scale = 10),
support = c(0, 50)
)
cat("Cpy (Weibull, p0=0.95) =",
cpy(d_weibull, LSL = 2, USL = 20, p0 = 0.95), "\n")
#> Cpy (Weibull, p0=0.95) = 0.9920777
cat("Cpy (Weibull, p0=0.90) =",
cpy(d_weibull, LSL = 2, USL = 20, p0 = 0.90), "\n")
#> Cpy (Weibull, p0=0.90) = 1.047193d_gamma <- pci_dist(
pdf = function(x, shape, rate) dgamma(x, shape = shape, rate = rate),
cdf = function(x, shape, rate) pgamma(x, shape = shape, rate = rate),
params = list(shape = 5, rate = 0.5),
support = c(0, 60)
)
cat("Cpy (Gamma, p0=0.95) =",
cpy(d_gamma, LSL = 2, USL = 25, p0 = 0.95), "\n")
#> Cpy (Gamma, p0=0.95) = 1.043152set.seed(42)
d <- pci_dist_normal(mean = 60, sd = 1)
ci_pct <- pci_ci("Cp", dist = d, n = 30, LSL = 57, USL = 63,
alpha = 0.05, B = 500, method = "percentile")
print(ci_pct)
#> Index : Cp
#> Estimate : 1.0024
#> 95% CI [percentile]: (0.7900, 1.3692)set.seed(42)
ci_bca <- pci_ci("Cp", dist = d, n = 30, LSL = 57, USL = 63,
alpha = 0.05, B = 500, method = "bca")
print(ci_bca)
#> Index : Cp
#> Estimate : 1.0024
#> 95% CI [bca]: (0.7862, 1.3657)Sweep \(\sigma\) from 0.5 to 2.0 with confidence bands at 90%, 95%, 97%, and 99%:
d <- pci_dist_normal(mean = 60, sd = 1)
grid_sigma <- pci_grid("Cp",
dist = d,
LSL = 57, USL = 63,
sigma_vals = seq(0.5, 2.0, by = 0.1),
mu = 60,
alpha_vals = c(0.10, 0.05, 0.03, 0.01),
n = 30, B = 500)
head(grid_sigma)d <- pci_dist_normal(mean = 60, sd = 1)
grid_cpy <- pci_grid("Cpy",
dist = d,
LSL = 57, USL = 63,
p0_vals = c(0.90, 0.95, 0.99),
alpha_vals = c(0.10, 0.05, 0.03, 0.01),
n = 30, B = 500)
head(grid_cpy)