andrey.boss

andrey.boss(data, *, score_func=None, lambda_value=1.0, seed=None, categorical=None, max_categories=10, knowledge=None)

Learn a CPDAG with BOSS (Best Order Score Search).

BOSS searches variable orders (permutations) using linear-Gaussian BIC (Bayesian information criterion). It converts the best-scoring order to a CPDAG, which represents an equivalence class of graphs. It is fast and accurate on linear-Gaussian data.

Parameters:

• data (ArrayLike) – Observed data, one row per sample. Coerced to float64. A DataFrame’s column names become the structure’s labels.

• score_func (str | None) – Local score objective: "local_score_BIC_from_cov", the covariance-based linear-Gaussian BIC; "local_score_BDeu" and "local_score_BIC_categorical" are the BDeu score and the multinomial BIC for categorical data coded by whole numbers, and "local_score_degen_gauss" and "local_score_cond_gauss" are the degenerate and conditional Gaussian BICs for mixed data, as for ges(). None reads the column types from categorical, or infers them, and picks "local_score_BDeu" when every column is categorical or "local_score_cond_gauss" when some are, each with a DataTypeWarning; otherwise it picks "local_score_BIC_from_cov". Pass "local_score_BIC_from_cov" to read whole numbers, such as counts, as numbers.

• lambda_value (float) – Multiplier on the BIC penalty of log(n_samples) per parameter: per parent for the linear-Gaussian BIC, per free parameter of the conditional probability table for the multinomial BIC. Larger values give sparser graphs; 1.0 is the standard BIC. Must be finite and non-negative. BDeu ignores it.

• seed (int | None) – Seed for the random order in which the search visits variables. None uses the ANDREY_SEED environment variable, or 0 when it is unset, so repeated calls return the same graph. Different seeds can reach different graphs.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data score: column labels (strings) or 0-based positions. None infers them: a DataFrame’s category, bool, and text columns, and every column of at most max_categories distinct whole numbers, named in a DataTypeWarning. [] marks none. The other scores take no column types and reject it. See Data types in the getting-started guide.

• max_categories (int) – The most distinct whole numbers a column may hold to be inferred as categorical, at least 2. With a continuous score, columns within it get a one-time DataTypeWarning that they may be category codes. With score_func=None, the columns it infers pick the default score.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. Each variable’s parents exclude the forbidden ones and include the required ones, and every causal order the search visits puts the cause of each required edge, and each earlier tier, first. The result directs every edge the knowledge decides, then Meek’s rules direct the edges that follow, as Knowledge describes.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "BOSS" and metadata["score"] the score of the DAG read off the final order; for the BIC, the sum over variables of n * log(residual variance) + lambda_value * log(n) * n_parents, for BDeu, -2 times the log of the structure prior times the marginal likelihood, and for the multinomial BIC, -2 times the log-likelihood plus the penalty; lower is better. Directed marks (->) are oriented edges, undirected (--) the ones the equivalence class leaves open. Under a mixed-data setting, metadata["categorical"] lists the categorical columns it read, by name or 0-based position.

Return type:

StructureOutput

Raises:

• NotImplementedError – If score_func names no available score.

• TypeError – If categorical is not a list of column names and positions, or max_categories is not an int. Also if knowledge is neither None nor a Knowledge.

• ValueError – If lambda_value is negative or non-finite, or if data or categorical does not suit the score, as for ges(). With the BIC score, also if a column repeats another up to scale and shift given more than two rows or, with more rows than columns, is a linear function of other columns. One message names every group of copies and says to keep one column per group. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When score_func is None and the columns are all categorical or mixed, so the search scores with BDeu or the conditional Gaussian BIC; when a mixed-data score infers the categorical columns, or finds none and so runs as the linear-Gaussian form; or when a continuous setting meets columns that look like category codes. Each message warns once per process.

See also

andrey.grasp

Greedy relaxation of the sparsest permutation algorithm [2], searching the same space of orders with different moves.

andrey.ges

Searches over equivalence classes for the same kind of data.

andrey.fci

Allows unmeasured common causes and selection bias.

andrey.direct_lingam

Orients every edge when the noise is non-Gaussian.

Notes

When to use it. BOSS suits the same data as ges(): continuous data with linear effects, Gaussian noise, and no unmeasured common causes. It was designed for hundreds of highly connected variables [1]. When every column is categorical, the default score is "local_score_BDeu", and on mixed data "local_score_cond_gauss"; "local_score_BIC_categorical" and "local_score_degen_gauss" also score those data.

