andrey.rsmax2

andrey.rsmax2(data, *, restrict='mmpc', maximize='hc', alpha=0.05, indep_test=None, score_func=None, penalty_discount=1.0, ebic_gamma=1.0, equivalent_sample_size=1.0, max_iter=None, tabu_length=10, max_tabu=None, categorical=None, seed=None, max_depth=None, max_parents=None, knowledge=None)

Learn a CPDAG from observational data with a restrict phase and a maximize phase.

The restrict phase learns the pairs that may be adjacent with conditional-independence tests; the maximize phase runs a score search over DAGs whose edges join only those pairs [1] [5]. A small candidate set makes the score search faster on many variables.

Parameters:

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

• restrict (str) – The restrict phase. "mmpc" keeps a pair when each variable is in the other’s MMPC parents and children [2], as mmhc() does. "iamb" and "inter_iamb" keep a pair when each variable is in the other’s Markov blanket from IAMB or Inter-IAMB [3] and no subset of the smaller blanket separates the two, the skeleton of iamb(structure="skeleton"), which drops the pairs of parents of a common child. "pc" keeps the skeleton of pc()’s adjacency search [4]. See Notes.

• maximize (str) – The maximize phase: hc() or tabu().

• alpha (float) – Significance level of the restrict phase’s tests. Smaller values keep fewer pairs.

• indep_test (str | None) – Conditional-independence test of the restrict phase, any of pc()’s. It must read the same kind of data as score_func. None picks from the declared column types, as for pc(). metadata["indep_test"] names the test the fit used: the one picked from the declared column types when this is None.

• score_func (str | None) – Local score of the maximize phase, any of hc()’s. None picks from the declared column types, as for ges(). metadata["score_func"] names the score the fit used: the one picked from the declared column types when this is None.

• penalty_discount (float) – Multiplier on the BIC penalty, as for hc(). BDeu ignores it.

• ebic_gamma (float) – Weight on the extended-BIC term of the linear-Gaussian BIC: each variable’s score adds 2 * ebic_gamma * log C(n_variables - 1, n_parents), which grows with the number of parent sets a variable can choose from. Larger values give sparser graphs; 0.0 is the standard BIC. Must be finite and non-negative. The other scores ignore it. See Notes.

• equivalent_sample_size (float) – The equivalent sample size of the BDeu score: how many rows the uniform prior over each variable’s conditional distributions is worth. A larger value weighs the prior more against the data, and the graph can change with it. Must be finite and positive. Only "bdeu" uses it: passing it with another score raises ValueError.

• max_iter (int | None) – Maximum moves of the maximize phase, at least 1; None sets no limit.

• tabu_length (int) – With maximize="tabu", the number of most recently visited DAGs a move may not return to, at least 1, as for tabu(); hill climbing ignores it.

• max_tabu (int | None) – With maximize="tabu", the non-improving moves allowed without a better DAG, from 1 to tabu_length; None uses tabu_length. Hill climbing ignores it.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data test and score, as for pc(). A continuous setting raises ValueError on a declared column.

• seed (int | None) – None, the default, is the configured seed, ANDREY_SEED or 0 when it is unset, not a fresh random one, so a fit repeats. Seed for the random draws of indep_test, which only "rcot", "cmiknn", and "gpdc" make; the other tests ignore it. With those tests, metadata["seed"] records the seed used.

• max_depth (int | None) – The maximum conditioning-set size of every conditional-independence test in the restrict phase, at least 0: 0 runs only unconditional tests, and None sets no limit. A smaller cap runs fewer tests and can keep pairs that a larger conditioning set would remove.

• max_parents (int | None) – A cap on every node’s parents, at least 0: the score search takes no move that gives a node more than max_parents parents, so no DAG the returned CPDAG represents does (every DAG of a class has the same largest number of parents). A cap can stop the search short of the lowest-scoring graph. None sets no cap.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. The maximize phase considers the pairs the restrict phase keeps and the pairs of the required edges, less the pairs the knowledge makes non-adjacent, and then honors the knowledge as hc() does. The restrict phase reads the data alone.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "rsmax2", metadata["score"] the score of the DAG found, computed as for andrey.ges() (lower is better). Under a mixed-data setting, metadata["categorical"] lists the categorical columns it read. With a randomized indep_test, metadata["seed"] is the seed used.

Return type:

StructureOutput

Raises:

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

• ValueError – If an argument has the wrong type, lies outside its range, or is not one of its choices; the message names the method, the argument, what it accepts, and the value. Also if score_func is not a score maximize takes, if score_func and indep_test read different kinds of data, if max_tabu is greater than tabu_length, or if data or categorical does not suit the test or score, as for pc() and hc(). Also if the required edges of knowledge give a node more than max_parents parents, or if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

• SearchLimitWarning – When the maximize phase takes a given max_iter moves.

• PerformanceWarning – With restrict="pc", when a conditioning pass needs at least ten million tests, as for pc().

• DataTypeWarning – As for pc(), when a mixed-data setting finds no declared categorical column or a setting reads columns that look like category codes as numbers; once per process.

See also

andrey.mmhc

restrict="mmpc" with maximize="hc".

Notes

