Unequal-Probability Designs and PPS Sampling

In the preceding units, sampling designs selected elements or primary clusters with equal probability. Under simple random sampling, every unit in the target universe possesses an identical inclusion probability $\pi_i = n / N$. In stratified random sampling with proportional allocation, every element within the universe shares an identical unconditional selection probability. While equal probability designs possess intuitive geometric simplicity and yield self-weighting datasets where unweighted sample averages directly estimate population means, they become extraordinarily inefficient when units vary dramatically in physical size, demographic scale, economic turnover, or administrative capacity.

Consider conducting an establishment survey to estimate total industrial manufacturing output in a country. The target population comprises a small handful of massive multinational conglomerate factories employing tens of thousands of workers alongside hundreds of thousands of micro-enterprises and artisan workshops employing one or two persons. If an equal probability simple random sample of $n = 500$ establishments is selected from this population, the probability that any specific mega-factory is included is tiny. If one mega-factory happens to be selected into the sample, its enormous production volume is multiplied by the massive design weight $N / n$, causing the estimated national total to explode. If, by chance, no mega-factories are selected, the national total is severely underestimated. The sampling distribution of the estimator exhibits extreme skewness and astronomical sampling variance.

The foundational breakthrough of unequal probability sampling—pioneered independently by Morris Hansen and William Hurwitz (1943) in the United States and P. C. Mahalanobis (1944) in India—was the recognition that units should be selected with probabilities proportional to their importance or size. If an auxiliary measure of size $x_i > 0$ is known for every unit $i$ in the sampling frame (such as baseline employee count, historical revenue, school enrollment, or farm acreage), and if this auxiliary size is approximately proportional to the survey variable $y_i$ ($y_i \approx \beta x_i$), then selecting units with probabilities proportional to $x_i$ transforms the estimation problem.

Under Probability Proportional to Size (PPS) selection, large units that dominate the population total are sampled with high probabilities (or even with certainty, $\pi_i = 1$), while tiny units are sampled with appropriately smaller probabilities. When an observation $y_i$ is divided by its selection probability $\pi_i$, the ratio $y_i / \pi_i$ becomes approximately constant across all units in the population:

Because the sampling variance of an estimator depends on the variability of these ratios across the sampled units, making the ratio nearly invariant drives the sampling variance toward zero. Unequal probability sampling transforms auxiliary information from a passive descriptive tool into an active, variance-crushing sampling mechanism.

To implement probability proportional to size sampling, the survey statistician must translate a vector of positive auxiliary sizes $\mathbf{x} = (x_1, x_2, \dots, x_N)^\top$ into an operational random selection protocol. Let:

denote the total aggregate measure of size in the population. The relative size (or single-draw selection probability) of unit $i$ is defined as:

Three classical algorithms are deployed to draw samples according to these probabilities:

1. The Cumulative Total Method

The cumulative total method, introduced by Hansen and Hurwitz (1943), operates by constructing a sequential partition of the interval $[1, X]$:

1. Compute the cumulative sums of the measures of size:

where $T_N = X$.

1. Associate with unit $i$ the discrete half-open integer interval:

The length of interval $I_i$ is exactly $T_i - T_{i-1} = x_i$.

1. Draw a random integer $R$ uniformly distributed on the set $\{1, 2, \dots, X\}$.
2. Select unit $i$ if $R \in I_i$, which occurs if and only if $T_{i-1} < R \le T_i$.

The probability of selecting unit $i$ on any single draw is:

To select a sample of size $n$ with replacement (PPSWR), the procedure is repeated independently $n$ times.

2. Lahiri's Rejection Method

When the population frame contains hundreds of thousands of units, computing and storing the cumulative array $T_i$ can be computationally cumbersome. In 1951, D. B. Lahiri devised an ingenious rejection sampling algorithm that selects units with exact PPS probabilities without ever computing cumulative totals:

1. Determine the maximum measure of size across the population:

1. Draw a random integer $i$ uniformly from $\{1, 2, \dots, N\}$.
2. Draw an independent random integer $R$ uniformly from $\{1, 2, \dots, M_{max}\}$.
3. If $R \le x_i$, accept unit $i$ into the sample. If $R > x_i$, reject the pair $(i, R)$ and repeat steps 2 and 3 until a candidate is accepted.

To verify that Lahiri's algorithm generates exact PPS selection probabilities, we calculate the conditional probability of selecting unit $i$ given that an acceptance occurs:

The overall probability of an acceptance on any single trial is:

Applying Bayes' theorem:

Lahiri's method delivers exact PPS probabilities while requiring only random access to individual $x_i$ values.

3. Systematic PPS Sampling (Cumulative Measure Method)

To select a sample of size $n$ without replacement across a continuous cumulative size scale, Goodman and Kish (1950) formulated systematic PPS selection:

1. Compute the sampling interval based on the aggregate size:

1. Select a single random start $r$ uniformly distributed on the real interval $(0, k]$.
2. Form the sequence of $n$ selection pointers:

1. Select unit $i$ into the sample if any selection pointer falls within its cumulative interval:

Under this scheme, the expected number of pointers hitting unit $i$ is:

Provided that no individual unit exceeds the interval length ($x_i \le k$, or equivalently $n p_i \le 1$), a unit can be struck by at most one pointer. Under this condition, the design is strictly without replacement, and the first-order inclusion probability is exact:

Units whose size exceeds the interval ($x_i > k$, meaning $n p_i > 1$) would be selected multiple times. In survey practice, such dominant units are identified prior to sampling, removed from the frame, designated as certainty units (or self-representing strata) with $\pi_i = 1$, and enumerated completely. Systematic PPS is then executed on the remaining non-certainty units.

When sampling is conducted with replacement, the draws are mutually independent, which vastly simplifies mathematical derivations. Although sampling with replacement is slightly less efficient than sampling without replacement, the Hansen–Hurwitz estimator serves as the theoretical benchmark for all unequal probability designs.

Formulation of the Estimator

Let $n$ independent draws be made from the population with selection probabilities $p_i = x_i / X$. Let $y_{(k)}$ and $p_{(k)}$ denote the survey attribute and selection probability of the unit selected on the $k$-th draw ($k = 1, 2, \dots, n$).

For each draw $k$, we define the single-draw estimator of the population total $Y$:

The Hansen–Hurwitz Estimator of the population total $Y$ is the arithmetic average of these $n$ independent single-draw estimators:

To estimate the population mean per element $\bar{Y} = Y / M_0$ (where $M_0$ is the known total number of elements):

Proof of Design Unbiasedness

On any single draw $k$, the random variable $\hat{Y}_k$ takes the value $\frac{y_i}{p_i}$ with probability $p_i$ ($i = 1, 2, \dots, N$). The expected value of $\hat{Y}_k$ is:

Therefore, every single-draw estimator is strictly unbiased.

Taking the expectation of the Hansen–Hurwitz estimator:

This proves that the Hansen–Hurwitz estimator is strictly design-unbiased under PPSWR.

Derivation of the Exact Sampling Variance

Because the $n$ draws are mutually independent, the variance of the average is the variance of a single draw divided by $n$:

The variance of a single draw $\hat{Y}_1$ is:

Expanding the square:

Since $\sum y_i = Y$ and $\sum p_i = 1$:

Therefore, the exact variance of the Hansen–Hurwitz estimator is:

The Pairwise Algebraic Identity

The variance of $\hat{Y}_{HH}$ can also be expressed in a deeply instructive pairwise form:

To prove this identity:

Since $\sum p_j = 1$ and $\sum y_i = Y$:

Dividing by $n$ completes the proof.

This pairwise formulation highlights the governing principle of PPS sampling: the variance depends directly on the differences between the ratios $\frac{y_i}{p_i}$ and $\frac{y_j}{p_j}$. If the survey variable $y$ is strictly proportional to the auxiliary measure $x$ ($y_i \propto x_i \implies y_i / p_i = \text{constant}$), then every difference $\left( \frac{y_i}{p_i} - \frac{y_j}{p_j} \right)$ is identically zero, and the sampling variance vanishes completely!

Unbiased Estimation of Variance

In practice, the true population variance $V(\hat{Y}_{HH})$ is unknown. Because the $n$ values $\hat{Y}_k = \frac{y_{(k)}}{p_{(k)}}$ are independent and identically distributed random variables with mean $Y$ and variance $n V(\hat{Y}_{HH})$, the sample variance among them provides an immediate, strictly unbiased estimator of the variance.

Define the sample variance:

Because $s_{HH}^2$ is the standard sample variance of $n$ independent observations, classical statistical theory guarantees that:

Therefore, the strictly unbiased estimator of the variance of the Hansen–Hurwitz estimator is:

Notice the profound operational elegance of this result: computing the estimated variance requires no knowledge of population totals or second-order joint inclusion probabilities. It requires only the $n$ observed values and their single-draw selection probabilities $p_{(k)}$.

