---
title: "Statistical Foundations of IQCC"
author: "Flavio Barros and contributors"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Statistical Foundations of IQCC}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
library(IQCC)
```

# Overview

IQCC focuses on control charts whose classical normal limits may be poorly
calibrated because the monitored statistic is discrete, skewed, bounded, or
strongly non-normal. The package combines three ideas:

1. compute limits on the original measurement scale;
2. evaluate actual false-alarm probabilities whenever possible;
3. validate formulas against independent calculations and published examples.

# Binomial p chart

For

\[
X\sim\operatorname{Binomial}(n,p),\qquad \hat p=X/n,
\]

\[
E(\hat p)=p,\qquad
\operatorname{Var}(\hat p)=\frac{p(1-p)}{n}.
\]

The normal limits use

\[
p\pm z\sqrt{\frac{p(1-p)}{n}}.
\]

The first Cornish-Fisher correction adds

\[
\frac{(z^2-1)(1-2p)}{6n}
\]

to the normal quantile on the original p scale.

```{r pchart-example}
pchart_table <- do.call(
  rbind,
  lapply(c("normal", "cf1", "cf2"), function(method) {
    lim <- pchart_limits(p = 0.015, n = 20, type = method)
    risk <- pchart_alpha_risk(
      p = 0.015,
      n = 20,
      lcl = lim$lcl,
      ucl = lim$ucl
    )
    data.frame(
      method = method,
      lcl = lim$lcl,
      center = lim$center,
      ucl = lim$ucl,
      alpha = risk,
      arl0 = ifelse(risk == 0, Inf, 1 / risk)
    )
  })
)
pchart_table
```

For unequal subgroup sizes, IQCC estimates the in-control proportion using the
pooled estimator

\[
\hat p_0=\frac{\sum_i X_i}{\sum_i n_i}.
\]

# Poisson u chart

Let

\[
X\sim\operatorname{Poisson}(\lambda n),\qquad U=X/n.
\]

Then

\[
E(U)=\lambda,\qquad
\operatorname{Var}(U)=\frac{\lambda}{n}.
\]

```{r uchart-example}
u_table <- do.call(
  rbind,
  lapply(c("normal", "cf1", "cf2"), function(method) {
    lim <- uchart_limits(lambda = 0.15, n = 20, type = method)
    risk <- uchart_alpha_risk(
      lambda = 0.15,
      n = 20,
      lcl = lim$lcl,
      ucl = lim$ucl
    )
    data.frame(
      method = method,
      lcl = lim$lcl,
      center = lim$center,
      ucl = lim$ucl,
      alpha = risk,
      arl0 = ifelse(risk == 0, Inf, 1 / risk)
    )
  })
)
u_table
```

# Double-sampling np chart

Let

\[
D_1\sim\operatorname{Binomial}(n_1,p),\qquad
D_2\sim\operatorname{Binomial}(n_2,p),
\]

independently. Define

\[
a=\lfloor WL\rfloor,\qquad
b=\lfloor UCL_1\rfloor+1,\qquad
c=\lfloor UCL_2\rfloor.
\]

The decision rule is:

- accept at stage 1 when \(D_1\leq a\);
- signal at stage 1 when \(D_1\geq b\);
- continue when \(a<D_1<b\);
- accept at stage 2 when \(D_1+D_2\leq c\).

Thus

\[
ARL(p)=\frac{1}{1-P_T(p)},
\qquad
ASS(p)=n_1+n_2P_2(p).
\]

```{r dsnp-example}
prob <- dsnp_prob_accept(
  p = 0.005,
  n1 = 34,
  n2 = 162,
  wl = 1.5,
  ucl1 = 2.5,
  ucl2 = 4.5
)

