Numerical validation against the methodological literature

Validation protocol

Unit tests measure whether code satisfies its contracts; scientific validation answers a different question: whether the implemented statistics reproduce the results that motivated them. IQCC’s validation records, for every numerical fit, the provenance of the published value, the parametrization used, the calculated value, and the ratio of the absolute discrepancy to an explicit tolerance.

Discrepancies are judged against a tolerance of half a unit in the last published decimal place by default. A small number of cells whose published value rounds from an independent calculation to within one unit of the last decimal instead receive an explicitly documented one-unit exception (for example the p-chart, DS-np, R-chart and generalized-variance borderline cells listed in the corresponding tests). Every exception is recorded per cell in the provenance, and no tolerance is widened beyond one unit of the last printed decimal. Tolerances are propagated through transformations that scale them linearly. Two failure modes are kept distinct: numerical error (which is controlled and machine-checkable) and Monte Carlo error (which is statistical and shrinks only with more draws). When a publication’s parametrization is ambiguous and cannot be settled from the source text, the package stops rather than silently choosing a convention.

Every result below is classified with one of these evidence labels:

Label	Meaning
`published_table`	Value transcribed from a published table, figure, or equation
`independent_derivation`	Result recomputed through a separate mathematical route
`exact_discrete`	Exact sum or CDF of the discrete distribution
`numerical_exact_distribution`	Exact mathematical distribution evaluated numerically with controlled error
`monte_carlo`	Simulated estimate with explicit seed and uncertainty
`property_test`	Monotonicity, invariance, API equivalence, or another property
`partial_reproduction`	Source insufficient for a complete reconstruction

The vignette never runs testthat; it calls the public API and builds its own evidence tables from the same fixtures used by the test suite.

Corrected p chart

Source and objective. Joekes, S. and Barbosa, E. P. (2013). An improved attribute control chart for monitoring non-conforming proportion in high quality processes. Control Engineering Practice, 21, 407-412. doi:10.1016/j.conengprac.2012.12.005.

Parametrization. Binomial data with in-control proportion p, subgroup size n, nominal two-sided risk alpha = 0.0027. Limits are computed on the original proportion scale by the normal approximation, the first Cornish-Fisher correction (cf1), or the second (cf2); the actual false alarm is evaluated by exact binomial tail sums. The pooled estimator replaces the unweighted mean when subgroup sizes differ.

The shared published fixtures live in inst/extdata/validation/p_chart_joekes_barbosa_2013.csv; the vignette reads them with system.file() and recomputes the limits through the public API.

p_fixtures <- read.csv(
  system.file("extdata", "validation",
              "p_chart_joekes_barbosa_2013.csv",
              package = "IQCC"),
  stringsAsFactors = FALSE
)
p_out <- do.call(rbind, lapply(seq_len(nrow(p_fixtures)), function(i) {
  f <- p_fixtures[i, ]
  lim <- pchart_limits(p = f$p, n = f$n, type = f$method)
  data.frame(
    table = f$table,
    method = f$method,
    n = f$n,
    p = f$p,
    published_ucl = f$published_ucl,
    calculated_ucl = round(lim$ucl, 4),
    ucl_ratio = round(abs(lim$ucl - f$published_ucl) / f$tol_ucl, 3),
    evidence_type = f$evidence_type
  )
}))
p_out
#>     table method  n     p published_ucl calculated_ucl ucl_ratio
#> 1 Table 2 normal 20 0.015        0.0965         0.0965     0.785
#> 2 Table 2    cf1 20 0.015        0.1612         0.1612     0.096
#> 3 Table 2    cf2 20 0.015        0.1303         0.1303     0.385
#> 4 Table 3 normal 20 0.004        0.0463         0.0463     0.823
#> 5 Table 3    cf1 20 0.004        0.1125         0.1125     0.533
#> 6 Table 3    cf2 20 0.004        0.0533         0.0533     0.269
#>     evidence_type
#> 1 published_table
#> 2 published_table
#> 3 published_table
#> 4 published_table
#> 5 published_table
#> 6 published_table

The exact binomial risk evaluation reproduces the published conditional risks of Tables 2 and 3 with tolerance ratio at most one. The maximum observed ucl_ratio across the six fixture cells is 0.82, within the half-unit-in-last-decimal convention of the source.

