andrey.bootstrap

andrey.bootstrap(method, data, *, n_resamples=100, seed=None, threshold=0.5, strict=False, resampling='subsample', sample_fraction=0.5, block_length=None, keep_fits=False, **method_kwargs)

Refit a method on resamples of the rows; return each edge’s frequency and an averaged graph.

Each resample is a subsample, sample_fraction of the rows of data drawn without replacement, or a bootstrap resample, as many rows drawn with replacement (resampling), and method is fitted on it. The result counts, for every pair of variables, the fits that give each endpoint-mark pair (->, <-, --, <->, o->, …) and builds the averaged graph: the vote of the fits, made a valid graph of the method’s kind [1] [2] [3]. A method with a row-aligned input (cdnod’s domain_index, gies’s interventions) resamples that input with the rows. A time-series method on one series subsamples one stretch of consecutive rows (and bootstraps blocks of consecutive rows), longitudinal_lingam resamples whole units, and multi_group_direct_lingam resamples each dataset on its own.

Parameters:

• method (str | Callable[..., Any]) – A method’s name, such as "pc", or the callable itself, such as andrey.pc.

• data (ArrayLike | Sequence[ArrayLike]) – The data, as the method takes it. multi_group_direct_lingam and longitudinal_lingam take a sequence of arrays.

• n_resamples (int) – The number of resamples and fits. n_boot, the former name, works with a DeprecatedArgumentWarning for one more release.

• seed (int | None) – The seed of the resamples, a non-negative integer. None uses ANDREY_SEED, or 0 when it is unset. Resample b depends on seed and b alone. A method that takes a seed fits resample b with a seed derived from this one and b, so the method’s own seed cannot be set here.

• threshold (float) – A pair of the averaged graph must be adjacent in more than this share of the fits, or in all of them, in [0.5, 1]. Its most frequent edge must also be in more fits than no edge is, which needs more than half of them, so a value below 0.5 would change nothing.

• strict (bool) – What happens when the vote of the fits has no valid graph of the method’s kind: False returns the valid resample graph nearest the others with a warning, and True raises ValueError.

• resampling (str) – How a resample is drawn. "subsample" draws round(sample_fraction * n_samples) distinct rows, kept in their original order; for a time-series method on one series it is one stretch of consecutive rows at a random start [4] [5]. "bootstrap" draws n_samples rows with replacement; for a time-series method on one series it draws blocks of block_length consecutive rows (a moving block bootstrap).

• sample_fraction (float) – The share of the rows in a subsample, in (0, 1], so that round(sample_fraction * n_samples) is at least 2; for longitudinal_lingam, the share of the units, and for multi_group_direct_lingam, the share of each dataset’s rows. Not used by "bootstrap".

• block_length (int | None) – For resampling="bootstrap" of a time-series method on one series, the number of consecutive rows in a block, from 1 to n_samples. None uses max(ceil(n_samples ** (1/3)), 10 * tau_max) for a method with a tau_max, and ceil(n_samples ** (1/3)) for the others, at most n_samples. A block break puts rows that were not neighbors side by side, which weakens a link that many steps long or longer, so a block should be at least ten times the longest lag. Not for any other resampling.

• keep_fits (bool) – Keep every fit in the fits property.

• **method_kwargs (Any) – The method’s other arguments: its keyword parameters, and a row-aligned input such as domain_index=.

Returns:

The frequencies and the averaged graph; for multi_group_direct_lingam, one per group. metadata records the resampling, the sample_fraction of a subsample, the block_length of a block bootstrap, the tie_break rule in effect (andrey.config), and averaged_graph: "vote", "nearest_resample" when the averaged graph is the valid resample graph nearest the others, or "nearest_resample_invalid" when no resample is a valid graph of the method’s kind.

Return type:

BootstrapOutput | list[BootstrapOutput]

Raises:

• ValueError – If method is not a method of Andrey, an argument is out of range, a row-aligned input has a different number of rows from data, block_length is given for a resampling that does not use blocks, strict=True and the vote has no valid graph of the method’s kind, or a fit raises it: the error then names the resample in a note.

• NotImplementedError – For gin, whose number of latent nodes varies from fit to fit.

• RuntimeError – If a worker process dies.

Warns:

AndreyWarning – Each distinct warning the fits raised, once; and, with strict=False, when the vote has no valid graph of the method’s kind, so the averaged graph is the valid resample graph nearest the others, or, when no resample is valid, the resample graph nearest the others.

Notes

