---
title: 'Pivotal Sampling: Design & Analysis'
url: https://www.emergentmind.com/topics/pivotal-sampling
type: topic
---

# Pivotal Sampling: Design & Analysis

Pivotal sampling is a family of fixed-size unequal-probability sampling designs without replacement in which inclusion-probability mass is repeatedly transferred within pairs until all probabilities are rounded to \(0\) or \(1\), while prescribed first-order inclusion probabilities are preserved. In survey sampling it appears as a particular case of balanced sampling, and in its ordered form it can be represented through microstrata, cross-border units, or an equivalent clustered population. The modern literature establishes several core facts: ordered pivotal sampling has exact marginals, admits explicit second-order inclusion probabilities through an equivalence with Deville’s systematic sampling, dominates multinomial with-replacement sampling in variance, yields weak consistency and martingale central-limit theorems for the Horvitz–Thompson estimator, and supports conservative variance estimation when unbiased HT variance estimation is unavailable. Subsequent work developed the local pivotal method for well-spread sampling in discrete and continuous spaces and adapted pivotal ideas to leverage-score-based active learning [1211.5442; 1609.02688; 1510.08895; 2305.02446; 2310.04966].

## 1. Finite-population formulation and place within balanced sampling

In the standard finite-population formulation, one considers a population \(U=\{1,2,\dots,N\}\) together with prescribed first-order inclusion probabilities \(\pi=(\pi_1,\dots,\pi_N)\) such that \(0<\pi_i<1\) and \(\sum_{i=1}^N \pi_i = n\), where \(n\) is the fixed sample size. Balanced-sampling methods use auxiliary information at the design stage to reduce the variance of Horvitz–Thompson estimation; in its simplest form, when only the fixed-size constraint is imposed, balanced sampling reduces to pivotal sampling. When the units are arranged in a fixed order, for example after sorting by an auxiliary variable, the resulting design is ordered pivotal sampling (OPS), which combines features of systematic sampling and the cube method [1211.5442].

For a study variable \(y\), the Horvitz–Thompson estimator of the population mean in the discrete setting is
\[
\hat\mu_{HT}=\frac1N\sum_{i\in S}\frac{y_i}{\pi_i},
\]
and its variance has the standard form
\[
\Var(\hat\mu_{HT})=\frac1{N^2}\sum_{i=1}^N\sum_{j=1}^N(\pi_{ij}-\pi_i\pi_j)\frac{y_i}{\pi_i}\frac{y_j}{\pi_j},
\]
where \(\pi_{ij}=\Pr\{i\in S,\;j\in S\}\). This representation makes the role of dependence explicit: negative values of \(\pi_{ij}-\pi_i\pi_j\) reduce variance whenever nearby or similar units have positively correlated responses [2305.02446].

A recurrent theme across the literature is that pivotal sampling is not merely a device for satisfying unequal marginal probabilities. It is also a mechanism for inducing controlled dependence: in ordered versions through microstratum structure, in local versions through nearest-neighbor contests, and in tree-based versions through hierarchical pairing. This suggests that pivotal sampling is best understood as a design principle for coupling exact marginals with repulsion among selected units.

## 2. Pairwise updating and ordered construction

The elementary pivotal step acts on a pair \((i,j)\) with current probabilities \((\pi_i,\pi_j)\). If \(\pi_i+\pi_j<1\), the update is
\[
(\pi_i',\pi_j')=
\begin{cases}
(0,\;\pi_i+\pi_j) & \text{w.p. }\dfrac{\pi_j}{\pi_i+\pi_j},\\[6pt]
(\pi_i+\pi_j,\;0) & \text{w.p. }\dfrac{\pi_i}{\pi_i+\pi_j},
\end{cases}
\]
whereas if \(\pi_i+\pi_j\ge 1\),
\[
(\pi_i',\pi_j')=
\begin{cases}
(1,\,\pi_i+\pi_j-1) & \text{w.p. }\dfrac{1-\pi_j}{2-\pi_i-\pi_j},\\[6pt]
(\pi_i+\pi_j-1,\,1) & \text{w.p. }\dfrac{1-\pi_i}{2-\pi_i-\pi_j}.
\end{cases}
\]
Each update forces one unit toward exclusion or inclusion while preserving the pair’s total inclusion mass. Repetition of this rule eventually leaves all probabilities in \(\{0,1\}\) and yields a fixed-size sample [2305.02446].