Practical orientation. Use cf2 when n p (1 - p) is small enough that the discrete and skewed binomial cannot be trusted to a normal limit; the exact risk function pchart_alpha_risk() reports what the chosen limits actually achieve.

Exact range chart

Source and objective. Barbosa, E. P., Gneri, M. A. and Meneguetti, A. (2013). Range Control Charts Revisited: Simpler Tippett-like Formulae, Its Practical Implementation, and the Study of False Alarm. Communications in Statistics - Simulation and Computation, 42(2), 247-262. doi:10.1080/03610918.2011.639967.

Parametrization. The relative range W = R / sigma has the Studentized range distribution. Exact probability limits use the quantiles qtukey(alpha/2, n, Inf) and qtukey(1 - alpha/2, n, Inf). The publication’s Tables 1a and 1b report the false-alarm risk of the normal-approximation chart whose limits are d2 +/- z d3 with z = 3 (alpha0 = 0.0027) and z = 3.09 (alpha0 = 0.0020); the source text states that the second table substitutes 3.09 for 3.

# Representative Table 2 rows (all 54 quantiles are covered by the suite;
# three rows are reproduced here as an executable summary).
r_probs <- c(0.0010, 0.00135, 0.9980, 0.9973, 0.9990, 0.99865)
r_published <- rbind(
  `2` = c(0.00177, 0.00239, 4.37025, 4.24261, 4.65351, 4.53274),
  `5` = c(0.36739, 0.39653, 5.23478, 5.12314, 5.48375, 5.37740),
  `10` = c(1.08458, 1.12634, 5.74143, 5.63772, 5.97331, 5.87416)
)
r_out <- data.frame(
  n = rep(c(2L, 5L, 10L), each = length(r_probs)),
  probability = rep(r_probs, 3),
  published = as.vector(t(r_published)),
  calculated = as.vector(mapply(
    qtukey, p = rep(r_probs, 3),
    nmeans = rep(c(2, 5, 10), each = length(r_probs)),
    MoreArgs = list(df = Inf)
  ))
)
r_out$tolerance_ratio <- round(
  abs(r_out$calculated - r_out$published) / 1e-5, 3
)
r_out
#>     n probability published  calculated tolerance_ratio
#> 1   2     0.00100   0.00177 0.001772454           0.245
#> 2   2     0.00135   0.00239 0.002392814           0.281
#> 3   2     0.99800   4.37025 4.370248438           0.156
#> 4   2     0.99730   4.24261 4.242608150           0.185
#> 5   2     0.99900   4.65351 4.653507531           0.247
#> 6   2     0.99865   4.53274 4.532742813           0.281
#> 7   5     0.00100   0.36739 0.367392008           0.201
#> 8   5     0.00135   0.39653 0.396528087           0.191
#> 9   5     0.99800   5.23478 5.234783685           0.369
#> 10  5     0.99730   5.12314 5.123140141           0.014
#> 11  5     0.99900   5.48375 5.483753781           0.378
#> 12  5     0.99865   5.37740 5.377402382           0.238
#> 13 10     0.00100   1.08458 1.084582650           0.265
#> 14 10     0.00135   1.12634 1.126343057           0.306
#> 15 10     0.99800   5.74143 5.741432384           0.238
#> 16 10     0.99730   5.63772 5.637723511           0.351
#> 17 10     0.99900   5.97331 5.973306528           0.347
#> 18 10     0.99865   5.87416 5.874157502           0.250

The max_tolerance_ratio of the R-chart suite is 0.378 for the printed five-decimal quantiles; the constants d2 and d3 are additionally validated by an independent tail-integral route, and the three-sigma risk inflation of Table 1a is reproduced without presenting the percentage columns as primary fixtures.

Double-sampling np chart

Source and objective. Joekes, S., Smrekar, M. and Barbosa, E. P. (2015). Extending a double sampling control chart for non-conforming proportion in high quality processes to the case of small samples. Statistical Methodology, 23, 35-49. doi:10.1016/j.stamet.2014.09.003.