Subsampling and the bootstrap. A bootstrap resample repeats some rows and leaves out others, so a conditional-independence test or a score sees more evidence than the resample holds. A method built on such tests, such as pc, keeps weak and false edges in more resamples than it does on data, and an orientation that needs two variables to stay non-adjacent, a collider, is found in fewer. A subsample holds distinct rows only, so it has no such excess. Banerjee, Andrews, and Kummerfeld [3] compared both schemes with two score searches on linear Gaussian data: no scheme was best everywhere, the bootstrap at the full size performed poorly, and subsampling half of the rows performed about as well as a bootstrap scored at the effective sample size, about half the size. The default draws half of the rows, as stability selection does [4], and it changes only the draw, so it means the same for every method.

On linear data with 20 and 50 variables and 1,000 and 5,000 rows (pc, ges, and direct_lingam; 10 datasets per cell), the averaged graph of half-size subsamples had no more errors (SHD) than the bootstrap’s in any of the 10 cells, and fewer in 5; for ges with 50 variables and 1,000 rows, 5.0 fewer on average (95% interval 2.7 to 7.9). For ges with 20 variables and 1,000 rows, the averaged graph had 18 wrong edge types instead of 32, and 9 of the 10 averaged graphs were valid CPDAGs instead of 3. The adjacency frequencies had a lower calibration error in 9 of the 10 cells, and a fit took 0.56 to 0.95 of the time, since it has half the rows. The one cell where the point estimate was worse, pc with 20 variables and 1,000 rows, was 0.43 errors above the bootstrap’s over 30 datasets (interval -0.23 to 1.13). A larger sample_fraction gives fits closer to the fit on data, so the resamples agree more and the frequencies say less about how much the fit varies.

For a time-series method on one series, a subsample is one stretch of consecutive rows, which keeps the lags inside it [5]. With pcmci, pcmciplus, lpcmci, and varma_lingam on simulated series with 5 and 10 variables and 500 and 1,000 rows (10 series per cell, 14 cells; 30 series per cell, 12 cells, for the first three), a stretch of half of the series had no more errors (SHD) than the block bootstrap in any cell, and fewer, by an interval below 0, in 9 of the 14 cells of 10 series. Those blocks were ceil(n_samples ** (1/3)) rows, shorter than ten times the longest lag.

The averaged graph. Each pair of variables takes its most frequent state in the fits, no edge counted as one of the states: the pair gets its most frequent edge when more fits have that edge than have no edge, which needs the pair adjacent in more than half of them, and threshold can ask for more. When two edges tie, the pair takes the least committal one, --, or o-o for a pag. This vote is the graph with the least mean structural Hamming distance to the fits. It is the Highest ensemble of Tetrad [6], with a tie between an edge and no edge going to no edge, as in tigramite’s majority vote [7].

A vote of CPDAGs or PAGs is often not a valid CPDAG or PAG: it can direct an edge that no v-structure compels. So the vote is then made a valid graph of the method’s kind by the published steps. A dag becomes a consistent extension, a DAG with the vote’s skeleton, v-structures, and directed edges [8]: the vote itself when it has no directed cycle, with a tie edge pointing as the columns’ tie_ranks() on data decide (tie_break, andrey.config: ranks computed from the data, or the column order under "position"), so reordering the columns relabels the graph. A cpdag becomes the CPDAG [9] of a consistent extension, so an edge that no v-structure compels is undirected. A pag becomes the PAG of the MAG that Zhang’s construction builds from it [10]. With knowledge, the edges it directs stay directed: a cpdag becomes its maximally oriented PDAG, and a pag is oriented as fci orients one under knowledge: the endpoints the knowledge sets come first, and Zhang’s orientation rules then run from them until no endpoint changes, as they do in a fit [10], [11]. A time-series result’s lag graphs keep the vote, as tigramite’s does [7]: a lag-0 graph is oriented with the lagged edges, which the steps for one graph do not see.

The vote has no consistent extension or MAG when, for example, its directed edges form a cycle or its undirected edges a chordless cycle of four or more nodes. With strict=False, the averaged graph is then the resample graph nearest the others, a graph the method returned: the fit whose state of each pair the most fits share, summed over the pairs, which is the fit with the least mean structural Hamming distance to the others (Tetrad’s Median [6]). Only a resample that is a valid graph of the method’s kind can be that fit: a cpdag with a consistent DAG extension, or a PAG that passes the steps above, because a method can return a graph labeled cpdag that no DAG extends (pc leaves conflicting orientations undirected). A warning says so, and metadata["averaged_graph"] is "nearest_resample"; it is "vote" otherwise. When no resample is valid, a warning says so, the averaged graph is the resample nearest the others, and metadata["averaged_graph"] is "nearest_resample_invalid". strict=True raises in both cases instead, as bnlearn’s cextend does by default [12]. Either way, frequency and table keep the whole vote.