In OPS, the pairing is not arbitrary. Define cumulative sums \(V_k=\pi_1+\cdots+\pi_k\) with \(V_0=0\). For each integer \(i=1,\dots,n-1\), there is a unique cross-border index \(k_i\) such that \(V_{k_i-1}\le i < V_{k_i}\), and one sets
\[
a_i=i-V_{k_i-1}, \qquad b_i=V_{k_i}-i.
\]
The units between consecutive cross-borders form microstrata \(U_i=\{k_{i-1}+1,\dots,k_i\}\). OPS then proceeds sequentially down the ordered list: a residual “jumper” enters a microstratum, accumulates or contests probability mass against the next cross-border unit, one unit is declared selected, and the remainder is carried forward. Exactly one sampled unit emerges from each microstratum, and exactly \(n\) units are rounded to \(1\) [1211.5442; 1510.08895].

An equivalent construction uses a clustered population \(U^c\) of size \(2n-1\). Each cross-border unit becomes a singleton cluster, while the non-cross-border units between successive cross-borders are aggregated into odd-indexed clusters with cluster probabilities \(\phi=(\phi_1,\dots,\phi_{2n-1})\). Ordered pivotal sampling is then applied to the clustered population, after which each selected cluster is “broken” back into one original unit using probabilities proportional to the original \(\pi\)-mass [1609.02688].

## 3. Design characterization and inclusion probabilities

OPS exactly respects the prescribed first-order inclusion probabilities:
\[
\Pr(k\in S_{op})=\pi_k.
\]
A more delicate question is the form of \(\pi_{kl}\). For OPS this question becomes tractable because OPS with parameter \(\pi\) and Deville’s systematic sampling with parameter \(\pi\) induce the same sampling design. This equivalence is established by matching the two procedures as two-stage designs on the same cluster population and showing that their step-by-step conditional selection probabilities coincide [1211.5442].

Through that characterization, explicit second-order inclusion probabilities become available. Let \(k\in U_i\) and \(l\in U_j\) with \(i<j\), and define
\[
c_\ell=\frac{a_\ell b_\ell}{(1-a_\ell)(1-b_\ell)}, \qquad
c(i,j)=\prod_{\ell=i}^{j-1} c_\ell, \qquad c(i,i)=1.
\]
Then, for example, if \(k\) and \(l\) are both non-cross-border units in the same microstratum \(U_i\), \(\pi_{kl}=0\). If \(k\in U_i\) and \(l\in U_j\) are non-cross-border with \(i<j\), then
\[
\pi_{kl}=\pi_k\pi_l[1-c(i,j)].
\]
Additional cases cover one or two cross-border units and involve the residual quantities \(b_{i-1}\) and \(b_{j-1}\) [1211.5442].

The cluster representation gives a complementary structural result. Both OPS and multinomial sampling admit the same two-stage decomposition on \(U^c\): first one selects \(n\) clusters from \(U^c\), either by ordered pivotal sampling or by multinomial sampling, and then, independently for each selected cluster \(u_i\), one draws exactly one original unit \(k\in u_i\) with probability \(\phi_{i,k}=\pi_k/\phi_i\). Writing \(Y_i=\sum_{k\in u_i} y_k\) for the cluster total, both induced estimators on the clustered population take the form
\[
\hat t=\sum_{i\in S^c} \frac{Y_i}{\phi_i},
\]
so their difference is entirely a first-stage design effect [1609.02688].

A frequent source of confusion concerns the relation between pivotal and systematic procedures. The precise statement supported by the literature is that **ordered** pivotal sampling is equivalent to **Deville’s** systematic sampling as a sampling design. That does not collapse the broader class of pivotal methods into ordinary systematic sampling, nor does it eliminate the role of the pairwise rounding mechanism.

## 4. Estimation, variance dominance, and asymptotic theory

Under without-replacement pivotal sampling on the original population, the usual Horvitz–Thompson estimator of the total is
\[
\hat t_{HT}=\sum_{k\in S}\frac{y_k}{\pi_k}.
\]
Via the cluster coupling, one obtains \(\hat t_{HT}=\hat t_{ops}\). The central comparison in the variance literature is with multinomial with-replacement sampling. Using the two-stage decomposition,
\[
V(\hat t_{ops})=
V\Bigl(\sum_{i\in S^c}\frac{Y_i}{\phi_i}\Bigr)+
E\Bigl[\sum_{i\in S^c}\sum_{k\in u_i}\Bigl(\check y_k-\frac{Y_i}{\phi_i}\Bigr)^2\Bigr],
\]
\[
V(\hat t_{ms})=
V_{ms}\Bigl(\sum_{i\in S^c}\frac{Y_i}{\phi_i}\Bigr)+
E\Bigl[\sum_{i\in S^c}\sum_{k\in u_i}\Bigl(\check y_k-\frac{Y_i}{\phi_i}\Bigr)^2\Bigr],
\]
where \(\check y_k=y_k/\pi_k\). Since the second terms coincide, the variance difference is entirely due to the first stage. Theorem 1 then states that, for every \(\pi\) (equivalently \(\phi\)),
\[
V(\hat t_{ops}) \le V(\hat t_{ms}),
\]
so any implementation of pivotal sampling is more efficient than multinomial sampling [1609.02688].