Parametrization. Two-stage sampling with first and second sample sizes n1, n2, warning line WL, first-stage limit UCL1, and combined second-stage limit UCL2. Acceptance, signal, and continuation probabilities are built from independent binomial components; ASS uses the published full-second-sample convention of the article, and a separate curtailed convention never inflates ASS.

dsnp_fixture <- data.frame(
  row = c(
    "Table 2, p0 = 0.5%, n = 80",
    "Table 2, p0 = 1%, n = 40",
    "Table 3, p0 = 0.5%, n = 80",
    "Table 4, p0 = 0.5%, n = 60"
  ),
  p0 = c(0.005, 0.010, 0.005, 0.005),
  p1 = c(0.0075, 0.0150, 0.010, 0.0075),
  n1 = c(64L, 32L, 65L, 50L),
  n2 = c(271L, 137L, 267L, 242L),
  wl = 1.5,
  ucl1 = c(3.5, 3.5, 3.5, 2.5),
  ucl2 = c(5.5, 5.5, 5.5, 4.5),
  published_arl0 = c(371.19, 371.55, 372.08, 200.04),
  published_arl1 = c(64.00, 63.67, 21.42, 51.35)
)
dsnp_fixture$calculated_arl0 <- vapply(
  seq_len(nrow(dsnp_fixture)),
  function(i) {
    f <- dsnp_fixture[i, ]
    round(dsnp_arl(f$p0, f$n1, f$n2, f$wl, f$ucl1, f$ucl2)$arl, 2)
  },
  numeric(1)
)
dsnp_fixture$calculated_arl1 <- vapply(
  seq_len(nrow(dsnp_fixture)),
  function(i) {
    f <- dsnp_fixture[i, ]
    round(dsnp_arl(f$p1, f$n1, f$n2, f$wl, f$ucl1, f$ucl2)$arl, 2)
  },
  numeric(1)
)
dsnp_fixture$tolerance_ratio <- pmax(
  round(abs(dsnp_fixture$calculated_arl0 - dsnp_fixture$published_arl0) / 0.005, 3),
  round(abs(dsnp_fixture$calculated_arl1 - dsnp_fixture$published_arl1) / 0.005, 3)
)
dsnp_fixture[, c("row", "published_arl0", "calculated_arl0",
                 "published_arl1", "calculated_arl1", "tolerance_ratio")]
#>                          row published_arl0 calculated_arl0 published_arl1
#> 1 Table 2, p0 = 0.5%, n = 80         371.19          371.19          64.00
#> 2   Table 2, p0 = 1%, n = 40         371.55          371.55          63.67
#> 3 Table 3, p0 = 0.5%, n = 80         372.08          372.08          21.42
#> 4 Table 4, p0 = 0.5%, n = 60         200.04          200.04          51.35
#>   calculated_arl1 tolerance_ratio
#> 1           64.00               0
#> 2           63.67               0
#> 3           21.42               0
#> 4           51.35               0

Nine published design rows across Tables 2, 3 and 4 are reproduced (ASS0, ARL0, ARL1) with tolerance ratio at most one; independent binomial component oracles and an exhaustive small-sample enumeration confirm the production probabilities. dsnp_design() recovers a published plan on an explicitly documented local grid; no global optimum beyond the reported grid is claimed.

u chart

Source and objective. No dedicated published table specific to the CF u-chart has been located by the validation work. The Poisson and Cornish-Fisher derivation is documented in paper/statistical-foundations.md, and the CF1/CF2 limits are validated as an independent derivation from Poisson cumulants. The shared derivation fixtures live in inst/extdata/validation/u_chart_cf_derivation.csv.

