Survey Sampling Algorithms

Fast implementations of survey sampling algorithms for single-stage probability sampling from finite populations, written in C. Provides equal probability methods (simple random sampling, systematic, Bernoulli), unequal probability methods (conditional Poisson / maximum entropy, Sampford, Brewer, systematic PPS, Pareto, sequential Poisson, Poisson, Chromy's minimum replacement, multinomial), balanced sampling via the cube method, and spatially balanced sampling via the local pivotal method and spatially correlated Poisson sampling. All sampling functions return design objects carrying sample indices, inclusion probabilities, and design metadata. Generics compute joint inclusion probabilities, pairwise expectations, and sampling covariances for variance estimation. Mostly based on algorithms from Tillé (2006, ).

sondage

Fast survey sampling algorithms for R. Sampling functions return design objects with generics for extracting inclusion probabilities, joint inclusion probabilities, and variance estimation quantities.

For without-replacement designs, the stored pik vector is the design-defining target inclusion probability vector. For methods with exact first-order guarantees, this equals the true first-order inclusion probabilities. For order-sampling methods such as sps and pareto, the stored vector remains the target pik, while the true finite-population first-order inclusion probabilities are only approximately equal to that target.

Installation

# From GitLab
pak::pkg_install("gitlab::dickoa/sondage")

Usage

library(sondage)

# Use built-in US state data
data(state)
states <- as.data.frame(state.x77)

# Compute inclusion probabilities from population size
pik <- inclusion_prob(states$Population, n = 10)

# Draw a sample (Conditional Poisson Sampling)
s <- unequal_prob_wor(pik, method = "cps")
states[s$sample, ]
#>              Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
#> California        21198   5114        1.1    71.71   10.3    62.6    20 156361
#> Georgia            4931   4091        2.0    68.54   13.9    40.6    60  58073
#> Michigan           9111   4751        0.9    70.63   11.1    52.8   125  56817
#> Mississippi        2341   3098        2.4    68.09   12.5    41.0    50  47296
#> Missouri           4767   4254        0.8    70.69    9.3    48.8   108  68995
#> Nebraska           1544   4508        0.6    72.60    2.9    59.3   139  76483
#> New York          18076   4903        1.4    70.55   10.9    52.7    82  47831
#> Pennsylvania      11860   4449        1.0    70.43    6.1    50.2   126  44966
#> Washington         3559   4864        0.6    71.72    4.3    63.5    32  66570
#> Wisconsin          4589   4468        0.7    72.48    3.0    54.5   149  54464

# Joint inclusion probabilities for variance estimation
pikl  <- joint_inclusion_prob(s)
delta <- sampling_cov(s)                # pi_ij - pi_i * pi_j
chk   <- sampling_cov(s, weighted = TRUE) # 1 - pi_i * pi_j / pi_ij

# Equal probability sampling
s <- equal_prob_wor(nrow(states), 10)
states[s$sample, ]
#>                Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
#> Minnesota            3921   4675        0.6    72.96    2.3    57.6   160  79289
#> Colorado             2541   4884        0.7    72.06    6.8    63.9   166 103766
#> South Carolina       2816   3635        2.3    67.96   11.6    37.8    65  30225
#> Utah                 1203   4022        0.6    72.90    4.5    67.3   137  82096
#> Missouri             4767   4254        0.8    70.69    9.3    48.8   108  68995
#> Wisconsin            4589   4468        0.7    72.48    3.0    54.5   149  54464
#> Rhode Island          931   4558        1.3    71.90    2.4    46.4   127   1049
#> Tennessee            4173   3821        1.7    70.11   11.0    41.8    70  41328
#> Vermont               472   3907        0.6    71.64    5.5    57.1   168   9267
#> Mississippi          2341   3098        2.4    68.09   12.5    41.0    50  47296

# PPS with minimum replacement (Chromy)
hits <- expected_hits(states$Population, n = 10)
s <- unequal_prob_wr(hits, method = "chromy")

# Balanced sampling (cube method)
pik <- inclusion_prob(states$Population, n = 10)
x <- matrix(states$Income)
s_bal <- balanced_wor(pik, aux = x)
s_bal
#> Balanced WOR [cube] (n=10, N=50): 5 10 13 14 18 24 25 32 38 44

# Controlled selection e.g to keep the number of sampled states per region
# within the integers adjacent to its expectation
B <- sapply(levels(state.region), function(g) as.double(state.region == g))
S <- colSums(B * pik)
s_ctrl <- balanced_wor(
  pik,
  bounds = list(B = B, lower = floor(S), upper = ceiling(S))
)
table(state.region[s_ctrl$sample])
#> 
#>     Northeast         South North Central          West 
#>             2             3             3             2

# Batch sampling for simulations (design object with matrix $sample)
sim <- unequal_prob_wor(pik, method = "cps", nrep = 1000)
dim(sim$sample)   # 10 x 1000
#> [1]   10 1000
inclusion_prob(sim) # generics still work
#>  [1] 0.17026107 0.01719095 0.10418188 0.09937783 0.99839394 0.11967728 0.14600534 0.02727003 0.38983426
#> [10] 0.23224269 0.04088150 0.03829108 0.52736187 0.25023432 0.13474880 0.10738457 0.15952261 0.17925688
#> [19] 0.04983021 0.19414000 0.27383066 0.42911441 0.18467321 0.11025758 0.22451854 0.03513548 0.07272008
#> [28] 0.02778811 0.03824398 0.34537328 0.05388068 0.85135243 0.25626292 0.03000174 0.50560237 0.12787242
#> [37] 0.10757297 0.55858818 0.04384870 0.13262937 0.03207408 0.19654203 0.57634431 0.05665949 0.02223049
#> [46] 0.23459761 0.16762355 0.08473020 0.21613500 0.01770903

Sampling functions

Equal probability without replacement (equal_prob_wor):

• equal_prob_wor(N, n, method = "srs") - Simple random sampling
• equal_prob_wor(N, n, method = "systematic") - Systematic sampling
• equal_prob_wor(N, n, method = "bernoulli") - Bernoulli sampling (random size)

Equal probability with replacement (equal_prob_wr):

• equal_prob_wr(N, n, method = "srs") - Simple random sampling with replacement

Unequal probability without replacement (unequal_prob_wor):

• unequal_prob_wor(pik, method = "cps") - Conditional Poisson / maximum entropy
• unequal_prob_wor(pik, method = "sampford") - Exact Sampford PPS with exact joint probabilities
• unequal_prob_wor(pik, method = "brewer") - Brewer’s method
• unequal_prob_wor(pik, method = "systematic") - Systematic PPS
• unequal_prob_wor(pik, method = "poisson") - Poisson sampling (random size)
• unequal_prob_wor(pik, method = "sps") - Sequential Poisson sampling (order sampling)
• unequal_prob_wor(pik, method = "pareto") - Pareto sampling (order sampling)

Unequal probability with replacement (unequal_prob_wr):

• unequal_prob_wr(hits, method = "chromy") - PPS with minimum replacement
• unequal_prob_wr(hits, method = "multinomial") - Multinomial PPS

Balanced sampling without replacement (balanced_wor):

• balanced_wor(pik, aux, method = "cube") - Cube method (Deville & Tillé, 2004)
• balanced_wor(pik, aux, strata, method = "cube") - Stratified cube (Chauvet, 2009)
• balanced_wor(pik, bounds = list(B, lower, upper)) - Cube with inequality constraints (Tripet & Tillé, 2026) providing controlled matrix rounding, minimum group sizes
• balanced_wor(pik, spread, method = "lpm2") - Spatially balanced, well-spread sampling with the local pivotal method 2 (Grafström, Lundström & Schelin, 2012)
• balanced_wor(pik, spread, method = "scps") - Spatially correlated Poisson sampling with Grafström’s (2012) maximal-weight strategy

Design queries

• inclusion_prob(x, n) - Compute inclusion probabilities from size measures
• inclusion_prob(s) - Extract the stored design-defining pik vector from a WOR design
• expected_hits(x, n) - Compute expected hits from size measures
• expected_hits(s) - Extract expected hits from a WR design
• joint_inclusion_prob(s) - Joint inclusion probabilities (WOR)
• joint_inclusion_prob(s, sampled_only = TRUE) - n x n submatrix for sampled units only (scales to large N)
• joint_expected_hits(s) - Pairwise expectations E(n_i n_j) (WR)
• joint_expected_hits(s, sampled_only = TRUE) - Submatrix for selected units only
• sampling_cov(s) - Sampling covariance matrix
• sampling_cov(s, weighted = TRUE) - Check quantities for SYG variance estimator
• sampling_cov(s, sampled_only = TRUE) - Covariance for sampled units only

Joint and covariance matrices preserve population-unit names. When the probability vector is unnamed, sampled-only matrices use population indices as row and column names so their units remain identifiable.

Method comparison

†For WOR methods, design marginals are first-order inclusion probabilities \pi_k. For WR methods, design marginals are expected hits E(N_k).

*For sps and pareto, inclusion_prob(s) returns the stored design-defining target vector; the true finite-population first-order inclusion probabilities are only approximately equal to that target, with the discrepancy vanishing asymptotically.

**For sps and pareto, the high-entropy approximation is built from the stored target pik vector. HE = high-entropy approximation.

Choosing a method

• Use cps when exactness matters more than speed: it is the maximum-entropy fixed-size unequal-probability design, with exact first and second-order inclusion probabilities.
• Use brewer when you want exact first-order inclusion probabilities with lower computational cost than cps, and can work with approximate second-order quantities in sondage.
• Use sampford for the SAS-familiar fixed-size PPS design with exact first- and second-order inclusion probabilities. Its C kernel uses complement sampling and a bounded-rejection/non-rejective hybrid.
• Use systematic when very fast sampling and ordering or implicit stratification are central, and structural zeros in some joint inclusion probabilities are acceptable.
• Use poisson when a random sample size is acceptable and independent selection is desirable.
• Use sps or pareto as fast high-entropy order-sampling alternatives when approximate first- and second-order quantities are acceptable.
• Use cube when balancing on auxiliary variables is more important than exact second-order inclusion probabilities.
• Use lpm2 or scps when the study variable is spatially structured and a well-spread sample matters more than joint inclusion probabilities. LPM2 makes local pairwise competitions and SCPS distributes each decision over the nearest feasible neighbours. Variance is then usually estimated with local-neighbourhood estimators.

Custom methods

register_method() lets you plug any unequal-probability, balanced, or spatially balanced sampling algorithm into sondage’s dispatchers and generics. Methods registered with type = "wor" or type = "wr" dispatch through unequal_prob_wor() / unequal_prob_wr(), and balanced methods (type = "balanced") through balanced_wor(), where they declare which design inputs they use (supports_aux, supports_strata, supports_spread):

Here type = "balanced" names the dispatcher family, which includes spatially balanced designs. Spread-only methods such as LPM2 and SCPS do not exactly balance auxiliary totals; their supports_aux = FALSE metadata makes that distinction enforceable rather than silently ignoring aux.

# A simple randomized pivotal sampler written in R
random_pivotal_sample <- function(pik, n = NULL, prn = NULL, ...) {
  p <- pik
  tol <- 1e-06
  active <- which(p > tol & p < 1 - tol)
  while (length(active) >= 2L) {
    ij <- sample(active, 2L)
    i <- ij[1L]
    j <- ij[2L]
    total <- p[i] + p[j]
    if (total < 1) {
      if (runif(1) < p[i] / total) p[c(i, j)] <- c(total, 0)
      else p[c(i, j)] <- c(0, total)
    } else {
      if (runif(1) < (1 - p[j]) / (2 - total)) {
        p[c(i, j)] <- c(1, total - 1)
      } else {
        p[c(i, j)] <- c(total - 1, 1)
      }
    }
    active <- which(p > tol & p < 1 - tol)
  }
  sort(which(p > 0.5))
}

register_method("random_pivotal", type = "wor", sample_fn = random_pivotal_sample)

pik <- inclusion_prob(1:8, n = 3)
s <- unequal_prob_wor(pik, method = "random_pivotal")
s
#> Unequal prob WOR [random_pivotal] (n=3, N=8): 6 7 8

unregister_method("random_pivotal")

See vignette("custom-methods") for more examples, including a custom balanced method with stratification support, a spatially balanced method using spread, and how to provide a joint_fn for variance estimation.

Why not sampling?

The sampling package (Tillé and Matei) is the reference toolkit for survey sampling in R, and it is more comprehensive than sondage. Many of the algorithms here follow the methods it established, and sondage would not exist without it.

What sondage adds is speed. The sampling algorithms are written in C, so they usually scale better to large populations. Every sampling function also returns a design object with S3 generics for inclusion probabilities, joint inclusion probabilities, and variance quantities, so the results are easy to carry into downstream work.

The two packages are complementary rather than competing. With register_method() you can plug any unequal probability algorithm from sampling into the sondage dispatchers and generics, so they can be used together.

References

Brewer, K.R.W. and Donadio, M.E. (2003). The High Entropy Variance of the Horvitz-Thompson Estimator. Survey Methodology, 29(2), 189-196.

Chauvet, G. (2009). Stratified balanced sampling. Survey Methodology, 35, 115-119.

Chromy, J.R. (1979). Sequential sample selection methods. Proceedings of the Survey Research Methods Section, American Statistical Association, 401-406.

Chromy, J.R. (2009). Some generalizations of the Horvitz-Thompson estimator. Proceedings of the Survey Research Methods Section, American Statistical Association.

Deville, J.C. and Tillé, Y. (2004). Efficient balanced sampling: the cube method. Biometrika, 91(4), 893-912.

Grafström, A. (2012). Spatially correlated Poisson sampling. Journal of Statistical Planning and Inference, 142(1), 139-147.

Grafström, A., Lundström, N.L.P. and Schelin, L. (2012). Spatially balanced sampling through the pivotal method. Biometrics, 68(2), 514-520.

Tripet, A. and Tillé, Y. (2026). Balanced sampling with inequalities: application to category bounding, matrix rounding, and spread sampling. Journal of the American Statistical Association, 121(553), 796-806.

Tillé, Y. (2006). Sampling Algorithms. Springer.