When not to. When unmeasured common causes are plausible, use fci(). When the noise is non-Gaussian and every edge should be oriented, use direct_lingam().

Choosing seed. The search visits variables in a random order and stops when no move improves the order. Different seeds can therefore reach different graphs; comparing a few shows whether a result depends on the seed.

How it works. BOSS searches variable orders for the best DAG with the lowest total score [1]. Each variable uses a grow-shrink tree to cache its best parent set for any set of variables allowed to precede it. Each sweep moves every variable to the position that most improves the total score. Sweeps repeat until no move improves it. The final parent sets form a DAG, which is reduced to its CPDAG.

A move or parent change must improve the score by more than 1e-6, a fixed absolute tolerance shared with grasp(). Scores within this tolerance of the best count as ties. The earliest order position wins; among candidate parents, the column with the lowest tie rank wins. Under the linear-Gaussian BIC each comparison is decided as full local scores decide it: a fast update whose rounding-error bound could change the comparison is scored in full.

Column order. Each column gets a tie rank computed from its own data, as for hc(), and the search runs on the columns sorted by it: it starts from that order, and the seed shuffles the order in which it visits them. The column-order guide says what the rank keeps under a new column order and on another machine.

Categorical data. "local_score_BDeu" reads every column as categorical, each distinct value one category, and scores a variable given its parents with the BDeu score of ges(): a Dirichlet prior with an equivalent sample size of 1 and a structure prior of one expected parent per variable, on the scale of the BIC. "local_score_BIC_categorical" reads the columns the same way and scores a variable with the multinomial BIC of ges(): -2 times the log-likelihood of its counts in each configuration of its parents, plus lambda_value * log(n) for each free parameter. Under either score each column gets its tie rank as a categorical column, and each candidate parent set is scored by full local scores.

Mixed data. "local_score_degen_gauss" [3] and "local_score_cond_gauss" [4] score data with continuous and categorical columns, as ges() describes, and read the categorical columns as pc() does. Each candidate parent set is scored by full local scores, without the Schur update the BIC search uses.

Limitations. Measured on simulated linear-Gaussian data from sparse random graphs (mean degree 2), from 20 to 800 variables with ten rows per variable, at lambda_value=4.0: its median errors (SHD) stayed below the empty graph’s at every size, 24 against 800 at 800 variables. The median fit took 0.031 s at 20 variables and 292 s at 800, on 16 cores. On the 223-variable ANDES network with 1,000 simulated rows, at lambda_value=1.0, its median errors passed the empty graph’s (506 against 338). With fewer rows than columns the BIC score is unreliable. A parent set of n_samples - 1 or more variables fits its child exactly, so the search never selects one, nor any parent set that leaves its child a residual variance of at most 1e-12 times its variance. Against exact_search() on small simulated categorical networks, with BDeu or the multinomial BIC, BOSS stopped at a DAG that scores worse than the optimum on some draws; comparing a few seeds shows whether a result depends on the seed. Neither categorical score is benchmarked with BOSS yet. The mixed-data scores are experimental and not benchmarked yet; their tests run BOSS on 5 variables with 3,000 rows.

References

[1] (1,2)

Andrews, B., Ramsey, J., Sanchez Romero, R., Camchong, J., and Kummerfeld, E. (2023). Fast Scalable and Accurate Discovery of DAGs Using the Best Order Score Search and Grow Shrink Trees. NeurIPS 36. https://doi.org/10.52202/075280-2794

[2]

Lam, W.-Y., Andrews, B., and Ramsey, J. (2022). Greedy Relaxations of the Sparsest Permutation Algorithm. UAI, PMLR 180:1052-1062. https://proceedings.mlr.press/v180/lam22a.html

[3]

Andrews, B., Ramsey, J., and Cooper, G. F. (2019). Learning high-dimensional directed acyclic graphs with mixed data-types. Proceedings of Machine Learning Research 104, 4-21.

[4]

Andrews, B., Ramsey, J., and Cooper, G. F. (2018). Scoring Bayesian networks of mixed variables. International Journal of Data Science and Analytics 6(1), 3-18.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import boss
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = rng.standard_normal(500)
>>> z = x + y + 0.5 * rng.standard_normal(500)
>>> boss(pd.DataFrame({"x": x, "y": y, "z": z})).structure.oriented_edges()
[('x', 'z', 'directed'), ('y', 'z', 'directed')]