API Reference

The workflow keeps three things apart, and the boundaries are the point: a ranking of entities by importance, a probability design over them, and a realized draw. Only the last can be estimated from.

Study

The user-facing workflow: data, then a plan, then a draw, then an estimate.

Four objects, and the boundaries between them are the point:

Study the data. Counts, features, and the influence decomposition. Ranking an ordering of entities by how much labeling them would help. LabelDesign a probability for every entity. A plan, not yet a draw. LabelSample a realized draw, carrying the probabilities delivered. Estimate coefficients and their uncertainty.

Only a LabelSample reaches an estimator. A Ranking cannot, and that is deliberate: the top entities by influence are an excellent operational priority list, but they are not a probability sample, and no probability attached to them afterward makes them one. Version 2 blurred these – calibrate_weights invented inclusion probabilities from the size of whatever set it was handed – and the result was a workflow the documentation called Horvitz-Thompson while it was neither unbiased nor design-based.

class fewlab.study.Ranking(items, scores, loss, coverage)

Bases: object

An ordering of entities by contribution, for operational use only.

A ranking answers “if I can only look at a few, which ones?”. It does not answer “what is B, and how sure am I?” – for that the selection has to be random, because unbiasedness comes from every entity having had a chance to be chosen. Passing a Ranking where a LabelSample belongs raises rather than silently producing a number that looks like an estimate.

Parameters:

• items (Index)

• scores (Series)

• loss (str)

• coverage (float)

items

Entities in decreasing order of contribution.

Type:

pandas.Index

scores

The contribution scores b_j that produced the order.

Type:

pandas.Series

loss

The estimand the scores were computed for.

Type:

str

coverage

Share of total contribution mass the ranked entities carry.

Type:

float

to_design(budget, **kwargs)

Point at the supported route from a ranking to an estimable design.

Parameters:

• budget (int) – Labeling budget.

• **kwargs (object) – Forwarded to the design constructor.

Returns:

Never returns.

Raises:

NotImplementedError – Always. Build the design from Study.design, which starts from the scores rather than from a truncated list.

Return type:

LabelDesign

class fewlab.study.Study(counts, X, *, ridge=None, ensure_full_rank=True)

Bases: object

Counts, features, and everything derived from them.

Parameters:

• counts (CountMatrix)

• X (DataFrame)

• ridge (float | None)

• ensure_full_rank (bool)

counts

Count matrix, units by labelable entities.

Type:

fewlab.counts.CountMatrix

X

Unit feature matrix.

Type:

pandas.DataFrame

influence

The exact decomposition B = sum_j h_j a_j’.

Type:

fewlab.influence.Influence

diagnostics

Conditioning and degeneracy detail.

Type:

dict

classmethod from_long(frame, X, *, unit, entity, count, entities=None, **kwargs)

Build a study from long-format usage data.

This is the natural entry point for a real panel, which arrives as one row per unit-entity pair rather than as a dense matrix. A browsing panel of 1,134 users by 34,078 domains is 0.27% dense: 309 MB dense, a couple of megabytes as triples.

Parameters:

• frame (DataFrame) – Long-format counts.

• X (DataFrame) – Unit feature matrix indexed by unit identifier.

• unit (str) – Column naming the unit.

• entity (str) – Column naming the labelable entity.

• count (str) – Column holding the count.

• entities (Index | None) – Optional full entity frame, when some entities should belong to the population despite having no recorded activity.

• **kwargs (object) – Forwarded to the constructor.

Returns:

The study.

Return type:

Study

Examples

>>> import pandas as pd
>>> from fewlab import Study
>>> long = pd.DataFrame(
...     {"user": ["a", "a", "b"], "site": ["x", "y", "x"],
...      "visits": [3.0, 1.0, 5.0]}
... )
>>> X = pd.DataFrame({"const": [1.0, 1.0], "z": [0.0, 1.0]},
...                  index=["a", "b"])
>>> Study.from_long(long, X, unit="user", entity="site",
...                 count="visits").influence.n_items
2

scores(*, loss='coefficient', contrast=None, label_norms=1.0)

Compute each entity’s contribution b_j to the design risk.

Parameters:

• loss (str) – “coefficient” or “prediction”. See contribution_scores.

• contrast (ndarray | None) – Optional (p, r) matrix targeting L’B only.

• label_norms (ndarray | Series | float) – Squared label norms a_j’ a_j; 1.0 is exact for a one-hot taxonomy.

Returns:

Contribution scores indexed by entity.

Return type:

Series

rank(budget=None, **score_kwargs)

Order entities by contribution, for operational triage.

Parameters:

• budget (int | None) – How many to return, or None for all.

• **score_kwargs (object) – Forwarded to scores.

Returns:

The ranking. It cannot be estimated from; see design.

Return type:

Ranking

design(budget, *, method='pps', cost=None, **score_kwargs)

Build a probability design for a labeling budget.

Parameters:

• budget (float) – Total labeling budget. With uniform cost this is the expected number of entities labeled.

• method (Literal['pps', 'core_plus_tail', 'srs']) – “pps” for the risk-minimizing square-root rule, “core_plus_tail” for a certainty core with a uniform tail, or “srs” for the uniform benchmark.

• cost (Series | ndarray | float | None) – Per-entity labeling cost. Only “pps” uses it.

• **score_kwargs (object) – Forwarded to scores.

Returns:

The design.

Raises:

ValueError – If method is unrecognized.

Return type:

LabelDesign

design_for_risk(risk_limit, *, cost=None, **score_kwargs)

Buy the least expected cost for a label-sampling risk ceiling.

Use independent draws via draw(design, fixed_size=False). A single contrast uses a variance ceiling (SE squared); the default coefficient loss bounds the sum of coefficient variances. With unit label norms, risk is exact for a one-hot taxonomy and bounds an unknown binary label. The ceiling excludes regression uncertainty and realized-sample error. A contrast design protects only that contrast, not other coefficients.

Parameters:

• risk_limit (float) – Finite, nonnegative Poisson risk ceiling.

• cost (Series | ndarray | float | None) – Positive costs, or None for minimum expected label count.

• **score_kwargs (object) – Forwarded to scores, including loss, contrast, and label_norms. Estimated norms give a planning criterion.

Returns:

The design, with achieved risk and expected cost in diagnostics.

Return type:

LabelDesign

draw(design, *, random_state=None, fixed_size=True)

Realize a design into an actual set of entities to label.

Parameters:

• design (LabelDesign) – The probability design to realize.

• random_state (int | Generator | None) – Seed or generator.

• fixed_size (bool) – Draw exactly budget entities using the cube method, balancing on the influence projections so the sample reproduces their population totals. False draws independent Bernoulli trials, giving a random sample size but the simplest variance.

Returns:

The realized sample, carrying the delivered probabilities.

Raises:

ValueError – If a fixed-size draw is asked for a non-integer budget.

Return type:

LabelSample

calibrate(sample)

Calibrate a realized sample so it reproduces the influence totals.

The estimator weights a sampled entity by 1 / pi_j, which is unbiased but leaves the realized sample’s influence total sum_{j in S} g_j / pi_j scattered around its population value sum_j g_j. Calibration adjusts the weights by the smallest chi-square distance that makes that equality hold exactly, which removes the part of the error attributable to an unlucky draw.

This matters for a Poisson draw and much less for a cube draw, because Study.draw(fixed_size=True) already balances on G at selection time and so arrives calibrated. Use it when you took the simpler sampler and want the balance back.

Returned as shrinkage factors rather than weights so they compose with estimate: the calibrated weight for entity j is alpha_j / pi_j. Note that calibration is asymptotically unbiased, not exactly so – the adjustment depends on the realized sample – which is why the resulting estimate reports itself as no longer strictly unbiased.

Parameters:

sample (LabelSample) – The realized draw to calibrate.

Returns:

Shrinkage factors indexed by the sampled entities, to pass as alpha to estimate.

Return type:

Series

estimate(sample, labels, *, proxy=None, alpha=1.0, variance='poisson')

Estimate B from a realized sample.

Parameters:

• sample (LabelSample) – The realized draw. A Ranking is rejected here on purpose.

• labels (Series | DataFrame) – Observed labels for the sampled entities. Arbitrary q-vectors are accepted, so soft labels – a coder distribution, a probability vector over classes – need no special handling.

• proxy (DataFrame | None) – Optional cheap prediction for every entity, enabling the difference estimator.

• alpha (Series | float) – Shrinkage on the sampled correction; 1.0 is unbiased.

• variance (str) – Which variance estimator to use. The default is exact for a Poisson draw, with no universal bound for other designs. Pass “balanced” after a cube draw to account for the auxiliaries used in balancing.

Returns:

The estimate with both variance components.

Raises:

TypeError – If handed a Ranking or a bare list of entities.

Return type:

Estimate

Influence

Exact influence decomposition for entity-level labels.

The problem this package solves has an exact finite-population identity, and this module computes it. There is no covariance model to approximate and no distributional assumption to make.

Setup. Rows are users, columns are entities that carry an expensive label – a website, a merchant, a name. C_ij counts how often user i used entity j, T_i is the row total, and V = C / T holds the row shares. Entity j carries an unobserved label vector a_j, and the outcome of interest is the user-level share Y = V A: user i’s share of activity on entities of each label class. Regressing that outcome on user features X gives

B = (X’X)^-1 X’ Y = (X’X)^-1 X’ V A = H A, H = (X’X)^-1 X’ V

and therefore, writing h_j for the j-th column of H,

B = sum_j h_j a_j’

exactly, with no error term. Every entity’s entire contribution to the estimand is the rank-one matrix h_j a_j’, so h_j – computable from counts and features alone, before a single label is bought – is what makes one entity worth more to label than another.

class fewlab.influence.Influence(H, G, XtX, v_norm2, items, ridge)

Bases: object

Per-entity influence on the downstream regression.

Parameters:

• H (ndarray)

• G (ndarray)

• XtX (ndarray)

• v_norm2 (ndarray)

• items (list[str])

• ridge (float | None)

H

(p, m) matrix (X’X)^-1 X’ V. Column j is h_j, entity j’s exact contribution direction, so that B = sum_j h_j a_j’.

Type:

numpy.ndarray

G

(p, m) matrix X’ V. Used as the balancing auxiliary by the cube sampler, where preserving sum_j g_j is the calibration target.

Type:

numpy.ndarray

XtX

(p, p) cross-product X’X, including any ridge term applied.

Type:

numpy.ndarray

v_norm2

(m,) squared norms ||v_j||^2 of the share columns. Needed only when the label enters as a covariate; see contribution_scores(loss=”regressor”).

Type:

numpy.ndarray

items

Entity identifiers, in the column order of H and G.

Type:

list[str]

ridge

Ridge value added to X’X, or None if the solve was exact.

Type:

float | None

property n_items: int

Number of labelable entities.

Returns:

Column count of H.

property n_features: int

Number of user features.

Returns:

Row count of H.

fewlab.influence.influence(counts, X, *, ridge=None, ensure_full_rank=True)

Compute the exact influence decomposition B = sum_j h_j a_j’.

Parameters:

• counts (DataFrame | CountMatrix) – Count matrix, users as rows and labelable entities as columns. Either a dense DataFrame or a CountMatrix in coordinate form, which is what a real usage panel should use. Rows whose total is zero are dropped, since they contribute no shares.

• X (DataFrame) – User feature matrix aligned with counts.index.

• ridge (float | None) – Ridge term added to X’X. None selects it automatically when ensure_full_rank is set and X’X is ill-conditioned.

• ensure_full_rank (bool) – Whether to apply the automatic ridge heuristic. Has no effect when ridge is given explicitly.

Returns:

The influence decomposition for these counts and features.

Raises:

ValueError – If X does not cover every unit in counts, or if every unit has a zero row total so no shares exist.

Return type:

Influence

Note

With a non-zero ridge the operator is H = (X’X + rI)^-1 X’ V, so the estimand becomes the ridge coefficient rather than OLS. The identity B = sum_j h_j a_j’ holds either way – it needs only that B is a fixed linear map applied to A.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import influence
>>> rng = np.random.default_rng(0)
>>> counts = pd.DataFrame(rng.poisson(5, (200, 40)))
>>> X = pd.DataFrame(np.c_[np.ones(200), rng.normal(size=(200, 2))])
>>> infl = influence(counts, X)
>>> infl.H.shape
(3, 40)