When to use. The restrict-then-maximize search suits i.i.d. observational data with no unmeasured common causes, in which every independence follows from the graph (faithfulness) [2]. Its default Fisher-Z test and BIC score assume continuous data with linear effects and Gaussian noise; with categorical or mixed columns, the defaults follow the declared column types, as for pc() and ges(). It is a hybrid, a constraint-based skeleton from the tests restrict names followed by the score search maximize names over only the skeleton’s pairs, designed for many variables with few candidate neighbors each [1]. It returns the CPDAG of the DAG the maximize phase ends with. When unmeasured common causes are plausible, use gfci(); when the noise is non-Gaussian and every edge should be oriented, use direct_lingam().

How it works. The restrict phase follows the sparse candidate idea of [1]: tests narrow the pairs a score search must consider. Each choice finds every true adjacency when its tests are exact. The skeletons of MMPC, IAMB, Inter-IAMB, and PC then hold only the true adjacencies. A Markov blanket also holds the other parents of a common child, so IAMB and Inter-IAMB test each pair of their moral graph against subsets of the smaller blanket and drop the pairs that separate. The maximize phase starts from the empty graph, or from the required edges, with every pair outside the candidates forbidden in both directions, and returns the CPDAG of the DAG it ends with. mmpc(), iamb(), and inter_iamb() describe their searches; "pc" runs pc()’s adjacency search without orienting.

Ties go to the column whose rank, computed from the data, comes first, as for hc(), so the result does not depend on the column order unless two columns are identical.

Choosing the penalty. The maximize phase scores as hc() does: the linear-Gaussian BIC charges a variable penalty_discount * log(n_samples) per parent plus the extended-BIC term 2 * ebic_gamma * log C(n_variables - 1, n_parents) [6], and ges() describes the choice. ebic_gamma=0.0 with penalty_discount=1.0 is the standard BIC. With restrict="pc" and maximize="tabu", on simulated linear-Gaussian data with 20 to 400 variables and 200 to 4,000 rows, the default ebic_gamma=1.0 made fewer or as many errors (SHD) as the standard BIC on average in 63 of 65 settings of graph, variables, and rows, and more in two: 38.0 against 36.3 on the 46-variable ECOLI70 network with 1,000 simulated rows, where one of three draws moved by 7, and 39.2 against 38.8 with 50 variables, mean degree 4, scale="unitless", and 1,000 rows. With 200 variables and 2,000 rows it made 42 errors against 65 (200 true edges). The default restrict="mmpc" with maximize="hc" is mmhc()’s search, whose numbers are given there.

Missing values. A NaN in data raises ValueError: its score phase scores every variable over the same rows, and score terms over different rows do not add up to the score of one dataset. Drop the incomplete rows first, which suits values missing completely at random, or use pc() or fci(), whose tests delete rows per test.

Limitations. Measured on simulated linear-Gaussian data with 20, 50, and 100 variables, 1,000 and 5,000 rows, and mean degree 2 and 4, with Fisher-z tests at alpha 0.05 and the BIC, on one core, with restrict="mmpc" and "iamb" and maximize="hc". In every setting the mean CPDAG errors stay below the empty graph’s; the closest is at 20 variables and mean degree 4, with 29 to 31 pairs wrong for "mmpc" and 21 to 24 for "iamb", against 40 for the empty graph. A fit took 0.04 to 3.5 s with "mmpc" and 0.07 to 8.8 s with "iamb", the longest at 100 variables, 5,000 rows, and mean degree 4. An edge the restrict phase misses is never added, and the maximize phase stops at a local optimum. MMPC tests each candidate against every subset of the current members, so its cost grows exponentially with the size of a variable’s parents-and-children set. IAMB and Inter-IAMB condition on the whole current blanket, so each of their tests needs more rows as the blankets grow.

References

[1] (1,2,3)

Friedman, N., Nachman, I., and Pe’er, D. (1999). Learning Bayesian network structure from massive datasets: the “sparse candidate” algorithm. UAI, 206-215.

[2] (1,2)

Tsamardinos, I., Brown, L. E., and Aliferis, C. F. (2006). The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning 65(1), 31-78.

[3]

Tsamardinos, I., Aliferis, C. F., and Statnikov, A. (2003). Algorithms for large scale Markov blanket discovery. FLAIRS, 376-380.

[4]

Colombo, D., and Maathuis, M. H. (2014). Order-independent constraint-based causal structure learning. JMLR 15, 3921-3962.

[5]

Scutari, M. (2010). Learning Bayesian networks with the bnlearn R package. Journal of Statistical Software 35(3). https://doi.org/10.18637/jss.v035.i03

[6]

Chen, J., and Chen, Z. (2008). Extended Bayesian information criteria for model selection with large model spaces. Biometrika 95(3), 759-771.

Examples

>>> import warnings
>>> import numpy as np
>>> from andrey import rsmax2
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = x + 0.1 * rng.standard_normal(500)
>>> with warnings.catch_warnings():
...     warnings.simplefilter("ignore")
...     out = rsmax2(np.column_stack([x, y]), restrict="pc", maximize="tabu")
>>> out.structure.kind
'cpdag'