# Average Sample Size For Curtailed Double Sampling Plan

I am currently implementing double sampling. The goal is to accept or reject a lot by only sampling a few instances from it. In a two-stage testing procedure (double sampling) one first draws a sample of size $$n_1$$ and compares the number $$K_1$$ of items “of interest” (e.g., “non compliant”) with two integers: $$c_1$$, $$c_2$$. If $$K_1 \le c_1$$, we accept the null hypotheses ($$H_0$$) that the batch is of acceptable quality; and if $$X > c_2$$, we accept the alternative hypotheses ($$H_1$$) that the batch quality is unacceptable. However, if $$c_1 < X \le c_2$$, we draw a second sample of size $$n_2$$. We stop in case we find more than $$c_2$$. You can play with some plans in here.

My problem is now finding the average sample size (ASN) with partial curtailment, that is we look at the first sample as a whole and then stop looking at the first sample once it is clear that we reject it.

I found this paper that computes it as

$$ASN = n_1 + \sum_{i=c_1 + 1}^{c2} P(n_1, i) [ n_2 [ P_1(n_2, c_2 - i) + \frac{c_2 - i + 1}{p} P_2(n + 1, c_2 - i + 2) ]]$$

where $$P(n_1, i)$$ is the probability of observing exactly $$j$$ defectives in a sample of size $$n_1$$, $$P_1(n_2, c_2 − i)$$ is the probability of observing $$c_2 − i$$ or fewer defectives in a sample of size $$n_2$$, and $$P_2(n_2 + 1, c_2 − i + 2)$$ is the probability of observing $$c_2 − i + 2$$ defectives in a sample of size $$n_2 + 1$$.

From when I read it, then I would think $$P$$, $$P_1$$ and $$P_2$$ are all PMFs, e.g. of Bernoulli $$p$$, but when I implement it, then I get very different curves for most of the plans compared to here.

I also had a look at

The Average Sample Number for Truncated Single and Double Attributes Acceptance Sampling Plans Author(s): C. C. Craig Source: Technometrics, Vol. 10, No. 4 (Nov., 1968), pp. 685-692

and at the end it looks like $$P_1$$ and $$P_2$$ are CDFs, but I do not manage to get the same results as their tables. Can you help me understand what I need to put in as $$P$$, $$P_1$$ and $$P_2$$?

• The factor of $1/p$ looks fishy, because as $p$ decreases, the ASN ought to decrease, but that factor makes it explode. Your source doesn't provide a derivation--it only gives a reference but omits the bibliography. As far as I can tell, $P, P_1,$ and $P_2$ are all the same function and are explicitly referred to as Binomial probabilities, of the kinds appearing in formulas (5) through (8).
– whuber
Commented Aug 9, 2022 at 22:05
• @whuber Thank you for the comment! The second source has a derivation, but I find the P_2 most difficult to understand as they use a negative binomial or so to "simplify" from what I understood. Regarding $p$, I think that it does not explode because the terms before go close to zero, as having a low chance for defects would lead to having the second step not occur, leaving us with only sampling $n_1$ instances. Commented Aug 9, 2022 at 22:25

This comes down to an analysis of "curtailment" (which is a simplified version of sequential sampling). Let's focus on that first. Because I think your reference is incorrect, I will be careful and show the details.

Curtailment occurs during the second stage of double sampling of a process that has an (assumed) probability $$q \gt 0$$ at each step of yielding a "defective" result. In this stage we plan to sample for as many as $$n$$ steps, but we will stop sampling at any step where the cumulative number of defectives reaches a threshold $$h \gt 1.$$ What is the expected sample size for this stage?

It helps to be formal because we need to track two random variables:

• Let $$X(i)$$ be the number of defectives observed on or before step $$i.$$
• Let $$N$$ be the number of steps taken to first observe $$h$$ defectives with unlimited sampling. It is related to $$X(i)$$ by $$N = \min\{i\mid X(i) \ge h\}.$$
• Write $$N\wedge n = \min(N, n)$$ for the number of samples limited by $$n.$$