fewlab.influence.contribution_scores(infl, *, loss='coefficient', contrast=None, label_norms=1.0)

Compute each entity’s contribution b_j to the design risk.

Under Poisson sampling with inclusion probabilities pi_j, the exact risk of the Horvitz-Thompson estimator under the quadratic loss tr[W (Bhat - B) (Bhat - B)’] is

sum_j (1 - pi_j) / pi_j * b_j, b_j = (h_j’ W h_j) * (a_j’ a_j)

so b_j is the entire statistical value of labeling entity j. It is the quantity the optimal design maximizes value against, and this function is where the choice of estimand enters. There is no target-free answer: the entity worth labeling to pin down a coefficient is not always the entity worth labeling to predict a fitted value.

Parameters:

• infl (Influence) – Influence decomposition from influence.

• loss (str) – Which estimand to optimize for. “coefficient” sets W = I and targets the mean squared error of the coefficients themselves. “prediction” sets W = X’X and targets fitted values. “regressor” targets the attenuation of a coefficient on the constructed share, which depends on the exposure orthogonal to the controls rather than aligned with them. Ignored when contrast is given.

• contrast (ndarray | None) – Optional (p, r) matrix L targeting only L’B, giving W = L L’ and b_j = ||L’ h_j||^2 (a_j’ a_j). Use it when a single coefficient or a specific comparison is the estimand.

• label_norms (ndarray | Series | float) – The squared label norms a_j’ a_j. The default of 1.0 is exact for a complete one-hot taxonomy and an upper bound for an unknown binary indicator. Supply an array for known squared label norms; estimated norms give a planning criterion, not a guarantee.

Returns:

The scores b_j, indexed by entity.

Raises:

ValueError – If loss is unrecognized, or contrast or label_norms has a shape inconsistent with infl.

Return type:

Series

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import contribution_scores, influence
>>> rng = np.random.default_rng(0)
>>> counts = pd.DataFrame(rng.poisson(5, (200, 40)))
>>> X = pd.DataFrame(np.c_[np.ones(200), rng.normal(size=(200, 2))])
>>> b = contribution_scores(influence(counts, X))
>>> bool((b > 0).all())
True

fewlab.influence.direction_spectrum(infl)

Report how many distinct directions the influence vectors span.

Normalize every h_j to unit length and take the singular values of the result, scaled so the largest is 1. A spectrum like [1, 0.001, 0, 0] means every entity influences the regression in essentially the same direction, differing only in magnitude.

That case is worth detecting because it collapses part of the design problem: when all h_j are parallel, b_j is proportional to ||h_j||^2 under every W, so loss=”coefficient”, loss=”prediction”, and any contrast all produce the same ranking and the same probabilities. Choosing the estimand stops mattering. Probability-proportional sampling still beats simple random sampling – the magnitudes still differ – but a benchmark run on such a population cannot distinguish one criterion from another, and will report every method as equivalent no matter what it is measuring.

Read it against a noise floor, not a fixed cutoff. Counts are integers, and sampling noise alone injects spurious directions with singular values of roughly 1/sqrt(typical count): on an exactly rank-one population the second value measures 0.17 at a mean count of 0.5, 0.05 at 5, and 0.015 at 100. A population is concentrated in one direction when the second value sits at that floor, and genuinely multi-directional when it stands well above it.

Parameters:

infl (Influence) – Influence decomposition from influence.

Returns:

Singular values in decreasing order, normalized so the first is 1.

Return type:

Series

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import direction_spectrum, influence
>>> rng = np.random.default_rng(0)
>>> activity = np.exp(rng.normal(size=300))
>>> popularity = np.exp(rng.normal(size=50))
>>> counts = pd.DataFrame(rng.poisson(np.outer(activity, popularity)))
>>> X = pd.DataFrame(np.c_[np.ones(300), rng.normal(size=(300, 2))])
>>> spectrum = direction_spectrum(influence(counts, X))
>>> bool(spectrum.iloc[1] < 0.2)  # rank-one counts: one direction
True

Designs

Probability designs: which entities to label, and with what probability.

Two goals, and neither is traded for the other. Unbiasedness is a hard constraint: every entity that contributes to the estimand must have a strictly positive chance of being labeled, and must be weighted by the reciprocal of that chance. Minimum standard error is the objective, optimized subject to it.

The square-root rule below satisfies the constraint automatically – it assigns pi_j > 0 to every entity with b_j > 0 without needing a floor – and minimizes the exact design risk among all Poisson designs of the same expected cost. An entity with b_j = 0 contributes nothing to the estimand, so leaving it unlabeled costs no bias at all, and it is the one case where pi_j = 0 is correct rather than merely convenient.

class fewlab.designs.LabelDesign(probabilities, scores, risk, method, diagnostics=<factory>)

Bases: object

A probability design over labelable entities.

This is a plan, not a draw. It fixes the probability with which each entity would be labeled; realizing it produces a LabelSample, which is the only thing an estimator will accept.

Parameters:

• probabilities (Series)

• scores (Series)

• risk (float)

• method (str)

• diagnostics (dict)

probabilities

First-order inclusion probabilities, indexed by entity.

Type:

pandas.Series

scores

The contribution scores b_j the design was built from.

Type:

pandas.Series

risk

Exact expected loss sum_j (1 - pi_j) / pi_j * b_j under Poisson sampling with these probabilities.

Type:

float

method

Name of the rule that produced the probabilities.

Type:

str

diagnostics

Method-specific detail, such as the certainty-core size.

Type:

dict

property budget: float

Expected number of labels purchased.

Returns:

Sum of the inclusion probabilities.

property certain: Index

Entities labeled with probability one.

Returns:

Index of the certainty core.

fewlab.designs.design_risk(scores, pi)

Compute the Poisson design risk sum_j (1 - pi_j) / pi_j * b_j.

This is exact when entities are drawn as independent Bernoulli trials and scores include true squared label norms. Fixed-size and balanced designs involve joint inclusion probabilities; this expression is a planning criterion for them, not a universal upper bound.

Parameters:

• scores (Series | ndarray) – Contribution scores b_j.

• pi (Series | ndarray) – Inclusion probabilities, aligned with scores.

Returns:

The expected loss. Infinite if any entity with b_j > 0 has pi_j = 0, which is the arithmetic signature of a biased design.

Return type:

float

fewlab.designs.min_cost_for_risk(scores, risk_limit, *, cost=None)

Minimize expected labeling cost for a Poisson risk ceiling.

Solve min sum(c_j pi_j) subject to sum(b_j (1 - pi_j) / pi_j) <= risk_limit. The solution is pi_j = min(1, lambda sqrt(b_j / c_j)). For an uncapped tail T, lambda = sum_T sqrt(b_j c_j) / (risk_limit + sum_T b_j). Sorting the ratios identifies the certainty core, including an empty core.

The risk is exact for independent Bernoulli draws when the scores include the true squared label norms. Unit norms are exact for a complete one-hot taxonomy and conservative for an unknown binary indicator. Estimated norms give a planning criterion. A scalar contrast targets its variance; joint coefficient loss targets the sum of coefficient variances. This ceiling excludes regression sampling uncertainty and does not bound realized error.

Parameters:

• scores (Series) – Finite, nonnegative contribution scores indexed by unique entity.

• risk_limit (float) – Finite, nonnegative ceiling in squared estimand units.

• cost (Series | ndarray | float | None) – Positive finite labeling costs. None means unit costs. A Series must contain exactly the score entities and is aligned by index.

Returns:

A design for Study.draw(…, fixed_size=False). Its budget is the expected label count; diagnostics also report expected cost and risk. Zero-score entities have zero probability. Guarantees apply only to the scored estimand, not to other coefficients or to balanced draws.

Raises:

• ValueError – If scores, costs, or the ceiling are invalid.

• ArithmeticError – If finite probabilities satisfying the ceiling cannot be represented numerically.

Return type:

LabelDesign

Examples

>>> import pandas as pd
>>> from fewlab import min_cost_for_risk
>>> design = min_cost_for_risk(pd.Series([1.0] * 10), risk_limit=40.0)
>>> round(design.budget, 6), round(design.risk, 6)
(2.0, 40.0)

fewlab.designs.poisson_pps(scores, budget, *, cost=None, pi_min=0.0)

Build the risk-minimizing design for a fixed labeling budget.

Minimizing sum_j (1 - pi_j) / pi_j * b_j subject to sum_j c_j pi_j = B and pi_j <= 1 gives

pi_j* = min(1, lambda * sqrt(b_j / c_j))

for the lambda that exhausts the budget. The cap at one is not a separate rule for “big” entities: it falls out of the constraint, and the set of entities it binds on is the certainty core. In a heavy-tailed population – website visits, merchant transactions, surname frequencies – that core absorbs most of the budget, which is the whole reason the design pays.

lambda is solved exactly rather than by search. Sorting the ratios descending, a proposed core leaves a tail with lambda = (B - sum_core c_j) / sum_tail sqrt(b_j c_j). Start with an empty core and promote an entity only when its proposed probability exceeds one.

Parameters:

• scores (Series) – Contribution scores b_j from contribution_scores.

• budget (float) – Total labeling budget. With uniform cost this is the expected number of entities labeled.

• cost (Series | ndarray | float | None) – Per-entity cost c_j of acquiring one label. None means uniform cost, under which the rule reduces to pi_j proportional to sqrt(b_j).

• pi_min (float) – Optional floor applied to entities with b_j > 0. The default of zero is deliberate: the square-root rule already gives every contributing entity a strictly positive probability, so no floor is needed for unbiasedness, and any floor moves the design away from minimum risk. Set it only to bound the largest weight.

Returns:

The design, carrying its probabilities and its exact risk.

Raises:

ValueError – If the budget is not positive or exceeds what the contributing entities can absorb.

Return type:

LabelDesign

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import poisson_pps
>>> b = pd.Series(np.arange(1, 21, dtype=float) ** 2)
>>> design = poisson_pps(b, budget=5)
>>> round(design.budget, 6)
5.0
>>> bool((design.probabilities > 0).all())
True

fewlab.designs.core_plus_tail(scores, budget)

Build the best certainty-core-plus-uniform-tail design.

Label the H highest-contribution entities with certainty and give every remaining entity the same probability p_H = (B - H) / (m - H). With S_H = sum_{j>H} b_j the Poisson risk collapses to a closed form,

R(H) = (m - B) / (B - H) * S_H

so the optimal core size is found by evaluating at most B cutoffs rather than guessed at. Promoting one more entity into the core pays exactly when b_{H+1} > S_H / (B - H).

R(H) is a planning criterion, and the distinction matters. It is the exact risk when the tail is drawn as independent Bernoulli trials at p_H. For a fixed-size simple random tail, the true variance instead uses contributions centered on their tail mean and a finite-population correction. Depending on the contributions it can exceed the Poisson criterion. The core size selected here is optimal for the Poisson design only.

This design is never better than poisson_pps, which is free to vary the tail probabilities. It is worth having because a uniform tail is a simple random sample: it is easy to explain to a labeling vendor, easy to audit, and its variance can be estimated by anyone with an introductory survey course. Use it when transparency is worth more than the last few percent of efficiency.

Parameters:

• scores (Series) – Contribution scores b_j from contribution_scores.

• budget (int) – Number of entities to label. Must be a positive integer no larger than the number of contributing entities.

Returns:

The design, with the chosen core size and the full R(H) curve in its diagnostics.

Raises:

ValueError – If the budget is not a positive integer, or exceeds the number of entities that contribute to the estimand.

Return type:

LabelDesign

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import core_plus_tail
>>> b = pd.Series(1.0 / np.arange(1, 101, dtype=float) ** 2)
>>> design = core_plus_tail(b, budget=20)
>>> design.diagnostics["core_size"]
10
>>> round(design.budget, 6)
20.0

A certainty core is not automatic. It appears only when the head is concentrated enough that b_{H+1} exceeds the average remaining contribution; on a population that decays more gently the whole budget goes to the tail, and the design says so rather than forcing a core.

>>> gentle = pd.Series(np.arange(1, 101, dtype=float)[::-1] ** 3)
>>> core_plus_tail(gentle, budget=20).diagnostics["core_size"]
0

fewlab.designs.simple_random(scores, budget)