u_fixtures <- read.csv(
  system.file("extdata", "validation",
              "u_chart_cf_derivation.csv",
              package = "IQCC"),
  stringsAsFactors = FALSE
)
u_out <- do.call(rbind, lapply(seq_len(nrow(u_fixtures)), function(i) {
  f <- u_fixtures[i, ]
  lim <- uchart_limits(lambda = f$lambda, n = f$n, type = f$method,
                       alpha = f$alpha)
  data.frame(
    method = f$method,
    lambda = f$lambda,
    n = f$n,
    expected_ucl = round(f$expected_ucl, 6),
    calculated_ucl = round(lim$ucl, 6),
    ratio = round(abs(lim$ucl - f$expected_ucl) / 1e-10, 0),
    evidence_type = f$evidence_type
  )
}))
u_out
#>    method lambda  n expected_ucl calculated_ucl ratio          evidence_type
#> 1  normal   0.05  5     0.349998       0.349998     0 independent_derivation
#> 2     cf1   0.05  5     0.616660       0.616660     0 independent_derivation
#> 3     cf2   0.05  5     0.483330       0.483330     0 independent_derivation
#> 4  normal   0.10 10     0.399998       0.399998     0 independent_derivation
#> 5     cf1   0.10 10     0.533329       0.533329     0 independent_derivation
#> 6     cf2   0.10 10     0.499996       0.499996     0 independent_derivation
#> 7  normal   0.50 20     0.974338       0.974338     0 independent_derivation
#> 8     cf1   0.50 20     1.041004       1.041004     0 independent_derivation
#> 9     cf2   0.50 20     1.035733       1.035733     0 independent_derivation
#> 10 normal   1.40 50     1.901992       1.901992     0 independent_derivation
#> 11    cf1   1.40 50     1.928658       1.928658     0 independent_derivation
#> 12    cf2   1.40 50     1.927862       1.927862     0 independent_derivation

The exact Poisson false-alarm risk is evaluated directly (exact_discrete), and the historical formula behind the legacy CF alias is documented in the code contract rather than presented as a derived limit.

Generalized variance chart

Source and objective. Barbosa, E. P., Gneri, M. A. and Meneguetti, A. Improving Shewhart-type Generalized Variance Control Charts for Multivariate Process Variability Monitoring using Cornish-Fisher Quantile Correction, Meijer-G Function and Other Tools. Research report, IMECC-UNICAMP. No confirmed DOI; none is invented.

Parametrization. For a p-variate normal subgroup of size n, |S| is distributed as a constant times a product of independent chi-square variables. IQCC reproduces the exact p = 2 case in closed form, the published p = 3 exact upper quantile table, moment-matched normal limits, first/second-order Cornish-Fisher limits, and Bartlett-factor simulation.

gv_ref <- paste(
  "Barbosa, Gneri and Meneguetti (IMECC-UNICAMP report)"
)
# Table 1 (p = 3, exact upper): three representative rows.
gv_n <- c(4L, 8L, 15L)
gv_published_0027 <- c(5.370, 5.084, 3.772)
gv_calculated <- vapply(
  gv_n,
  function(n) gv_limits(n, 3, det_sigma = 1, alpha = 0.0027,
                        type = "exact", side = "upper")$ucl,
  numeric(1)
)
gv_t1 <- data.frame(
  table = "Table 1",
  n = gv_n,
  published = gv_published_0027,
  calculated = round(gv_calculated, 3),
  tolerance_ratio = round(abs(gv_calculated - gv_published_0027) / 5e-4, 3),
  evidence_type = "published_table"
)
gv_t1
#>     table  n published calculated tolerance_ratio   evidence_type
#> 1 Table 1  4     5.370      5.370               0 published_table
#> 2 Table 1  8     5.084      5.084               0 published_table
#> 3 Table 1 15     3.772      3.772               0 published_table

Unresolved parametrizations in the source report. Several aspects of the source cannot be pinned down from the text and are recorded honestly rather than forced:

The published two-sided CF risk column of Table 4 is not reproduced by either the first-order or the second-order Cornish-Fisher expansion. The order used in the column is not disambiguated in the source: order 1 overshoots for small n, and order 2 undershoots for every n. The production expansion matches the corresponding independent oracle, so the column is left as an unresolved order rather than imputed an error.
The one-sided “normal” columns of Tables 5 and 6 describe three-sigma limits (z = 3), whereas production’s one-sided normal chart uses qnorm(1 - alpha). The published values match the three-sigma oracle; this is a documented convention difference, not an error.
The published normal and CF columns of Table 6 are numerically consistent with limits at det_sigma = 0.5 (equivalently, half the production scale at det_sigma = 1). Because the source does not state which |Sigma| was used in that simulation, this is recorded as an implicit/unresolved scale convention rather than imputed an error to the source.

