andrey.mmhc

andrey.mmhc(data, *, 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, categorical=None, seed=None, max_depth=None, max_parents=None, knowledge=None)

Learn a CPDAG from observational data with max-min hill climbing (MMHC).

MMHC [1] first learns a skeleton with the max-min parents and children search (MMPC), which uses conditional-independence tests, then runs hill climbing over DAGs whose edges join only the pairs of that skeleton. It is rsmax2() with restrict="mmpc" and maximize="hc".

Parameters:

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

• alpha (float) – Significance level of the conditional-independence test. Smaller values give sparser skeletons, so the score search considers fewer pairs.

• indep_test (str | None) – Conditional-independence test of the MMPC phase, any of pc()’s: "fisherz", "student_t", or "normal_scores" for continuous data, "chisq" or "gsq" for categorical data, "degen_gauss" or "cond_gauss" for mixed data, or a nonlinear test. 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 hill-climbing phase, any of hc()’s: "bic" (Gaussian BIC) for continuous data, "bdeu" or "discrete_bic" for categorical data, "degen_gauss" or "cond_gauss" for mixed data. 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 hill-climbing moves, at least 1; None sets no limit.

• 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 MMPC 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 score search considers the pairs of the MMPC skeleton and of the required edges, less the pairs the knowledge makes non-adjacent, and then honors the knowledge as hc() does. The MMPC phase reads the data alone.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "mmhc", 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 and indep_test read different kinds of data, 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 search takes a given max_iter moves.

• 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.rsmax2

Other restrict and maximize phases.

andrey.mmpc

The MMPC phase alone, for one variable or all.

Notes

When to use. MMHC suits i.i.d. observational data with no unmeasured common causes, in which every independence follows from the graph (faithfulness) [1]. 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 conditional-independence tests followed by a score-based search over only the skeleton’s pairs, designed for many variables with few parents and children each [1]. It returns the CPDAG of the DAG the search stops at. 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 MMPC phase estimates each variable’s parents and children: a forward phase adds, one at a time, the candidate whose weakest association with the variable, over every conditioning set drawn from the members so far, is the strongest, and a backward phase drops each member some subset of the others separates (mmpc()). A pair enters the skeleton when each variable is in the other’s set, the symmetry correction of [1]. The hill-climbing phase starts from the empty graph, or from the required edges, and adds, removes, or reverses one edge at a time while the score improves, with every pair outside the skeleton forbidden in both directions. The result is the CPDAG of the DAG it stops at.

Ties between candidates in the MMPC phase and between moves in the hill-climbing phase 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. The column-order guide lists what reordering the columns does to each method.

Choosing the penalty. The hill-climbing 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) [2], and ges() describes the choice. ebic_gamma=0.0 with penalty_discount=1.0 is the standard BIC. 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 64 of 65 settings of graph, variables, and rows, and 0.2 more in the last (7.0 against 6.8 with 20 variables, mean degree 2, scale="unitless", and 1,000 rows); with 200 variables and 2,000 rows it made 49 errors against 73 (200 true edges). The gain is smaller than hill climbing’s, since the MMPC phase has already left out most false edges.

Missing values. A NaN in data raises ValueError: its hill-climbing 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. 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 against 40 for the empty graph. A fit took 0.05 to 1.8 s, the longest at 100 variables, 5,000 rows, and mean degree 4. An edge the MMPC phase misses is never added, and the hill-climbing 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.

References

[1] (1,2,3,4)

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.

[2]

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 mmhc
>>> 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 = mmhc(np.column_stack([x, y]))
>>> out.structure.kind
'cpdag'