This variance dominance has direct consequences for consistency. If there exists a constant \(A\) such that
\[
N^{-1}\sum_{k\in U}\pi_k\Bigl(\frac{y_k}{\pi_k}\Bigr)^2 \le A<\infty,
\]
then
\[
E\bigl[N^{-1}(\hat t_{HT}-Y)^2\bigr]=O(n^{-1}),
\]
and therefore \(\hat t_{HT}\to Y\) in probability as \(n\to\infty\). The argument combines the pivotal-versus-multinomial inequality with the classical \(O(n^{-1})\) variance rate for multinomial sampling under mild moment conditions [1609.02688].

Variance estimation is more subtle because many second-order inclusion probabilities under pivotal sampling are zero, so the usual unbiased HT variance estimator is generally undefined. A conservative alternative is the Hansen–Hurvitz variance estimator computed as if the design were with replacement:
\[
v_{HH}(\hat t_{HT})
= \frac{n-1}{n}\sum_{k\in S}\Bigl(\check y_k-\frac1n\sum_{\ell\in S}\check y_\ell\Bigr)^2.
\]
Its expectation satisfies
\[
E[v_{HH}(\hat t_{ops})]-V(\hat t_{ops})
=
\bigl[V(\hat t_{ms})-V(\hat t_{ops})\bigr]\ge 0,
\]
hence \(v_{HH}\) is conservative for \(\Var(\hat t_{HT})\). This resolves a practical difficulty created by vanishing \(\pi_{kl}\) values [1609.02688].

A separate asymptotic line of work proves central-limit theorems for OPS. Under the design-based approach, one writes
\[
\hat Y_{HT}-T=\xi_1+\cdots+\xi_n,
\]
where \(\{\xi_i\}\) is a martingale-difference array adapted to the filtration generated by the sequential microstratum-by-microstratum construction. Under conditions \((H1)\)–\((H3)\), including \(\max_i \pi_i\le f<1\), bounded fourth moments, and nonvanishing within-microstratum dispersion, one obtains
\[
(\hat Y_{HT}-T)/\sqrt{\Var_p(\hat Y_{HT})}\Rightarrow N(0,1).
\]
Under a model-assisted framework with
\[
y_k=\beta\,\pi_k+\varepsilon_k,
\]
allowing nonzero \(\Cov(\varepsilon_k,\varepsilon_\ell)\) for nearby units and suitable decay conditions, one similarly obtains
\[
(\hat Y_{HT}-T)/\sqrt{\Var_{m,p}(\hat Y_{HT})}\Rightarrow N(0,1).
\]
These results are especially relevant in spatial sampling, where correlation among neighboring units is expected rather than exceptional [1510.08895].

## 5. Ordering effects, spatial balance, and the local pivotal method

The effect of ordering is particularly transparent when \(\pi_k=n/N\) and \(N=np\). In that case the microstrata are non-overlapping sets of size \(p\), and OPS reduces to stratified simple random sampling with one draw per stratum. The corresponding variance formulas are
\[
V_{srs}(\hat t_{y\pi})=N^2(1-f)/n\times S_y^2,
\]
\[
V_{ops}(\hat t_{y\pi})=N^2(1-f)/n\times (1/n)\sum_{i=1}^n S_{y,i}^2,
\]
\[
V_{sys}(\hat t_{y\pi})=N^2(1-f)/n\times (1/n)S_Y^2.
\]
The worst-case design effects satisfy
\[
DMAX(ops)=\frac{N-1}{N-n}, \qquad DMAX(sys)=\frac{n(N-1)}{N-n}.
\]
Thus OPS can never do much worse than SRS, whereas systematic sampling can be poor when the order is unfavorable. In the same equal-\(\pi\) setting, the dispersion of the eigenvalues of the variance-covariance matrix \(\Delta=[\pi_{kl}-\pi_k\pi_l]\) satisfies
\[
\delta(srs)<\delta(ops)<\delta(sys),
\]
so OPS is more “egalitarian” than systematic sampling, though less so than SRS [1211.5442].

These ordering results connect directly to spatial balance. In longitudinal rotation designs the natural ordering is time; in spatial surveys one often orders units along a space-filling curve. OPS then produces well-spread samples and can exploit local positive correlation in \(y\) without incurring the extreme pathologies associated with pure systematic selection [1510.08895].

The local pivotal method (LPM) modifies the pairing rule so that contests occur between nearby points in an auxiliary space rather than along a fixed one-dimensional order. To sample from an arbitrary continuous distribution \(\widetilde U\subset \mathbb{R}^d\), one first draws a large i.i.d. master sample \(x_1,\dots,x_N\sim \widetilde U\), sets \(\pi_i=n/N\) in the equal-probability case, and then repeatedly selects a point \(i\) with \(0<\pi_i<1\), finds its nearest neighbor \(j\), and applies the same two-point pivotal update. LPM2 is a minor variant that only updates pairs that are each other’s nearest neighbor. The resulting subset of size \(n\) is “thin” and well spread; the repeated local pairing generates an automatic, data-driven stratification at a scale tuned by \(n\), and the sample is as “repulsive” as possible under the constraint \(\sum \pi_i=n\) [2305.02446].