Build the uniform-probability design, as a benchmark.

Every contributing entity gets the same probability B / m. This is what you do when you ignore the fact that entities differ in how much they matter to the estimand, and it is the baseline any efficiency claim should be measured against.

Parameters:

• scores (Series) – Contribution scores b_j, used for the risk calculation and to identify which entities contribute at all.

• budget (int) – Number of entities to label.

Returns:

The uniform design and its exact risk.

Raises:

ValueError – If the budget is not a positive integer, or exceeds the number of contributing entities.

Return type:

LabelDesign

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import poisson_pps, simple_random
>>> b = pd.Series(np.arange(1, 201, dtype=float) ** 2)
>>> srs = simple_random(b, budget=20)
>>> pps = poisson_pps(b, budget=20)
>>> pps.risk < srs.risk
True

fewlab.designs.robust_pps(scores, budget, *, prevalence, cost=None, pi_min=0.0001)

Build the design that minimizes worst-case risk under a prevalence cap.

poisson_pps needs the squared label norms, which for a binary trait means knowing which entities carry it. When that is unknown but the trait is known to be rare, the worst case is tractable. Suppose at most a fraction alpha of entities are positive, so at most M = ceil(alpha m) of them, and write t_j = (1 - pi_j)/pi_j * w_j for the label-free contribution. Then

max over label configurations of the risk = sum of the M largest t_j,

since an adversary spends its M positives on the largest terms. That maximum is convex in pi, despite being a maximum over a combinatorial set, by the epigraph identity

sum of the M largest t = min over tau of [ M tau + sum_j (t_j - tau)_+ ],

a minimum of jointly convex functions. Minimizing it under the budget is therefore a convex program, solved here by bisecting the budget multiplier inside a golden-section search on tau.

Two consequences follow in a line each and both are checked by the test suite rather than assumed. At alpha = 1 the sum of the m largest terms is the whole sum, so the program becomes poisson_pps with unit label norms and returns the same design. And on the worst case this design is never beaten by poisson_pps, since it minimizes exactly that quantity; it is beaten on the average case, which is what the insurance costs.

Parameters:

• scores (Series) – Label-free contributions w_j = h_j’ W h_j, from contribution_scores with the default unit label norms.

• budget (float) – Total labeling budget.

• prevalence (float) – Upper bound alpha on the fraction of entities carrying the trait, in (0, 1].

• cost (Series | ndarray | float | None) – Per-entity labeling cost, or None for uniform.

• pi_min (float) – Lower clip. Unlike poisson_pps, a small positive floor is the default here, because the worst-case objective is flat in entities the adversary will not choose and so nothing else keeps their probabilities away from zero.

Returns:

The design, with the fitted tau and the worst-case risk among its diagnostics.

Raises:

ValueError – If prevalence is outside (0, 1] or the budget is not positive.

Return type:

LabelDesign

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import robust_pps
>>> w = pd.Series(np.exp(np.linspace(0, 4, 200)))
>>> design = robust_pps(w, budget=40, prevalence=0.05)
>>> round(design.budget, 4)
40.0
>>> design.diagnostics["n_positive"]
10

Estimation

Design-based estimation from a realized label sample.

Nothing here accepts a set of entities on its own. An estimator needs the probability with which each entity could have been labeled, not merely which ones were, so every function takes a LabelSample carrying the probabilities the sampler actually delivered. That is the whole reason the type exists: a ranking of entities by importance is not a probability sample, and no probability assigned to it after the fact makes it one.

Two variances are reported, and they answer different questions:

Var(beta_hat) = sigma^2 (X’X)^-1 + sum_j (1-pi_j)/pi_j b_j

fixed by your n units fixed by your labeling budget

Only the second moves when you spend more on labels. Their ratio is a stopping rule: once the design term is small next to the regression term, further labels buy almost nothing.

class fewlab.estimator.LabelSample(items, probabilities, method, diagnostics=<factory>)

Bases: object

A realized draw from a probability design.

Parameters:

• items (Index)

• probabilities (Series)

• method (str)

• diagnostics (dict)

items

The entities that were labeled.

Type:

pandas.Index

probabilities

The probabilities the sampler delivered, over every entity in the frame. These may differ from the ones a design requested – a fixed-size sampler rescales them – and it is the delivered ones that make the estimator unbiased.

Type:

pandas.Series

method

How the draw was made.

Type:

str

diagnostics

Sampler detail.

Type:

dict

property weights: Series

Horvitz-Thompson weights 1 / pi_j for the labeled entities.

Returns:

Weights indexed by labeled entity.

Raises:

ValueError – If any labeled entity has zero delivered probability, which would make the weight infinite and the design invalid.

class fewlab.estimator.Estimate(coefficients, design_se, regression_se, total_se, fitted, design_cov, bias_bound, bias_estimate, diagnostics=<factory>)

Bases: object

Coefficients, both variance components, and any bias that was bought.

Parameters:

• coefficients (DataFrame)

• design_se (DataFrame)

• regression_se (DataFrame | None)

• total_se (DataFrame | None)

• fitted (DataFrame)

• design_cov (ndarray)

• bias_bound (float)

• bias_estimate (DataFrame | None)

• diagnostics (dict)

coefficients

(p, q) estimate of B, features by label class.

Type:

pandas.DataFrame

design_se

Label-sampling standard errors from the chosen covariance estimator. The Poisson covariance estimator is design-unbiased for independent draws; other designs require different estimators.

Type:

pandas.DataFrame

regression_se

Standard errors from having finitely many units, or None when the inputs needed to compute them were not supplied.

Type:

pandas.DataFrame | None

total_se

The two components combined in quadrature, or None.

Type:

pandas.DataFrame | None

fitted

(n, q) estimated unit-level shares Yhat.

Type:

pandas.DataFrame

design_cov

Full (pq, pq) design covariance of vec(B).

Type:

numpy.ndarray

bias_bound

Worst-case Frobenius norm of the bias, zero when unshrunk.

Type:

float

bias_estimate

Proxy-based point estimate of the bias, when a proxy was given and shrinkage applied.

Type:

pandas.DataFrame | None

diagnostics

Includes variance_ratio, the share of total variance the labeling budget is responsible for.

Type:

dict

summary()

Tabulate the estimate against its uncertainty.