While sampling with replacement provides mathematical tractability, drawing units with replacement is inherently wasteful because duplicate selections provide zero additional empirical information about the unobserved elements of the finite population. Sampling without replacement is always statistically more efficient.

In 1952, D. G. Horvitz and D. J. Thompson formulated the universal mathematical theory for estimating finite population parameters from any probability design without replacement, regardless of how complex the selection mechanism may be.

Inclusion Probabilities

In sampling without replacement, an element $i$ can appear at most once in the sample $s$. The probabilistic behavior of the sample is governed by its inclusion probabilities:

1. First-Order Inclusion Probability:

where $I_i$ is the sample inclusion indicator:

The fundamental design requirement is that $\pi_i > 0$ for all $i \in U$.

2. Second-Order (Joint) Inclusion Probability:

Notice that for $i = j$, $I_i^2 = I_i$, so $\pi_{ii} = \pi_i$.

For any sampling design with a fixed sample size $n$ ($\sum_{i=1}^N I_i = n$), the inclusion probabilities satisfy two strict conservation laws:

The Horvitz–Thompson Estimator

The Horvitz–Thompson (HT) Estimator of the population total $Y$ is defined as:

The inverse of the inclusion probability $d_i = \frac{1}{\pi_i}$ is the design weight (or expansion factor) of unit $i$. It represents the number of units in the target population represented by sampled unit $i$.

Proof of Design Unbiasedness:

Taking the expectation over the randomization distribution:

The Horvitz–Thompson estimator is strictly design-unbiased for any probability design, provided only that all $\pi_i > 0$.

The Horvitz–Thompson Variance Formula

To derive the variance of $\hat{Y}_{HT}$, we evaluate the covariance between the sample inclusion indicators $I_i$ and $I_j$:

Expressing the variance of the linear combination:

Substituting the covariance expressions yields the classical Horvitz–Thompson Variance Formula:

The Horvitz–Thompson Sample Variance Estimator

By substituting the unbiased estimators of the indicator products into the variance formula, Horvitz and Thompson derived the sample variance estimator:

Because $E\left[ \frac{I_i}{\pi_i} \right] = 1$ and $E\left[ \frac{I_i I_j}{\pi_{ij}} \right] = 1$, it is mathematically guaranteed that:

The Fatal Flaw of the HT Variance Estimator:

Despite its theoretical unbiasedness, $\hat{V}_{HT}$ suffers from a severe practical pathology: for many empirical samples, $\hat{V}_{HT}$ can evaluate to a negative number. A negative variance estimate is completely meaningless and invalidates the construction of confidence intervals. This pathology occurs because when the cross-product terms $\pi_{ij} - \pi_i \pi_j$ are negative (which is common under without-replacement designs where the selection of one unit depresses the joint selection probability of another), a sample containing large positive values $y_i y_j$ can cause the negative double sum to overpower the positive single sum.

To rescue variance estimation from the negative-variance pathology, A. R. Sen (1953) in India and F. Yates and P. M. Grundy (1953) in the United Kingdom independently discovered an alternative, vastly superior formulation of the variance and its sample estimator, applicable to all sampling designs with a fixed sample size $n$.

The Sen–Yates–Grundy Algebraic Identity

For any design where the sample size is fixed ($\sum_{i=1}^N I_i = n$), we previously established the joint probability identity:

Subtracting $\sum_{j \ne i}^N \pi_i \pi_j = \pi_i (n - \pi_i)$ from both sides:

Therefore, we have the critical relation:

Dividing by $\pi_i^2$:

Derivation of the SYG Variance Formula

Now substitute this identity into the first term of the Horvitz–Thompson variance:

Symmetrizing over all pairs $(i, j)$ with $i \ne j$:

Now consider the second term of the Horvitz–Thompson variance:

Adding the two components together:

Recognizing the perfect square inside the bracket, we arrive at the celebrated Sen–Yates–Grundy (SYG) Variance Formula:

The Sen–Yates–Grundy Sample Variance Estimator

Applying the sample inclusion indicators to the SYG formula yields the sample variance estimator:

Mathematical Superiority and Non-Negativity:

1. Guaranteed Non-Negativity: The squared difference $\left( \frac{y_i}{\pi_i} - \frac{y_j}{\pi_j} \right)^2$ is strictly non-negative. Therefore, $\hat{V}_{SYG}$ is guaranteed to be strictly non-negative whenever the design satisfies:

This condition is satisfied by all well-behaved probability sampling designs without replacement (such as Brewer's, Durbin's, and Sampford's methods).

2. Zero-Variance Property: If the survey variable is strictly proportional to size ($y_i \propto \pi_i$), every term $\left( \frac{y_i}{\pi_i} - \frac{y_j}{\pi_j} \right)$ is identically zero. Thus, the sample variance estimator $\hat{V}_{SYG}$ evaluates to exactly zero for every single sample, reflecting the true zero variance of the design. The naive Horvitz–Thompson estimator does not possess this property.

Executing PPS sampling without replacement for an arbitrary sample size $n$ requires specifying a selection protocol that achieves $\pi_i = n p_i$ while guaranteeing tractable, strictly positive joint inclusion probabilities $\pi_{ij} > 0$. Several landmark methods achieve this:

1. Brewer's Method for $n = 2$

For a sample size of $n = 2$, K. R. W. Brewer (1963) designed a protocol that guarantees exact $\pi_i = 2 p_i$:

1. The first unit is drawn with selection probability:

where $D = \frac{1}{2} \sum_{k=1}^N \frac{p_k (1 - p_k)}{1 - 2 p_k}$.

1. The second unit is drawn from the remaining $N - 1$ units with probability proportional to its original size:

Brewer proved that under this two-step selection:

• The first-order inclusion probability is exact: $\pi_i = 2 p_i$.
• The joint inclusion probability is:

• Furthermore, $\pi_{ij} < \pi_i \pi_j$ holds universally, guaranteeing that the Sen–Yates–Grundy variance estimator is strictly non-negative for every possible sample!

2. Durbin's Method for $n = 2$

J. Durbin (1967) proposed an alternative two-draw procedure:

1. The first unit is selected with probability $p_i$.
2. If unit $i$ is selected on the first draw, the second unit $j$ ($j \ne i$) is selected with conditional probability:

Durbin's method produces the exact same joint inclusion probabilities $\pi_{ij}$ as Brewer's method, confirming their underlying structural equivalence.

3. Systematic PPS Sampling and Variance Approximations

For general sample sizes $n > 2$, systematic PPS selection (arranging units along the cumulative size line and stepping by interval $k = X/n$) is universally favored in operational fieldwork due to its speed and simplicity.

However, systematic PPS has an analytical drawback: because the selection uses a single random start, many joint inclusion probabilities $\pi_{ij}$ are identically zero (units separated by non-interval distances can never appear in the same sample). When some $\pi_{ij} = 0$, an unbiased estimator of variance does not exist!

To solve this problem, survey statisticians deploy two robust variance approximations:

1. The Hartley–Rao Variance Approximation (1962):

Assuming the population list is randomly ordered prior to systematic selection, Hartley and Rao derived the asymptotic variance:

which behaves similarly to PPSWR variance reduced by a finite-population-like factor $\left( 1 - \frac{n-1}{n}\pi_i \right)$.

2. The Successive Difference Replication Estimator:

When the sample units are ordered as selected along the systematic path, $s = \{1, 2, \dots, n\}$, local successive differences eliminate linear trends and spatial gradients:

This estimator is widely used by national census bureaus and statistical agencies.

While PPS sampling without replacement using the Horvitz–Thompson estimator is statistically efficient, calculating exact second-order joint inclusion probabilities $\pi_{ij}$ for general sample sizes $n > 2$ is computationally intricate. Furthermore, as proven in previous sections, systematic PPS sampling suffers from zero joint inclusion probabilities ($\pi_{ij} = 0$), precluding strictly unbiased variance estimation.

In 1962, J. N. K. Rao, H. O. Hartley, and W. G. Cochran devised an ingenious, universally celebrated sampling procedure—the Rao–Hartley–Cochran (RHC) Random Group Method. The RHC method achieves three simultaneous breakthroughs:

1. It selects units without replacement with probabilities proportional to size.
2. It yields a variance that is strictly smaller than the with-replacement Hansen–Hurwitz variance.
3. It provides a simple, strictly unbiased variance estimator that requires zero knowledge of joint inclusion probabilities $\pi_{ij}$!

The RHC Selection Protocol

Let a finite population $U$ contain $N$ units with auxiliary sizes $x_i > 0$ and total size $X = \sum_{i=1}^N x_i$. Let $p_i = x_i / X$ denote the normalized size of unit $i$ ($\sum p_i = 1$). We wish to draw a sample of $n$ distinct units.

The RHC selection protocol executes in two distinct phases:

Phase 1: Random Group Partitioning

The finite population of $N$ units is partitioned at random into $n$ mutually exclusive and collectively exhaustive groups $G_1, G_2, \dots, G_n$ of predetermined sizes $N_1, N_2, \dots, N_n$, such that:

To maximize efficiency, the group sizes $N_g$ should be as equal as possible, such that each $N_g$ is either $\lfloor N/n \rfloor$ or $\lfloor N/n \rfloor + 1$.

Phase 2: Independent Intra-Group PPS Selection

For each random group $g \in \{1, 2, \dots, n\}$: Compute the aggregate group size:

Select exactly one unit from group $G_g$ with probability proportional to its size within that group:

Because one unit is drawn from each of the $n$ disjoint groups, the resulting sample $s$ contains exactly $n$ distinct units without replacement.

The Rao–Hartley–Cochran Estimator

Let $(g)$ denote the single unit selected from group $G_g$, with survey measurement $y_{(g)}$, auxiliary size $x_{(g)}$, and normalized size $p_{(g)} = x_{(g)} / X$.

The Rao–Hartley–Cochran Estimator of the population total $Y$ is:

Proof of Strict Design Unbiasedness

We evaluate the expectation using iterated conditional expectation: first conditioning on the random group partition (Phase 1), and subsequently taking expectation over all possible random partitions.

Let $E_2(\cdot \mid G_1, \dots, G_n)$ denote expectation over the Phase 2 selection conditional on fixed groups:

where $Y_g$ is the total of the survey variable in group $G_g$.

Summing across all $n$ groups:

Now taking the Phase 1 expectation $E_1$ over all possible random group partitions:

This proves that the RHC estimator is strictly design-unbiased.

Exact Sampling Variance Derivation

Applying the Law of Total Variance:

Because $E_2(\hat{Y}_{RHC} \mid G_1, \dots, G_n) = Y$ is a constant, the first term vanishes identically:

Therefore, the entire variance is governed by the second term:

Conditional on the groups, the selections within different groups are completely independent. Thus:

For a single group $g$:

Using the pairwise identity within group $G_g$:

Now take the expectation $E_1$ over all random partitions. Because the partition assigns units to groups completely at random without replacement:

Summing across all $n$ groups:

Substituting this expectation into the pairwise variance formula yields the exact Rao–Hartley–Cochran Variance Formula:

Recall our pairwise formulation for the with-replacement Hansen–Hurwitz variance:

Therefore, the exact relationship connecting RHC variance to Hansen–Hurwitz variance is:

Equal Group Size Simplification:

When $N$ is an exact multiple of $n$, each group has size $N_g = N / n$. Then:

Under equal group sizes, the RHC variance reduces to the remarkably simple expression:

Notice that because $\frac{n - 1}{N - 1} > 0$ for any sample size $n > 1$:

The Rao–Hartley–Cochran method is strictly more efficient than sampling with replacement, capturing a finite-population-like variance reduction of $\left( 1 - \frac{n - 1}{N - 1} \right)$.

Unbiased Estimation of Variance from Sample Data

The paramount practical virtue of the RHC design is its variance estimator. Rao, Hartley, and Cochran proved that an strictly design-unbiased estimator of $V(\hat{Y}_{RHC})$ is:

For equal group sizes $N_g = N / n$, this simplifies to:

This estimator is strictly non-negative for every sample, requires no matrix inversion, and eliminates all joint inclusion probability calculations.

Two additional landmark classical PPS schemes provide vital theoretical foundations:

1. The Midzuno–Sen Sampling Scheme

In 1952, H. Midzuno in Japan and A. R. Sen (1953) in India independently designed an unequal probability selection scheme that guarantees strictly positive joint inclusion probabilities $\pi_{ij} > 0$ for all pairs while achieving exact closed-form inclusion probabilities.

Protocol:

1. Draw 1: Select the first unit from the population with probability proportional to size:

2. Draws 2 through $n$: From the remaining $N - 1$ units, select $n - 1$ units with equal probability without replacement (SRSWOR).

First-Order Inclusion Probability:

A unit $i$ can be selected either on the first draw or during the subsequent SRSWOR draws:

Since $\sum_{j \ne i} p_j = 1 - p_i$:

Notice the structure of $\pi_i$: It is a convex linear combination of the PPS probability $p_i$ and the equal-probability fraction $\frac{n - 1}{N - 1}$. The second-order joint inclusion probability is:

Because $\pi_{ij} > 0$ for all pairs, strictly unbiased variance estimation under the Horvitz–Thompson and Sen–Yates–Grundy formulas is always guaranteed.

2. Murthy's Unordered Estimator

In sequential sampling without replacement (such as Durbin's or Brewer's methods), the units in the sample are drawn in a specific sequential order $(i_1, i_2, \dots, i_n)$. An estimator that depends upon the order of draws is called an ordered estimator (e.g., the Des Raj estimator).

In 1957, M. N. Murthy applied the Rao–Blackwell theorem to survey sampling, proving that conditioning an ordered estimator on the unordered sample set $s = \{i_1, \dots, i_n\}$ eliminates the arbitrary noise of draw order, strictly reducing sampling variance.

For a sample of size $n = 2$ drawn with initial selection probabilities $p_i$, Murthy's unordered estimator of the population total is:

3. The Des Raj Ordered Estimator (1956)

Prior to Murthy's unordered formulation, Des Raj (1956) formulated a sequential estimation theory for probability proportional to size sampling without replacement.

Suppose a sample of $n$ units is drawn sequentially without replacement:

• On Draw 1, unit $i_1$ is selected with probability $p_{i_1}$.
• On Draw 2, unit $i_2$ is selected from the remaining $N - 1$ units with conditional probability $\frac{p_{i_2}}{1 - p_{i_1}}$.
• On Draw $r$ ($r = 2, \dots, n$), unit $i_r$ is selected with conditional probability:

where $P_{r-1} = \sum_{k=1}^{r-1} p_{i_k}$ is the cumulative probability of the units already selected.

Des Raj constructed a sequence of $n$ unbiased estimators $t_1, t_2, \dots, t_n$ of the population total $Y$, one for each sequential draw:

• For Draw 1:

• For Draw 2:

• For Draw $r$:

Proof of Unbiasedness:

Condition on the first $r - 1$ draws:

Because the conditional expectation of $t_r$ is identically $Y$ regardless of the preceding draws, taking unconditional expectation confirms:

The Des Raj Estimator of the population total is the arithmetic average of these $n$ unbiased estimators:

Exact Unbiased Variance Estimator:

Des Raj proved that despite the complex sequential dependency, the sample variance among the $t_r$ values yields a strictly unbiased estimator of the variance of $\hat{Y}_{DR}$:

For $n = 2$:

Notice that for $n = 2$, this variance estimator is strictly non-negative!

For general sample sizes $n > 2$, selecting units without replacement such that $\pi_i = n p_i$ while guaranteeing strictly positive, well-behaved joint inclusion probabilities was solved by Michael R. Sampford in 1967. Sampford's Method is the standard, universal algorithm for PPSWOR of arbitrary sample size.

The Selection Protocol

Let $p_i = x_i / X$ with $\sum p_i = 1$, and assume $n p_i < 1$ for all $i \in U$. Define the modified selection probabilities:

Notice that $\sum_{i=1}^N \lambda_i = 1$.

Sampford's protocol proceeds as follows:

1. Draw 1: Select the first unit $i_1$ with probability proportional to size:

2. Draws 2 through $n$: Select the remaining $n - 1$ units independently with replacement with modified probabilities:

3. Acceptance / Rejection Step:

• If all $n$ drawn units $\{i_1, i_2, \dots, i_n\}$ are distinct (no duplicate units), accept the sample.
• If any unit appears more than once, reject all $n$ units, discard the trial, and restart the entire procedure from Draw 1.

Mathematical Properties of Sampford's Design

1. Exact First-Order Inclusion Probabilities:

Sampford proved that conditional on acceptance, the probability that unit $i$ is included in the sample satisfies:

2. Joint Inclusion Probabilities:

The second-order joint inclusion probabilities are given by:

where $L_m$ are symmetric polynomial functions of $\boldsymbol{\lambda}$, and $K_n$ is a normalizing constant.

3. Guaranteed Non-Negativity of the SYG Variance Estimator:

Sampford proved that for all pairs $i \ne j$:

This guarantees that the Sen–Yates–Grundy variance estimator:

is strictly non-negative for every possible sample drawn under Sampford's scheme.

While the Horvitz–Thompson estimator provides a universal framework for any probability sampling design, implementing unequal probability sampling without replacement (WOR) presents a fundamental technical challenge: for an arbitrary selection protocol, calculating the second-order joint inclusion probabilities $\pi_{ij}$ is often analytically intractable or yields probabilities where $\pi_{ij} > \pi_i \pi_j$, which causes the Sen–Yates–Grundy variance estimator to produce negative variance estimates.

To overcome this, early pioneers formulated specialized without-replacement selection protocols that guarantee exact first-order probabilities $\pi_i = n p_i$ while ensuring that joint probabilities $\pi_{ij}$ satisfy the condition $\pi_{ij} < \pi_i \pi_j$ everywhere, guaranteeing that the Sen–Yates–Grundy variance estimator is strictly non-negative. Two of the most celebrated methods in the literature are Brewer's scheme and Murthy's unordered estimator.

1. Brewer's Unequal Probability Selection Scheme

In 1963, K. R. W. Brewer introduced an elegant selection scheme for drawing a sample of $n = 2$ units without replacement such that the inclusion probability of unit $i$ is strictly proportional to its measure of size:

subject to the physical feasibility condition that $p_i < 1/2$ for all $i$.

Selection Algorithm for $n = 2$:

1. First Draw: Select the first unit $i$ from the population $U$ with revised selection probability:

where the normalizing denominator $D$ is the constant that ensures the probabilities sum to unity:

2. Second Draw: Having selected unit $i$ on the first draw, select the second unit $j$ ($j \ne i$) from the remaining $N - 1$ units with probability strictly proportional to its original size:

Proof of Exact Inclusion Probabilities:

The probability that unit $i$ is selected on the first draw is $p_i^{(1)}$. The probability that unit $i$ is selected on the second draw given that some other unit $j$ was selected first is:

The unconditional first-order inclusion probability of unit $i$ is the sum of its chances of being selected on either draw:

Now substitute the definition of $D = \frac{1}{2} \left( 1 + \sum_{k=1}^N \frac{p_k}{1 - 2 p_k} \right)$:

This proves that Brewer's method yields strictly exact first-order PPS inclusion probabilities $\pi_i = 2 p_i$ without any approximation error.

Joint Inclusion Probabilities $\pi_{ij}$:

The joint probability that both unit $i$ and unit $j$ are included in the sample is:

Brewer proved that under this formulation, the inequality $\pi_{ij} < \pi_i \pi_j = 4 p_i p_j$ holds universally for all distinct pairs $(i, j)$. Consequently, the Sen–Yates–Grundy variance estimator is guaranteed to be strictly positive for all realized samples.

2. Murthy's Unordered Estimator

In standard Horvitz–Thompson estimation, the order in which units are drawn is ignored, but the estimator relies heavily on knowing the joint inclusion probabilities $\pi_{ij}$. In 1956, Des Raj proposed an ordered estimator where observations are weighted according to their sequential draw probabilities. While Des Raj's estimator is strictly unbiased, it is inefficient because it depends on the order of the sample draws.

Recognizing this deficiency, V. K. Murthy (1957) applied the Rao–Blackwell theorem to Des Raj's estimator: by conditioning on the unordered set of sampled units $s = \{i_1, i_2, \dots, i_n\}$, Murthy eliminated the order dependence, producing a uniformly minimum variance unordered estimator that does not require calculating $\pi_{ij}$.

Murthy's Estimator for $n = 2$:

Consider drawing $n = 2$ units without replacement.

• Draw 1: Unit $i$ is selected with probability $p_i$.
• Draw 2: Unit $j$ is selected from the remaining units with conditional probability $p_{j|i} = \frac{p_j}{1 - p_i}$.

The joint probability of drawing unit $i$ first and unit $j$ second is:

The joint probability of drawing unit $j$ first and unit $i$ second is:

The unconditional probability of obtaining the unordered sample pair $s = \{i, j\}$ is:

By Rao–Blackwellization, Murthy's unordered estimator for the population total $Y$ based on the realized pair $(y_1, y_2)$ with size measures $(p_1, p_2)$ is:

Proof of Unbiasedness:

To prove $E(\hat{Y}_M) = Y$, evaluate the expectation over all possible unordered sample pairs $s = \{i, j\}$:

Expanding the double sum:

Since $\sum_{j \ne i}^N p_j = 1 - p_i$:

This rigorous proof demonstrates that Murthy's estimator is strictly design-unbiased.

The Exact Variance of Murthy's Estimator:

Murthy derived the closed-form population variance:

Notice the structure of this variance:

1. Like the Sen–Yates–Grundy variance, it depends on differences of ratios $(y_i / p_i - y_j / p_j)^2$. If $y_i$ is perfectly proportional to $p_i$, the variance vanishes identically!
2. The weight coefficient $p_i p_j \left( \frac{1 - p_i - p_j}{2 - p_i - p_j} \right)$ is strictly positive for all $p_i + p_j < 1$.

Murthy's Strictly Non-Negative Sample Variance Estimator:

From a realized sample pair $(y_1, y_2)$ with size probabilities $(p_1, p_2)$, the unbiased sample estimator of variance is:

Because $(1 - p_1) > 0$, $(1 - p_2) > 0$, and $(1 - p_1 - p_2) > 0$ under any standard sampling frame, Murthy's variance estimator is strictly non-negative for every single sample that can ever be drawn. It requires no matrix inversion and no computation of second-order inclusion probabilities, making it an extraordinarily reliable tool in empirical survey research.

When survey operations require complete flexibility or sequential data ingestion across distributed databases, statisticians deploy Bernoulli-type unequal probability mechanisms.

1. Poisson Sampling

In Poisson Sampling (Hajek, 1964), every unit $i \in U$ is subjected to an independent Bernoulli trial with predefined inclusion probability $\pi_i \in (0, 1]$. The sample inclusion indicators $I_1, I_2, \dots, I_N$ are mutually independent random variables:

The sample size $n_s = \sum_{i=1}^N I_i$ is a random variable with:

The Horvitz–Thompson estimator of the total is:

Because the indicators are mutually independent, $\text{Cov}(I_i, I_j) = 0$ for all $i \ne j$. Therefore, the exact sampling variance has zero cross-product terms:

The strictly unbiased sample variance estimator is:

Drawback of Poisson Sampling:

Because the sample size $n_s$ fluctuates randomly around $n$, the variance is inflated by the random sample size variation. To stabilize the sample size, methodologists condition on $\sum I_i = n$ (Conditional Poisson Sampling).

2. Pareto PPS Sampling (Rosén, 1997)

To achieve fixed sample sizes with unequal probabilities without rejection loops, Bengt Rosén (1997) formulated Pareto PPS Sampling, which is widely deployed by Statistics Sweden and Eurostat:

1. For each unit $i \in U$, draw an independent uniform random variable:

1. Compute the Pareto ranking key:

1. Sort the population by $Q_i$ in ascending order:

1. The sample consists of the units with the $n$ smallest values of $Q_i$.

Pareto sampling guarantees an exact fixed sample size $n$, achieves first-order inclusion probabilities that are asymptotically exact ($\pi_i \approx n p_i$), and exhibits outstanding variance stability.

The pinnacle of modern survey design is the integration of PPS sampling into self-weighting two-stage designs. This elegant architecture is deployed by virtually all major international household survey programs, including the Demographic and Health Surveys (DHS), the Multiple Indicator Cluster Surveys (MICS), and national living standards surveys.

The Self-Weighting Architecture

Consider a two-stage national household survey:

• Primary Sampling Units (PSUs) are census enumeration areas containing varying numbers of households $M_i$.
• Secondary Sampling Units (SSUs) are individual households.

If PSUs were selected with equal probability, large enumeration areas would have the same chance of selection as small enumeration areas, requiring complex weighting adjustments at the household level that inflate variance.

Instead, the statistician implements a two-stage design with the following protocol:

1. Stage 1: Select $n$ PSUs with Probability Proportional to Estimated Size (PPS):

where $M_0 = \sum_{i=1}^N M_i$ is the total number of households in the country.

2. Stage 2: Within each selected PSU $i$, select a fixed number of households $m$ (typically $m = 20$ or $25$) using Simple Random Sampling Without Replacement (SRSWOR):

Proof of the Self-Weighting Property

The unconditional overall inclusion probability for any household $j$ residing in PSU $i$ is:

Notice what happens algebraically: the cluster size $M_i$ in the numerator of the first stage cancels out identically with the cluster size $M_i$ in the denominator of the second stage!

The overall design weight for every household in the entire country is:

Why This Design is Universally Revered:

1. Self-Weighting Elegance: Because every household has the exact same design weight, the unweighted sample average $\bar{y} = \frac{1}{nm} \sum \sum y_{ij}$ is a strictly design-unbiased estimator of the national population mean. No design weights are required to compute point estimates.

2. Equal Fieldwork Workload: Because exactly $m$ households are interviewed in every selected cluster, field interview teams have an identical, predictable workload in every village or neighborhood, optimizing team travel and logistical planning.

3. Weight Variance Elimination: Kish's weight design effect $\text{Deff}_w = 1 + L = 1.00$, meaning zero variance is lost due to unequal weights.

By combining unequal probability selection at the first stage with equal sample takes at the second stage, unequal probability sampling achieves the operational dream: equal probability of selection for ultimate elements coupled with balanced field operations. In the next unit, we examine how auxiliary information can be harnessed not only during sample selection, but post-data-collection through Ratio, Regression, and Calibration Estimation.