Getting Started with ProcessCapabilityR

Shikhar Tyagi, Sumit Kumar, Vrijesh Tripathi

2026-07-23

Introduction

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.


1. Normal-Case Sanity Check

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   = 3

For 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 = 1

2. Off-Center Process

Now 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.333333

Note that \(C_p\) is unchanged (it ignores centering), while \(C_{pk}\) drops to 0.667 because the process is closer to the USL.


3. Taguchi Indices — Off-Target Example

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


4. Performance Indices

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

With \(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.


5. Generalized \(C_{py}\): Normal Distribution

\(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) = 1

With 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.002644

6. Non-Normal Example: Weibull Distribution

d_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.047193

7. Non-Normal Example: Gamma Distribution

d_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.043152

8. Bootstrap Confidence Intervals

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

9. Sensitivity Sweep: \(\sigma\)

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)
plot(grid_sigma, x_axis = "sigma")

10. Sensitivity Sweep: \(p_0\) for \(C_{py}\)

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)
plot(grid_cpy, x_axis = "p0")

References