Returns:

One row per feature and label class, with the coefficient, each standard error component, and the t-statistic against the total.

Return type:

DataFrame

fewlab.estimator.estimate(infl, sample, labels, *, counts=None, X=None, proxy=None, alpha=1.0, variance='poisson')

Estimate B from a realized label sample.

The plain Horvitz-Thompson estimator is

Bhat = sum_{j in S} (1/pi_j) h_j a_j’

which is design-unbiased for B = sum_j h_j a_j’ because each entity enters with probability pi_j and is weighted by its reciprocal.

Given a cheap proxy mu_j available for every entity, the difference estimator

Bhat = sum_j h_j mu_j’ + sum_{j in S} (1/pi_j) h_j (a_j - mu_j)’

is used instead. It is exactly unbiased whatever the proxy says, because the proxy term is a known constant and the sampled term corrects it without bias. A poor proxy costs variance and never bias, so there is no risk in supplying one.

Parameters:

• infl (Influence) – Influence decomposition for the frame.

• sample (LabelSample) – The realized draw, carrying delivered probabilities.

• labels (Series | DataFrame) – Observed labels for the sampled entities. A Series for a single label, a frame with one column per class for a taxonomy.

• counts (DataFrame | CountMatrix | None) – Count matrix, needed for fitted values and the regression variance component. Optional.

• X (DataFrame | None) – Feature matrix, needed for the regression variance component. Optional.

• proxy (DataFrame | None) – (m, q) predictions for every entity in the frame, with the same columns as labels. Optional.

• alpha (Series | float) – Shrinkage applied to the sampled correction. The default of 1.0 is the unbiased estimator. See trim_alpha before departing from it.

• variance (str) – Which variance estimator to use. “poisson” is exact for an independent Bernoulli draw; it has no universal conservatism guarantee for other designs. “hajek” accounts for a fixed sample size, and “balanced” also accounts for the auxiliaries a cube draw held fixed; use them only when the sample was drawn that way, since on a Poisson draw they understate.

Returns:

The estimate, its two variance components, and any bias incurred.

Raises:

ValueError – If shapes or indices are inconsistent, or a sampled entity has zero delivered probability.

Return type:

Estimate

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import Study
>>> rng = np.random.default_rng(0)
>>> counts = pd.DataFrame(rng.poisson(4, (150, 30)))
>>> X = pd.DataFrame(np.c_[np.ones(150), rng.normal(size=(150, 2))])
>>> study = Study(counts, X)
>>> sample = study.draw(study.design(budget=12), random_state=1)
>>> labels = pd.Series(rng.integers(0, 2, len(sample)), index=sample.items)
>>> est = study.estimate(sample, labels)
>>> est.coefficients.shape
(3, 1)

fewlab.estimator.trim_alpha(sample, weight_cap)

Build shrinkage factors that cap Horvitz-Thompson weights.

Weight trimming is the classical way to trade a little bias for a large variance reduction when a few units carry enormous weights. Capping the weight at w is exactly alpha_j = min(1, w * pi_j).

Parameters:

• sample (LabelSample) – The realized draw whose weights are to be capped.

• weight_cap (float) – The largest weight to allow.

Returns:

Shrinkage factors indexed by entity, ready to pass to estimate.

Raises:

ValueError – If the cap is not positive.

Return type:

Series

References

Potter, F. (1990). A study of procedures to identify and trim extreme sampling weights. Proceedings of the Survey Research Methods Section, American Statistical Association, 225-230.

Beaumont, J.-F., Haziza, D., & Ruiz-Gazen, A. (2013). A unified approach to robust estimation in finite population sampling. Biometrika, 100(3), 555-569.

Labels as covariates

When the purchased label enters the regression as a covariate.

Everything else in this package assumes the label builds the outcome: the unit level share y_i = sum_j V_ij a_j is regressed on features X, and the Horvitz-Thompson estimator of y carries its unbiasedness straight through to the coefficients.

That guarantee does not survive putting the constructed share on the other side of the equation. Suppose the target is gamma in

z_i = gamma * y_i + x_i’ delta + eps_i

and y is replaced by an unbiased estimate yhat = y + u. Ordinary least squares now regresses on a covariate measured with error, so gamma is attenuated rather than merely noisy, and the attenuation does not shrink as the panel grows: it shrinks only as the labeling budget grows. Unbiasedness of yhat is no protection at all, which is the trap this module exists to avoid.

What makes the problem tractable here is that the error is not a nuisance of unknown size. The design chose pi, so Var(u) is known up to the labels themselves, and the sample estimates it without bias. Subtracting it from the cross-product removes the first-order attenuation.

class fewlab.covariate.ShareEffect(coefficient, naive, attenuation, design_variance, diagnostics=<factory>)

Bases: object

The coefficient on a label-derived share, and how far it was attenuated.

Parameters:

• coefficient (float)

• naive (float)

• attenuation (float)

• design_variance (float)

• diagnostics (dict)

coefficient

The attenuation-corrected estimate of gamma.

Type:

float

naive

The uncorrected estimate, which regresses z on yhat directly.

Type:

float

attenuation

The factor by which the naive estimate is shrunk toward zero, in (0, 1]. The corrected estimate is the naive one divided by it.

Type:

float

design_variance

The estimated E[u’ M_X u] that was subtracted.

Type:

float

diagnostics

Sample size and the residual sum of squares involved.

Type:

dict

fewlab.covariate.estimate_share_effect(infl, counts, X, sample, labels, outcome, *, correct=True)

Estimate the coefficient on a label-derived share, correcting attenuation.

Residualizing on the controls, the naive slope is ytilde’ztilde / ytilde’ytilde. Its denominator is inflated by the design noise, since E[yhat’ M_X yhat] = y’ M_X y + E[u’ M_X u], and that inflation is what shrinks the coefficient. The correction subtracts an unbiased estimate of the second term,

sum_{j in S} (1 - pi_j) / pi_j^2 * a_j^2 * (v_j’ M_X v_j),

which is computable because the design fixed pi and the sample reveals a_j.

Parameters:

• infl (Influence) – Influence decomposition for the frame.

• counts (DataFrame | CountMatrix) – The count matrix.

• X (DataFrame) – Control matrix, including an intercept.