data.frame(
  pt = prob$pt,
  p_signal = prob$p_signal,
  arl0 = dsnp_arl(0.005, 34, 162, 1.5, 2.5, 4.5)$arl,
  arl1 = dsnp_arl(0.0075, 34, 162, 1.5, 2.5, 4.5)$arl,
  ass0 = dsnp_ass(0.005, 34, 162, 1.5, 2.5)$ass
)
```

# Generalized variance chart

Let \(S\) be the sample covariance matrix from a p-variate normal subgroup of
size \(n\). Then

\[
(n-1)S\sim W_p(n-1,\Sigma).
\]

The monitored statistic is

\[
G=|S|.
\]

Bartlett's decomposition gives

\[
G\overset{d}{=}
\frac{|\Sigma|}{(n-1)^p}
\prod_{j=1}^{p}Y_j,
\qquad
Y_j\sim\chi^2_{n-j},
\]

with independent factors.

```{r gv-example}
gv_table <- do.call(
  rbind,
  lapply(c("normal", "cf", "exact"), function(method) {
    lim <- gv_limits(
      n = 10,
      p = 2,
      det_sigma = 0.5320,
      type = method
    )
    risk <- gv_alpha_risk(
      n = 10,
      p = 2,
      det_sigma = 0.5320,
      type = method
    )
    data.frame(
      method = method,
      lcl = lim$lcl,
      center = lim$center,
      ucl = lim$ucl,
      alpha = risk$alpha,
      arl0 = risk$arl0
    )
  })
)
gv_table
```

```{r gv-figure}
set.seed(321)
x <- array(rnorm(10 * 8 * 2), dim = c(10, 8, 2))
cchart.GV(x, Sigma = diag(2), type = "exact", plot = TRUE)
```

For dimension two, the exact distribution simplifies to

\[
G\overset{d}{=}
\frac{|\Sigma|}{4(n-1)^2}
\left(\chi^2_{2n-4}\right)^2.
\]

For higher dimensions, IQCC offers simulation using the Bartlett factors. The
package does not claim a generic Meijer-G implementation. Published exact
quantiles are used only in explicitly supported dimension-three cases.

# Auxiliary trace chart

The generalized variance chart is sensitive to determinant changes in the
covariance matrix. A complementary diagnostic monitors the trace of the
standardized covariance matrix,

\[
T=(n-1)\operatorname{tr}(\Sigma_0^{-1}S).
\]

If the process is multivariate normal and \(\Sigma=\Sigma_0\), then

\[
T\sim\chi^2_{p(n-1)}.
\]

This gives an exact upper control limit for the trace chart:

```{r trv-limits}
trv_limits(n = 8, p = 2, alpha = 0.0027)
```

## Published Table 3

The report defines \(N=n-1\) in Section 3.3, but the numerical rows of its
Table 3 use the displayed \(N\) as subgroup size \(n\). The probability headings
also appear to be reversed: the values printed under 0.9980 match exact 0.9973
chi-square quantiles, while those under 0.9973 match exact 0.9980 quantiles.
IQCC keeps the mathematical tail convention rather than reproducing the
apparent label error.

```{r trv-published-table-3}
table3 <- data.frame(
  report_row = rep(c(3, 16, 30), each = 2),
  printed_heading = rep(c(0.9980, 0.9973), 3),
  nominal_probability = rep(c(0.9973, 0.9980), 3),
  published = c(20.07, 20.79, 75.88, 77.17, 128.18, 129.83)
)
table3$df <- 3 * (table3$report_row - 1)
table3$calculated <- vapply(
  seq_len(nrow(table3)),
  function(i) trv_limits(
    n = table3$report_row[i],
    p = 3,
    alpha = 1 - table3$nominal_probability[i]
  )$ucl,
  numeric(1)
)
table3$tolerance <- 0.03
table3$tolerance_ratio <- abs(
  table3$calculated - table3$published
) / table3$tolerance
table3
```

## Published Case B: partial reproduction

The report's three-dimensional bolt example uses \(n=15\), 30 Phase I
subgroups, and 40 Phase II subgroups. In Case B, the covariance parameters in
the last row of Table 8 are increased by 20% for Phase II. Tables 9 and 10
publish rounded Phase I mean covariance matrices whose generalized variances
remain close:

```{r trv-published-case-b}
case_a_sbar <- matrix(
  c(4.2366, 1.4773, 1.1929,
    1.4773, 6.1264, 2.3399,
    1.1929, 2.3399, 3.9335),
  nrow = 3, byrow = TRUE
)
case_b_sbar <- matrix(
  c(4.0112, 1.7865, 1.2484,
    1.7865, 5.8933, 1.7188,
    1.2484, 1.7188, 3.8909),
  nrow = 3, byrow = TRUE
)
data.frame(
  case = c("A", "B"),
  published_determinant = c(69.8438, 66.1893),
  determinant_from_rounded_matrix = c(det(case_a_sbar), det(case_b_sbar))
)
data.frame(
  figure = "9c",
  published_trv_ucl = 72.01,
  calculated_trv_ucl = trv_limits(n = 15, p = 3, alpha = 0.0027)$ucl
)
```

The published Case B figure shows that the generalized-variance chart does not
cross its exact or Cornish-Fisher limit, whereas the auxiliary `tr(V)` chart has
Phase II signals. This is a partial, qualitative reproduction: the report
prints only subgroup rows 1, 2, 3, and 30 plus summaries, not all 70 generated
subgroups or an RNG seed. Consequently, the published signal sequence cannot
be recomputed without inventing observations. The deterministic example below
is therefore retained and explicitly labeled synthetic.

The trace chart complements `cchart.GV()` because covariance changes can
preserve the determinant while changing the sum of standardized eigenvalues.
The following synthetic deterministic example has two subgroup covariance
matrices with the same determinant; `|S|` is unchanged, while `tr(V)`
increases.

```{r trv-same-determinant}
base <- sqrt(3 / 4) * rbind(
  c(1, 1), c(1, -1), c(-1, 1), c(-1, -1)
)
g0 <- base %*% chol(diag(2))
g1 <- base %*% chol(diag(c(4, 0.25)))

data.frame(
  generalized_variance = gv_stat(list(g0, g1)),
  trace_statistic = trv_stat(list(g0, g1), Sigma0 = diag(2))
)
```

```{r execution-marker, results='asis'}
cat("<!-- IQCC_EXECUTED_FOUNDATIONS -->\n")
```

# Validation strategy

IQCC combines algebraic property tests, independently coded numerical oracles,
and reproduction of published numerical examples. This is particularly
important because a plausible formula may still use a different tail
convention, integer threshold, or parameter-estimation rule than the intended
method.

# 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.
- 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*.