Even if we were to terminate sampling early ($$N \lt n$$), it would not change the results of our calculations if we went ahead and obtained a full sample of size $$n.$$ Breaking up the possible outcomes according to the values of $$X(n),$$ we may express the expected sample size as a weighted average of the expectations for each of those values:

$$E[N\wedge n] = \sum_{k = 0}^n E[N \wedge n\mid X(n) = k] \Pr(X(n) = k).\tag{*}$$

The probabilities are Binomial: $$\Pr(X(n)=k)$$ is the chance of seeing exactly $$k$$ defectives in any sample of size $$n.$$ A formula is

$$\Pr(X(n) = k) = \binom{n}{k} q^k(1-q)^{n-k} = \frac{n!}{k!(n-k)!} q^k(1-q)^{n-k}.$$

We must figure out the conditional expectations in $$(*),$$ which for brevity I will write as

$$g(n,h,k) = E[N \wedge n\mid X(n) = k].$$

That is, knowing there are $$k$$ defectives among the next $$n$$ observations, when will we stop sampling? Clearly, when $$k \lt h$$ we will obtain all $$n$$ samples. Otherwise, when $$k\ge h,$$ we will stop sampling by step $$n$$ (and likely sooner than that).

Consider the first observation among those $$n.$$ There are two possibilities: it is defective or not. The key idea is that because there are $$k$$ defectives among the next $$n$$ observations, the chance that the next observation is defective is $$k/n.$$ (The information given by the numbers $$(n,k)$$ does not determine where in the sequence $$1,2,\ldots, n$$ the defectives might occur: all subsets of $$k$$ of those positions are equally likely.) Consequently

$$g(n,h,k) = 1 + \frac{i}{n} g(n-1,h-1,k-1) + \frac{n-i}{n} g(n-1, h, k).$$

The initial $$1$$ counts the first observation and the next two terms use the law of conditional expectations to find the expected number of additional observations: $$g(n-1,h-1,k-1)$$ is the expected sample size with $$n-1$$ observations left to go; the threshold to terminate sampling has been reduced to $$h-1$$ because we just saw a defective; and there remain $$k-1$$ defectives among these observations. $$g(n-1, h, k)$$ is the expected sample size with $$n-1$$ observations left to go and we haven't yet seen a defective.

The (unique) solution to this recursion (consistent with obvious starting conditions) is

$$g(n, h, k) = \frac{h(n+1)}{k+1}.\tag{**}$$

To see why, note that the formula is correct in all edge cases; e.g., when $$k=n,$$ every one of the next $$n$$ observations is a defective and so we stop after seeing the first $$h = h(n+1)/(n+1)$$ of them. Then all we need to check is the recursion, finding

\begin{aligned} g(n, h, k) &= \frac{h(n+1)}{k+1}\\ &= 1 + \frac{k}{n}\frac{(h-1)n}{k} + \frac{n-k}{n} \frac{hn}{k+1}\\ &= 1 + \frac{k}{n} g(n-1,h-1,k-1) + \frac{n-k}{n} g(n-1, h, k). \end{aligned}

with straightforward algebra. Writing $$F(\ ;n,q)$$ for the Binomial$$(n,q)$$ distribution function we may now express $$(*)$$ in the closed form

\begin{aligned} E[N\wedge n] &= n\Pr(X(n)\lt h) + \sum_{k=h}^n \frac{h(n+1)}{k+1} \binom{n}{k}q^k(1-q)^{n-k}\\ &= nF(h-1;n,q) + h\sum_{k=h}^n \binom{n+1}{k+1}q^k(1-q)^{n-k}\\ &= nF(h-1;n,q) + h\sum_{k=h+1}^{n+1} \binom{n+1}{k}q^{k-1}(1-q)^{n+1-k}\\ &= nF(h-1;n,q) + \frac{h}{q}\sum_{k=h+1}^{n+1} \binom{n+1}{k}q^{k}(1-q)^{n+1-k}\\ &= nF(h-1;n,q) + \frac{h}{q}\left(1 - F(h;n+1,q)\right). \end{aligned}.

The first step breaks the expectation into the cases $$X(n)\lt h,$$ where the full sample will be obtained, and the remaining cases $$X(n)\ge h,$$ where early termination is possible. A little algebra and a change in the index from $$k$$ to $$k+1$$ yields the result.