On linear data with 20 and 50 variables and 500, 1,000, and 5,000 rows (pc, ges, fci, and direct_lingam; 20 datasets per cell; 100 resamples with resampling="bootstrap"), the averaged graph had fewer errors (SHD) than the rule of 0.11.0 (a pair adjacent in at least half the fits takes its most frequent edge) in 11 of the 24 cells, by an interval below 0, and more in none. It had fewer errors than the fit on all the rows in 14 cells and more in one, ges with 20 variables and 1,000 rows. The vote had a valid graph of the method’s kind in 111 of the 120 datasets of pc, 114 of fci, and all of ges and direct_lingam; the rule of 0.11.0 gave a valid CPDAG or PAG in 0 to 13 of the 20 datasets of a cell.

A frequency is a ranking of the edges, not a probability, and no error bound comes with the averaged graph. Resampling can also make a column constant; a fit that raises stops the run, and its row numbers count the rows of the resample.

Workers. The fits run on the worker processes andrey.config.n_jobs and ANDREY_N_JOBS set, as for every method: -1 uses every usable core, and 0 or 1 fits in this process. Each fit runs on one worker and one BLAS thread, and the result is the same at any worker count.

Limitations. Not benchmarked yet. Refitting is the cost: n_resamples times the method’s own, shared across the workers. A moving block bootstrap weakens a link that spans a block break, and a stretch holds only sample_fraction of the series, so each understates what a very short series can show.

References

[1]

Efron, B. (1979). Bootstrap methods: another look at the jackknife. Annals of Statistics, 7(1), 1-26.

[2]

Friedman, N., Goldszmidt, M., and Wyner, A. (1999). Data analysis with Bayesian networks: a bootstrap approach. Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, 196-205.

[3] (1,2)

Banerjee, R., Andrews, B., and Kummerfeld, E. (2025). An extensive simulation study evaluating the interaction of resampling techniques across multiple causal discovery contexts. arXiv:2503.15436.

[4] (1,2)

Meinshausen, N., and Buhlmann, P. (2010). Stability selection. Journal of the Royal Statistical Society: Series B, 72(4), 417-473.

[5] (1,2)

Politis, D. N., Romano, J. P., and Wolf, M. (1999). Subsampling. Springer.

[6] (1,2)

Tetrad, edu.pitt.dbmi.algo.resampling.ResamplingEdgeEnsemble and edu.cmu.tetrad.util.GraphSampling, github.com/cmu-phil/tetrad, development branch, 2026-08-13. Median is not in release 7.6.10.

[7] (1,2)

Debeire, K., Gerhardus, A., Runge, J., and Eyring, V. (2024). Bootstrap aggregation and confidence measures to improve time series causal discovery. Proceedings of the Third Conference on Causal Learning and Reasoning, PMLR 236.

[8]

Dor, D., and Tarsi, M. (1992). A simple algorithm to construct a consistent extension of a partially oriented graph. Technical Report R-185, Cognitive Systems Laboratory, UCLA.

[9]

Chickering, D. M. (2002). Learning equivalence classes of Bayesian-network structures. Journal of Machine Learning Research, 2, 445-498.

[10] (1,2)

Zhang, J. (2008). On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence, 172(16-17), 1873-1896.

[11]

Andrews, B., Spirtes, P., and Cooper, G. F. (2020). On the completeness of causal discovery in the presence of latent confounding with tiered background knowledge. Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics, PMLR 108, 4002-4011.

[12]

Scutari, M. (2010). Learning Bayesian networks with the bnlearn R package. Journal of Statistical Software, 35(3), 1-22.

Examples

>>> import numpy as np
>>> from andrey import bootstrap
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(300)
>>> data = np.column_stack([x, x + rng.standard_normal(300), rng.standard_normal(300)])
>>> out = bootstrap("pc", data, n_resamples=10, seed=0)
>>> out.n_resamples
10
>>> bool(out.frequency("absent")[0, 1] < 0.5)
True
>>> out.metadata["resampling"], out.metadata["sample_fraction"]
('subsample', 0.5)