For smooth integrands \(y:\mathbb{R}^d\to\mathbb{R}\), the continuous-population analysis gives
\[
\Var(\hat\mu_{LPM})=O\bigl(N^{-1}+n^{-1-2/d}\bigr)
\]
under regularity such as Lipschitz or bounded variation. The \(n^{-1-2/d}\) term is the stratification-type improvement over the \(n^{-1}\) i.i.d. rate. This suggests that local pivotal sampling operationalizes automatic stratification in spaces where manual stratification would be cumbersome or ill-defined [2305.02446].

## 6. Empirical evidence, implementations, and recent extensions

The variance comparison between OPS and multinomial sampling has been checked numerically on clustered populations. For \(n=3\), so \(|U^c|=5\), the probabilities \(\phi_i\) were varied over the grid \(\{0.05,0.10,\dots,0.95\}\) subject to \(\sum \phi_i=3\), yielding \(24\,396\) cases. For each case the second-order inclusion matrix \(B=(\phi_{ij}/(\phi_i\phi_j))_{i\ne j}\) was formed, its second largest eigenvalue \(\lambda\) was computed, and the supremum of the variance ratio \(V_{ops}/V_{ms}\) was identified with \(\lambda\). The observed range was \(\lambda\in[0.625,0.991]\). For \(n=5\), with \(|U^c|=9\) and a \(0.10\)-spaced grid, \(31\,998\) cases gave \(\lambda\in[0.666,0.975]\). In both settings the ratios remained at or below \(1\), in accord with the theoretical inequality [1609.02688].

In the continuous and Monte Carlo setting, the LPM literature reports concrete precision gains. For estimating \(\mu=\int_0^1 x\,dx=0.5\) with \(N=10^4\) and \(n=100\), the standard deviation is approximately \(0.028\) for the i.i.d. estimator and approximately \(0.004\) for LPM2. In European call option pricing with \(P_{EC}\approx 3.886\) and \(N=10^4\), the estimated standard deviation at \(n=100\) is approximately \(0.90\) under i.i.d. sampling and approximately \(0.31\) under LPM2. In nonlinear dynamical-system stability estimation, LPM2 reduces sampling noise by a factor \(3\)–\(4\) in a rainforest bistability example at \(n=50\), and in a 4D forced Jeffcott rotor example cuts the standard deviation from approximately \(0.02\) to approximately \(0.007\) at \(n=500\) [2305.02446].

The same work provides implementation hooks. In R, the CRAN package **BalancedSampling** supplies `lpm1(prob,X)` and `lpm2(prob,X)`, where `prob` is the vector of \(\pi_i\) and `X` is the auxiliary matrix. In MATLAB, the function
```matlab
[s,svar] = lpm2(prob,X,distfcn,ns)
```
supports custom distances and local-mean variance estimation through nearest neighbors [2305.02446].

A distinct recent extension places pivotal sampling inside active regression. Given a matrix \(A\in\mathbb{R}^{n\times d}\), leverage scores \(\tau_i=\|u_i\|_2^2=a_i^\top(A^\top A)^{-1}a_i\), and a target sample size \(k\), one sets raw marginals proportional to leverage scores, ceilings them at \(1\), and then applies a binary-tree-based pivotal algorithm in which nearby points compete early in the tree. The resulting design samples exactly \(k\) points, preserves the marginals, and has a negatively dependent joint law. Because pivotal sampling on a fixed tree is strongly Rayleigh, it satisfies one-sided \(\ell_\infty\)-independence with parameter \(D=O(1)\). This leads to a general active-regression guarantee with
\[
k=O(d\log d+d/\varepsilon),
\]
and, for polynomial regression, an improved bound of
\[
k=O(d/\varepsilon).
\]
In experiments on surrogate-modeling problems from parametric PDEs and uncertainty quantification, the pivotal method reduced the number of labels needed to reach a fixed error by up to \(50\%\) relative to Bernoulli leverage-score sampling [2310.04966].

Taken together, these results show that pivotal sampling has developed from a survey-sampling algorithm into a broader methodology for dependent subsampling with exact marginals, negative dependence, and spatial coverage. The survey-theoretic core remains the Horvitz–Thompson framework and the microstratum or clustered-population analysis; the newer developments demonstrate that the same pairwise-rounding principle is effective in continuous Monte Carlo designs and in leverage-score-guided learning problems.

Source: https://www.emergentmind.com/topics/pivotal-sampling