The Exact(sim.) column of Table 6 is Monte Carlo evidence (approximately 10^6 Wishart draws, no documented seed); it is reproduced only in order of magnitude and is never called exact. The generic product-of-chi-square quantile engine remains open as a separate research item.

The standardized first-order CF quantiles of Table 2 are reproduced for n = 15:30 within the harmonized tolerance policy (half-unit by default, one-unit exception only for the documented borderline cells).

Validation coverage and remaining gaps

coverage <- data.frame(
  family = c("p chart", "R chart", "DS-np", "u chart", "|S|"),
  primary_source = c(
    "Joekes & Barbosa (2013)",
    "Barbosa et al. (2013)",
    "Joekes et al. (2015)",
    "Poisson cumulant derivation",
    "IMECC-UNICAMP report"
  ),
  published_tables = c("2, 3", "1a, 1b, 2", "2-4", "-", "1, 2, 4, 5, 6*"),
  published_cells = c("12", "88", "27", "-", "deterministic"),
  independent_checks = c("binomial risk oracle", "Tippett CDF, d2/d3",
                         "component oracles, enumeration", "Poisson risk",
                         "moment oracle, p = 2 CDF"),
  max_tolerance_ratio = c(
    sprintf("%.2f", max(p_out$ucl_ratio)),
    sprintf("%.3f", max(r_out$tolerance_ratio)),
    sprintf("%.2f", max(dsnp_fixture$tolerance_ratio)),
    "exact",
    "0.97"
  ),
  status = c("reproduced", "reproduced", "reproduced",
             "derived, no published table",
             "reproduced with unresolved parametrizations"),
  remaining_issue = c("-", "-", "-", "-",
                      "Table 4 CF order and Table 6 scale convention")
)
coverage
#>    family              primary_source published_tables published_cells
#> 1 p chart     Joekes & Barbosa (2013)             2, 3              12
#> 2 R chart       Barbosa et al. (2013)        1a, 1b, 2              88
#> 3   DS-np        Joekes et al. (2015)              2-4              27
#> 4 u chart Poisson cumulant derivation                -               -
#> 5     |S|        IMECC-UNICAMP report   1, 2, 4, 5, 6*   deterministic
#>               independent_checks max_tolerance_ratio
#> 1           binomial risk oracle                0.82
#> 2             Tippett CDF, d2/d3               0.378
#> 3 component oracles, enumeration                0.00
#> 4                   Poisson risk               exact
#> 5       moment oracle, p = 2 CDF                0.97
#>                                        status
#> 1                                  reproduced
#> 2                                  reproduced
#> 3                                  reproduced
#> 4                 derived, no published table
#> 5 reproduced with unresolved parametrizations
#>                                 remaining_issue
#> 1                                             -
#> 2                                             -
#> 3                                             -
#> 4                                             -
#> 5 Table 4 CF order and Table 6 scale convention

* Table 6 normal and CF columns are reproduced on the implicit det_sigma = 0.5 scale (an unresolved source convention); the Exact(sim.) column is Monte Carlo. The status column is derived from the executed checks above, not from a narrative.

cat("<!-- IQCC_EXECUTED_NUMERICAL_VALIDATION -->\n")

References

Joekes, S. and Barbosa, E. P. (2013). An improved attribute control chart for monitoring non-conforming proportion in high quality processes. Control Engineering Practice, 21, 407-412.
Barbosa, E. P., Gneri, M. A. and Meneguetti, A. (2013). Range Control Charts Revisited: Simpler Tippett-like Formulae, Its Practical Implementation, and the Study of False Alarm. Communications in Statistics - Simulation and Computation, 42(2), 247-262.
Joekes, S., Smrekar, M. and Barbosa, E. P. (2015). Extending a double sampling control chart for non-conforming proportion in high quality processes to the case of small samples. Statistical Methodology, 23, 35-49.
Barbosa, E. P., Gneri, M. A. and Meneguetti, A. Improving Shewhart-type Generalized Variance Control Charts for Multivariate Process Variability Monitoring using Cornish-Fisher Quantile Correction, Meijer-G Function and Other Tools. Research report, IMECC-UNICAMP.

Flavio Barros and contributors

2026-08-12