• sample (LabelSample) – The realized label draw, carrying delivered probabilities.

• labels (Series) – Observed labels for the sampled entities. Scalar labels only; a covariate is one number per unit.

• outcome (Series) – The dependent variable z, indexed by unit.

• correct (bool) – Whether to apply the attenuation correction. False returns the naive estimate, for comparison.

Returns:

The estimate, the naive value, and the attenuation between them.

Raises:

ValueError – If the corrected denominator is non-positive, which means the design noise accounts for the whole of the observed variation and no coefficient is identified at this budget.

Return type:

ShareEffect

Note

The correction removes the first-order attenuation but is a ratio with a reduced denominator, so it overshoots slightly: on a synthetic panel where the naive estimate was attenuated by 32 percent, the corrected one ran about 11 percent high. Treat it as a large improvement rather than as exact, and obtain intervals by redrawing the design rather than from a closed form, which this function deliberately does not offer.

Counts

A count matrix that stays sparse, because real usage data is.

A browsing panel is the motivating case and it sets the shape: 1,134 users by 34,078 domains, 0.27% of cells non-zero. Dense, that is 309 MB of mostly zeros; as triples it is a few megabytes. Panels arrive in long form anyway – one row per user-entity pair – so this holds them that way and implements only the four operations the rest of the package needs.

class fewlab.counts.CountMatrix(rows, cols, values, index, columns)

Bases: object

Counts in coordinate form, with dense operations computed on demand.

Parameters:

• rows (ndarray)

• cols (ndarray)

• values (ndarray)

• index (Index)

• columns (Index)

rows

Row position of each stored count.

Type:

numpy.ndarray

cols

Column position of each stored count.

Type:

numpy.ndarray

values

The counts themselves.

Type:

numpy.ndarray

index

Unit identifiers, indexed by row position.

Type:

pandas.Index

columns

Entity identifiers, indexed by column position.

Type:

pandas.Index

property shape: tuple[int, int]

Matrix dimensions.

Returns:

Units by entities.

property nnz: int

Number of stored counts.

Returns:

Count of non-zero entries.

property density: float

Share of cells that are non-zero.

Returns:

Density in [0, 1].

classmethod from_long(frame, *, unit, entity, count, entities=None)

Build from one row per unit-entity pair.

Parameters:

• frame (DataFrame) – Long-format counts.

• unit (str) – Column naming the unit.

• entity (str) – Column naming the labelable entity.

• count (str) – Column holding the count.

• entities (Index | None) – Optional full entity frame. Supply it when some entities have zero recorded activity but should still belong to the population; they contribute nothing to the estimand and will be assigned zero probability, which costs no bias.

Returns:

The sparse count matrix.

Raises:

ValueError – If counts are negative.

Return type:

CountMatrix

classmethod from_dense(frame)

Build from a dense DataFrame.

Parameters:

frame (DataFrame) – Dense count matrix, units by entities.

Returns:

The sparse count matrix, keeping only non-zero cells.

Return type:

CountMatrix

row_totals()

Sum each unit’s counts across entities.

Returns:

Row totals, one per unit.

Return type:

ndarray

project(weights)

Compute weights’ C, the entity-wise projection of unit weights.

With weights = X / T this is G = X’ V, the only contraction of the count matrix the design needs, and it costs one pass over the stored entries rather than a dense multiply.

Parameters:

weights (ndarray) – (n, p) array of per-unit weights.

Returns:

The (p, m) projection.

Return type:

ndarray

column_sq_norms(scale)

Compute sum_i (scale_i * C_ij)^2 for every entity.

With scale = 1/T this is ||v_j||^2, the total squared share entity j accounts for across units. It is the piece the design needs when the label enters a regression as a covariate rather than as the outcome, because there the identifying variation is the part of v_j orthogonal to the controls rather than the part aligned with them.

Parameters:

scale (ndarray) – Per-unit multiplier, one entry per unit.

Returns:

One squared norm per entity.

Return type:

ndarray

dense_columns(positions)

Materialize a subset of columns.

Parameters:

positions (ndarray) – Column positions to extract.

Returns:

An (n, len(positions)) dense array.

Return type:

ndarray

select_rows(keep)

Drop units.

Parameters:

keep (ndarray) – Boolean mask over units.

Returns:

A matrix restricted to the kept units.

Return type:

CountMatrix

fewlab.counts.as_counts(counts)

Accept either representation.

Parameters:

counts (DataFrame | CountMatrix) – A dense DataFrame or an already-sparse matrix.

Returns:

The sparse representation.

Return type:

CountMatrix

Balanced sampling

Balanced sampling that delivers the inclusion probabilities it promises.

Everything fewlab estimates is weighted by 1 / pi, which is unbiased only if item j really is included with probability pi[j]. Two things make that hard to arrange at once:

• Fixed size. The sample must contain exactly budget items, so the inclusion indicators cannot be independent.

• Balance. The point of the exercise is to make the sample representative – sum_j (I_j / pi_j - 1) g_j near zero – which correlates the indicators further.

The obvious approaches sacrifice the first requirement to get the other two. A draw of budget items without replacement with probabilities proportional to pi (rng.choice(..., replace=False, p=pi)) does not have inclusion probabilities pi: successive draws condition on what came before, and the gap is large, not a rounding artefact. On the package’s own test fixture that draw gave items with pi == 1 an inclusion rate of 0.75. Repairing a sample by greedily swapping items to improve balance is worse still, since the swap decisions depend on the data and distort pi further – the same fixture fell to 0.45.

The cube method solves all three at once. Write the problem geometrically: the vector of inclusion probabilities is a point in the cube [0, 1]^m, and a sample is a vertex. Balance and fixed size are linear constraints, so they cut out an affine subspace through that point.

Flight phase. Repeatedly take a random step along a direction lying in the null space of the constraint matrix, as far as the cube permits, choosing between the two directions with probabilities that make the walk a martingale. Because the step stays in the null space the constraints hold exactly throughout; because the walk is a martingale the expectation is exactly the pi we started with; and because each step runs until it hits a face, at least one coordinate is settled to 0 or 1 every iteration. The whole flight is therefore O(m) steps.