We are ready to compute the expected sample size for double sampling with curtailment. The initial sample of $$n_1$$ observations will be completed, contributing $$n_1$$ to the expectation. If the number of defectives in that sample is $$c_1$$ or less or if it exceeds $$c_2,$$ sampling stops with a size of $$n_1.$$ The lower (red) and upper (purple) paths in the figure illustrate these scenarios.

Otherwise, additional sampling is needed. Sampling terminates either when the cumulative number of defectives exceeds $$c_2$$ (as in the teal path) or the total sample size reaches $$n_1 + n_2$$ (green path). For paths like the teal, $$N$$ (marked with a solid vertical line) can have values between $$n_1$$ and $$n_1 + n_2.$$

The expected extra effort is (as before) computed by partitioning the possibilities according to the value of $$X(n_1),$$ which ranges from $$c_1+1$$ through $$c_2.$$ Let this value be $$j.$$ The threshold to terminate this second stage of sampling with a maximum sample size of $$n=n_2$$ is $$h = c_2+1 - j$$ because $$j$$ defectives have already been observed. Applying the formula $$(**)$$ gives

\begin{aligned} &E[N] \\ &= n_1 + \sum_{j = c_1 + 1}^{c_2} E[N-n_1 \mid X(n_1) = j]\Pr(X(n_1) = j)\\ &= n_1 + \sum_{j = c_1 + 1}^{c_2} \left[nF(c_2-j;n_1,q) + \frac{h}{q}\left(1 - F(c_2+1 - j;n+1,q)\right)\right]\binom{n_1}{j}q^j(1-q)^{n_1-j}. \end{aligned}\tag{***}

Although this looks a lot like formula $$(13)$$ in your reference, it is not the same, because the $$F(\cdots)$$ values are not the Binomial probabilities described there: they are cumulative probabilities and cannot be equated with or converted into individual probabilities. (This makes me suspect the Technometrics article, which I haven't seen, is likely correct.)

Here is a plot of $$E[N]$$ as a function of $$q$$ for the sampling plan in the previous illustration ($$c_1=4,$$ $$c_2 =14,$$ $$n_1=30,$$ and $$n_2=40$$).

This is a general pattern: the sampling effort will be greatest when the rate of defectives is in a "sweet spot" enabling the sample path to stay out of the gray areas. Otherwise it's highly likely sampling will end after the first round of size $$n_1.$$ Even so, the expectation is substantially less than the maximum planned sample size of $$n_1+n_2.$$

To settle the question of whose formula (if either) is correct, we need a third approach. Let's simulate some samples in a way that does not depend on much if any analysis. Small errors will be most apparent in small problems, so here is an example with $$n_1=10,$$ $$n_2=10,$$ $$c_1=3,$$ $$c_2=7,$$ and $$q=0.5.$$

#
# The formula.
#
f <- function(n, h, q) n * pbinom(h-1, n, q) + h * (1 - pbinom(h, n+1, q)) / q
a <- Vectorize(function(c, n, q) {
j <- seq(c[1]+1, c[2])
n[1] + sum(dbinom(j, n[1], q) * f(n[2], c[2]+1-j, q))
}, "q")
#
# Simulation.
#
cc <- c(3, 7)
nn <- c(10,10)
q <- 0.5
set.seed(17)
N <- replicate(5e3, {
k1 <- rbinom(1, nn[1], q)         # Take the first sample
if(k1 <= cc[1] || k1 > cc[2]) nn[1] else {
pop <- rbinom(nn[2], 1, q)      # Take the second sample, one at a time
x <- cumsum(pop) + k1           # Compute the sample path
x[nn[2]] <- cc[2] + 1           # Force termination at the end
nn[1] + match(TRUE, x > cc[2])  # Find where sampling can first terminate
}
})
signif(c(Simulation = mean(N), Formula = a(cc, nn, q)), 4)


The output prints the average sample length in 5,000 simulated samples along with what this formula says:

Simulation    Formula
14.00      14.04


They are statistically indistinguishable. Repeated experiments of this nature bear out the correctness of formula $$(***).$$