Landing phase. The flight stalls once fewer coordinates remain unsettled than there are constraints. The remaining constraints are then dropped one at a time, least important first, and the flight resumes. Balance degrades gracefully; pi and the fixed size do not degrade at all, because the sample-size constraint is the last one standing and is never dropped.

Reference:

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

fewlab.cube.cube_sample(pi, constraints, rng)

Draw a balanced sample whose inclusion probabilities are exactly pi.

Parameters:

• pi (ndarray) – Target inclusion probabilities. Must sum to an integer for the sample size to be fixed; see scale_pi_to_budget().

• constraints (ndarray) – Balancing rows, most dispensable first. The walk preserves constraints @ state exactly for as long as it can, and drops rows from the top when it stalls, so the last row is the one that survives to the end and should be the sample-size constraint.

• rng (Generator) – Source of randomness.

Returns:

Boolean mask of the selected items.

Return type:

np.ndarray

Raises:

ValueError – If constraints has a column for every item is violated, that is, if its width does not match pi.

fewlab.cube.scale_pi_to_budget(pi, budget, *, max_iter=1000)

Rescale inclusion probabilities to sum to budget, respecting the cap.

sum(pi) is the expected sample size, so a fixed-size design of budget items is only coherent when the probabilities sum to budget. Scaling by a constant would push large entries above one, which is not a probability, so entries that would exceed one are pinned there and the remaining budget is redistributed over the rest. Repeating that to a fixed point is the standard construction, and it terminates because each pass pins at least one entry.

Parameters:

• pi (ndarray) – Non-negative target probabilities.

• budget (int) – Desired sample size, at most len(pi).

• max_iter (int) – Guard on the redistribution loop.

Returns:

Probabilities in [0, 1] summing to budget.

Return type:

np.ndarray

Raises:

ValueError – If budget is negative, exceeds the number of items, or exceeds the number of items with positive probability.

Examples: >>> import numpy as np >>> scaled = scale_pi_to_budget(np.array([0.9, 0.2, 0.2, 0.2]), 2) >>> float(scaled.sum()) 2.0 >>> bool((scaled <= 1.0).all()) True

Calibration

Weight calibration methods for optimal survey sampling.

This module implements GREG (Generalized Regression) calibration and related techniques for adjusting sampling weights to match known population totals.

fewlab.calibration.calibrate_weights(pi, g, selected, pop_totals=None, *, distance='chi2', ridge=1e-08, nonneg=True)

Compute calibrated weights via GREG/Deville-Särndal calibration.

Parameters:

• pi (Series) – Inclusion probabilities for all items (index = item names).

• g (ndarray) – Regression projections g_j = X^T v_j for all items (shape (p, m)).

• selected (Sequence[str] | Index) – Item identifiers drawn in the sample.

• pop_totals (ndarray | None) – Known population totals (shape (p,)); defaults to g.sum(axis=1).

• distance (str) – Calibration distance measure; currently only “chi2” is supported.

• ridge (float) – Ridge regularization parameter for numerical stability.

• nonneg (bool) – Whether to enforce non-negative calibrated weights.

Returns:

Calibrated weights indexed by the selected items.

Raises:

• NotImplementedError – If distance is not “chi2”.

• ValueError – If pop_totals has the wrong shape.

Return type:

Series

Notes

The closed-form solution for chi-square distance is w* = d_S + G_S^T (G_S G_S^T + ridge I)^{-1} (t - G_S d_S), where d_S are the base weights.

References

Deville, J.-C., & Särndal, C.-E. (1992). Calibration estimators in survey sampling. Journal of the American Statistical Association, 87(418), 376-382.

Fewest labels for a fixed standard error

Inclusion probabilities under row-wise standard-error constraints.

Instead of fixing a budget and minimising variance, this module fixes a per-row variance ceiling and minimises the expected number of labels.

class fewlab.rowse.RowSEResult(probabilities, expected_labels, max_violation, tolerance, iterations, best_iteration, feasible)

Bases: object

Probabilities and convergence detail from the row-SE optimizer.

Parameters:

• probabilities (Series)

• expected_labels (float)

• max_violation (float)

• tolerance (float)

• iterations (int)

• best_iteration (int)

• feasible (bool)

probabilities

Inclusion probabilities, indexed by entity.

Type:

pandas.Series

expected_labels

The objective, sum_j pi_j.

Type:

float

max_violation

Largest remaining constraint violation.

Type:

float

tolerance

Feasibility tolerance used.

Type:

float

iterations

Dual steps taken.

Type:

int

best_iteration

Step at which the returned design was found.

Type:

int

feasible

Whether every row-SE ceiling is met.

Type:

bool

to_series()

Return just the probabilities.

Returns:

The inclusion probabilities.

Return type:

Series

fewlab.rowse.row_se_min_labels(counts, eps2, *, pi_min=0.0001, max_iter=8000, tol=1e-06, return_result=False, raise_on_failure=False)

Minimize expected labels subject to row-wise SE limits.

The routine solves:

` minimize   sum_j pi_j subject to sum_j q_ij / pi_j <= eps2_i + sum_j q_ij,  q_ij = (c_ij / T_i)^2 `

Parameters:

• counts (DataFrame) – Non-negative count matrix with units as rows and items as columns.

• eps2 (ndarray | Series) – Row-wise squared standard-error tolerance; scalar applies to every row.

• pi_min (float) – Minimum allowable inclusion probability.

• max_iter (int) – Maximum optimization iterations.

• tol (float) – Convergence tolerance for constraint violations.

• return_result (bool) – If True, return a RowSEResult with diagnostics.

• raise_on_failure (bool) – If True, raise a ValidationError when constraints remain violated.

Returns:

Probability series if return_result is False (default) or a RowSEResult with diagnostics when return_result is True.

Raises:

ValidationError – If inputs are invalid or raise_on_failure is True and the constraints remain violated after optimization.

Return type:

RowSEResult | Series

See also

pi_aopt_for_budget: A-optimal probabilities for a fixed budget. items_to_label: Deterministic selection without SE constraints.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from fewlab import row_se_min_labels
>>> rng = np.random.default_rng(0)
>>> counts = pd.DataFrame(rng.poisson(10, (200, 20)))
>>> pi = row_se_min_labels(counts, eps2=0.3**2)
>>> bool(0 < pi.sum() <= counts.shape[1])
True