ReferenceDiscovery

Discovery

Causal discovery: learn causal structure from data. Import from ergodic.discovery.

Two ingredients drive the searches. A conditional-independence test (CITest) is a class constructed with data, which makes it the stateful oracle the constraint algorithms query: Linear (partial correlation), Discrete (G squared, natively categorical), Mixed (degenerate Gaussian for mixed columns), GCM (any regression Learner becomes a test), and DSeparationOracle as ground truth. A decomposable Score answers local(node, parents) for the score-based searches: BIC (Gaussian), DiscreteBIC, BDeu, and MixedBIC (degenerate Gaussian). Algorithms are thin consumers: PC wraps a test and returns a CPDAG, FCI wraps a test and returns a PAG (latent confounders allowed), MMPC wraps a test and returns an undirected skeleton screen, GES wraps a score and returns a CPDAG, AStar wraps a score and returns the exact optimum, and HillClimb is the fast greedy baseline. discover is the one-line facade, and bootstrap() reruns any method on resamples and returns edge frequencies (BootstrapResult), the first soft output. EdgeProbabilities holds a model's directed-edge probabilities (or bootstrap frequencies) over named variables, and decode() turns one into a DAG with an explicit, versioned decoder (DecodeResult). test_graph() runs the other direction: it tests the independences a candidate graph implies against data and combines them with Fisher's C (GraphTest). graph_posterior() turns edge probabilities, or data and a score (BGe for linear-Gaussian data), into a posterior over DAGs under that method (GraphSamples), exact by a subset recursion when it can be and by structure MCMC otherwise, and evaluate_graph() reads a candidate DAG, CPDAG, PAG or ADMG against it (GraphEvaluation); adjustment_probability() and average_effect() carry it to the effect level (AveragedEffect). Background knowledge (including a skeleton restriction through restrict_to) enters through DomainKnowledge.

Time series and panels get their own methods: PCMCI reads a TimeSeriesData or PanelData through a lagged design matrix (lagged_tabular()) and returns a window graph over nodes like "X[t-1]" (TimeSeriesDiscoveryResult), and granger() is the predictive baseline, a pair-level table on purpose rather than a graph.

Also exported from ergodic: DecodeResult, EdgeProbabilities, GraphTest, adjustment_probability, average_effect, decode, discover, evaluate_graph, granger, graph_posterior, test_graph.

BIC#

class BIC(data: TabularData | None = None, *, penalty: float = 1.0) -> None

The Gaussian BIC, decomposable and score-equivalent.

Nodes must be single numeric columns (continuous, datetime, or binary); scoring a wider categorical block raises TypeError lazily, at the first query against that node. Categorical parents are fine: their one-hot columns enter the regression and each column counts toward k.

BIC.penalty#

property BIC.penalty: float

The multiplier on the complexity term.

BIC.local_from_arrays#

def BIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The same local by direct least squares on (n, 1) y and (n, k) X.

Agrees with the covariance path to numerical precision; the tests hold the two to 1e-8.

BIC.is_bound#

property BIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

BIC.data#

property BIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

BIC.n#

property BIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

BIC.variables#

property BIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

BIC.bind#

def BIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

BIC.local#

def BIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

BIC.local_diff#

def BIC.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

BIC.enumerate_locals#

def BIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

BIC.graph_score#

def BIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

BIC.diagnostics#

def BIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

BIC.local_from_codes#

def BIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

BOSS#

class BOSS(
    score: Score | str = 'auto',
    *,
    use_bes: bool = True,
    restarts: int = 0,
    seed: int | None = None,
) -> None

Best order score search over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

BOSS.fit#

def BOSS.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Sweep permutations to convergence and return the class found.

DECODERS#

DECODERS: mappingproxy = mappingproxy({'greedy_v1': 'greedy pair-MAP: best-scoring edges first, skip any that closes a cycle', 'exact_map_v1': 'exact MAP of the factorized pair score, A* over topological orders', 'expected_count_v1': 'superseded by expected_count_v2: budget counts hard-forbidden mass', 'expected_count_v2': "the MAP order's best K edges, K the rounded expected edge count left open by the hard knowledge", 'soft_f1_v1': 'superseded by soft_f1_v2: budget counts hard-forbidden mass', 'soft_f1_v2': "the MAP order's best K edges, K from a label-free soft-F1 proxy over what the hard knowledge leaves open", 'pair_sampler_v1': "sample each pair's state, add edges in random order, drop cycle-closing ones", 'perturb_and_map_v1': 'Gumbel-perturb the pair states, then solve exact MAP per sample'})

The registered decoders: versioned id to a one-line summary (read-only).

FCI#

class FCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
) -> None

FCI over a conditional-independence test, returning a PAG.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the PAG face when the result is a valid complete PAG and the mixed face otherwise (noisy tests or knowledge-seeded marks). With a perfect oracle and no knowledge the output is exactly the canonical PAG: pag_of_mag(mag_of_dag(truth)).

FCI.fit#

def FCI.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run FCI and return the discovered graph with its provenance.

GCM#

class GCM(data: TabularData | None = None, *, learner: Any = None) -> None

The generalized covariance measure over a regression learner.

Regresses x and y on z with a fresh clone of the learner per endpoint, forms the elementwise residual product, and reports sqrt(n) * mean(R) / sd(R) against the standard normal, two-sided. Residuals are in-sample, as in the paper's basic form. With an empty conditioning set the endpoints are centered instead, so the statistic is a plain covariance test.

GCM.statistic#

def GCM.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

GCM.is_bound#

property GCM.is_bound: bool

Whether the test is bound to data and can answer queries.

GCM.data#

property GCM.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

GCM.n#

property GCM.n: int

Rows the oracle answers from (after dropping incomplete rows).

GCM.variables#

property GCM.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

GCM.bind#

def GCM.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

GCM.test#

def GCM.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

GCM.independent#

def GCM.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

GCM.diagnostics#

def GCM.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

GES#

class GES(
    score: Score | str = 'auto',
    *,
    max_degree: int | None = None,
) -> None

Greedy equivalence search over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

GES.fit#

def GES.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run the forward and backward phases and return the class found.

MMPC#

class MMPC(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = 3,
) -> None

The max-min parents-and-children screen over a conditional-independence test.

MMPC is a cheap skeleton screen to run before an expensive search: it returns the undirected adjacencies the data supports, and DomainKnowledge().restrict_to(result.graph) turns them into a support that PC, GES, or AStar then searches inside, with far fewer tests or score calls. Adjacencies are kept under the AND rule (each endpoint must sit in the other's parents-children set), the stricter of the two symmetry choices, so the screen errs toward sparsity.

Knowledge enters in three modes. Prune: a pair forbidden in both directions is never a candidate. Protect: required and confounded pairs enter every CPC untested, survive the backward sweep, and pass the AND rule unconditionally. Validate: the finished skeleton is checked at the adjacency level only, because a skeleton has no directions, so is_consistent (which demands required edges be directed) cannot apply; orientation constraints (tiers, single-direction forbids, roots, sinks) do not constrain an undirected adjacency and are left to the search that runs inside the support.

fit(data, knowledge) returns a DiscoveryResult whose graph holds one tail-tail edge per surviving pair, with the mixed face.

MMPC.fit#

def MMPC.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run MMPC and return the undirected skeleton with its provenance.

PC#

class PC(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
    colliders: str = 'standard',
) -> None

PC-stable over a conditional-independence test.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the oriented PDAG is a valid CPDAG and the mixed face otherwise.

PC.fit#

def PC.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run PC-stable and return the discovered graph with its provenance.

PCMCI#

class PCMCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_lag: int = 1,
    max_cond_size: int | None = None,
) -> None

PCMCI+ over a conditional-independence test, for time-series and panel data.

fit(data) takes TimeSeriesData or PanelData (panel windows stack per entity) and returns a TimeSeriesDiscoveryResult whose graph lives on the window nodes v[t-k].

fit(data, knowledge) reads a DomainKnowledge over the process names, not the window nodes, so a constraint on (x, y) applies at every lag, the contemporaneous slice included:

• A directional forbid x -> y means the link (x, y, lag) is never screened or MCI-tested at any lag and the lag-zero orientation x -> y is banned; the contemporaneous edge itself may survive when y -> x stays permitted (only the direction dies). A pair forbidden in both directions is an adjacency ban: the contemporaneous pair is removed without testing (and records no sepset) and both lagged directions are pruned. An explicit forbid(x, x) removes the self-memory link x[t-1] -> x[t].
• Tiers and order: a later variable never causes an earlier one at any lag, pruned exactly as a forbid; additionally a surviving contemporaneous edge between ordered processes is seeded earlier -> later before the collider phase, and a data-driven arrowhead that disagrees follows the orientation-conflict convention already in this file (the edge stays undirected and is counted in diagnostics["orientation_conflicts"]).
• Roots take no incoming cross-variable link at any lag and sinks no outgoing one; self-memory x[t-1] -> x[t] survives both, since a process can carry its own past without being caused by anything else. A skeleton restriction (restrict_to) bans pairs outside its support the same way, self-memory exempt.
• required (x, y) names no lag, so it cannot be forced; it is validated after the run instead, and fit raises ValueError when no compatible link (a lagged x -> y, an oriented lag-zero x -> y, or a still-undirected lag-zero edge) survived.
• confound has no analogue in a window graph (it carries no bidirected edges) and raises NotImplementedError.

diagnostics["n_knowledge_pruned"] counts what the knowledge took off the table: lagged links never screened, contemporaneous pairs removed outright, and lag-zero direction bans, one count each.

PCMCI.fit#

def PCMCI.fit(
    data: TimeSeriesData | PanelData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> TimeSeriesDiscoveryResult

Run PCMCI+ and return the window graph with its provenance.

knowledge constrains processes, so every constraint applies at every lag; the class docstring lists the mappings. confound raises here (a window graph has no bidirected edges), and a required edge that survived at no lag raises after the run.

RFCI#

class RFCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
) -> None

RFCI over a conditional-independence test: FCI's marks at PC's price.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the PAG face when the output happens to be a valid complete PAG and the mixed face with diagnostics["pag_valid"] = 0.0 otherwise. Unlike FCI, the mixed face is an expected honest outcome here even under a perfect oracle: the RFCI skeleton can keep edges only the possible-d-sep phase would remove, and the result then stands for more than one Markov equivalence class. Domain knowledge enters exactly as in FCI: prune (a pair forbidden in both directions is never tested, unless it is confounded), protect (required and confounded pairs are never tested or removed), orient through ancestry (the same seeding helpers), validate at the end.

RFCI.fit#

def RFCI.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run RFCI and return the discovered graph with its provenance.

AStar#

class AStar(
    score: Score | str = 'auto',
    *,
    max_parents: int | None = None,
    node_limit: int = 20,
) -> None

Exact structure search: the highest-BIC DAG, found with A*.

AStar.fit#

def AStar.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Search for a highest-scoring DAG and return it with its provenance.

AdjustmentGroup#

class AdjustmentGroup(
    adjustment_set: frozenset[str] | None,
    mass: float,
    mass_error: float,
    n_dags: int,
    value: float,
    std_error: float | None,
    estimate: EffectEstimate | None = None,
    bound: bool = False,
) -> None

The DAGs that share one adjustment set, and the effect estimated with it.

AveragedEffect#

class AveragedEffect(
    treatment: str,
    outcome: str,
    groups: tuple[AdjustmentGroup, ...],
    mean: float,
    mean_mc_error: float,
    variance: float | None,
    within_variance: float | None,
    between_variance: float,
    identified_mass: float,
    zero_effect_mass: float,
    unidentified_mass: float,
    truncated_mass: float,
    exact: bool,
    diagnostics: dict[str, Any],
) -> None

An effect averaged over a posterior over DAGs; see average_effect().

AveragedEffect.std#

property AveragedEffect.std: float | None

The mixture's standard deviation, sqrt(variance).

AveragedEffect.to_frame#

def AveragedEffect.to_frame() -> pd.DataFrame

The groups as a dataframe, one row per group.

Columns: the set, its mass, the mass's error, whether that error is the rule-of-three bound, the DAG count, the estimate and its standard error.

The set reads "(empty)" for no adjustment and "(no directed path)" for the zero-effect group.

BDeu#

class BDeu(data: TabularData | None = None, *, ess: float = 10.0) -> None

The Bayesian Dirichlet equivalent uniform score, decomposable.

Every column must be discrete (categorical, ordinal, or binary): the score runs on integer codes, so binding to data with a continuous column raises TypeError. Score-equivalent for any fixed ess, but the value of ess shapes the result; see the module docstring.

BDeu.ess#

property BDeu.ess: float

The equivalent sample size of the Dirichlet prior.

BDeu.local_from_codes#

def BDeu.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The Dirichlet log marginal likelihood from integer codes.

BDeu.is_bound#

property BDeu.is_bound: bool

Whether the score is bound to data and can be evaluated.

BDeu.data#

property BDeu.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

BDeu.n#

property BDeu.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

BDeu.variables#

property BDeu.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

BDeu.bind#

def BDeu.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

BDeu.local#

def BDeu.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

BDeu.local_diff#

def BDeu.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

BDeu.enumerate_locals#

def BDeu.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

BDeu.graph_score#

def BDeu.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

BDeu.diagnostics#

def BDeu.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

BDeu.local_from_arrays#

def BDeu.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

BGe#

class BGe(
    data: TabularData | None = None,
    *,
    alpha_mu: float = 1.0,
    alpha_w: float | None = None,
) -> None

The BGe score: the linear-Gaussian marginal likelihood, decomposable and equivalent.

Every node must be a single numeric column (continuous, datetime or binary); binding to data with a categorical column of more than two levels raises TypeError. marginal_likelihood is True: summed over a DAG the locals are the log marginal likelihood of the data under that DAG, so graph_posterior() reads it as an exact posterior under the prior, not an approximation. See the module docstring for the formula.

BGe.alpha_mu#

property BGe.alpha_mu: float

The prior's equivalent sample size for the mean.

BGe.alpha_w#

property BGe.alpha_w: float

The Wishart degrees of freedom (resolved at binding; 0 before).

BGe.t#

property BGe.t: float

The prior scale's diagonal, t = alpha_mu (alpha_w - n - 1) / (alpha_mu + 1).

BGe.local_from_arrays#

def BGe.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The same local from (N, 1) y and (N, l) parent columns X.

Uses the bound score's variable count n and prior, and the given rows' own scatter matrix, so it agrees with local() on the bound data.

BGe.is_bound#

property BGe.is_bound: bool

Whether the score is bound to data and can be evaluated.

BGe.data#

property BGe.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

BGe.n#

property BGe.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

BGe.variables#

property BGe.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

BGe.bind#

def BGe.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

BGe.local#

def BGe.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

BGe.local_diff#

def BGe.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

BGe.enumerate_locals#

def BGe.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

BGe.graph_score#

def BGe.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

BGe.diagnostics#

def BGe.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

BGe.local_from_codes#

def BGe.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

BootstrapResult#

class BootstrapResult(
    method: str,
    n_runs: int,
    n_failed: int,
    edges: Mapping[frozenset[str], float],
    marks: Mapping[str, float],
    variables: tuple[str, ...],
) -> None

Edge frequencies from rerunning one discovery method on bootstrap resamples.

edges maps an unordered node pair to the fraction of successful runs in which the pair was adjacent, whatever the marks. marks maps a canonical edge glyph (glyph(), for example "X --> Y" or "X <-> Y") to its frequency, so one pair can spread its mass over several glyphs; the glyph frequencies for a pair sum to its adjacency frequency. Both mappings hold only pairs and glyphs that appeared at least once, with keys in sorted order. Frequencies are over the n_runs successful replicates, not over what was attempted. A temporal method keys them on the window nodes, "X[t-1] --> Y[t]".

BootstrapResult.frequency#

def BootstrapResult.frequency(a: str, b: str) -> float

The adjacency frequency of the pair, 0.0 when it never appeared.

BootstrapResult.summary#

def BootstrapResult.summary(threshold: float = 0.5) -> MixedGraph

The stability graph: pairs at or above threshold, one edge each.

A qualifying pair gets its modal glyph, the orientation it wore most often (ties broken by the sorted glyph string, so the result is deterministic). The graph keeps the mixed kind: a vote over graphs carries no equivalence-class guarantee, even when every replicate was a CPDAG. Nodes are all the data's variables, so a variable that never earned a stable edge shows up isolated. threshold must be in (0, 1].

CIResult#

class CIResult(
    statistic: float,
    p_value: float,
    dof: float | None = None,
) -> None

One conditional-independence query: the statistic, its p-value, and dof.

CITest#

class CITest(data: TabularData | None = None) -> None

Base for conditional-independence tests; see the module docstring.

Subclasses set name and needs, may take options in __init__ (set them before calling super().__init__(data), since binding may read them), and implement statistic. The optional _prepare hook runs once after binding, for sufficient statistics.

CITest.is_bound#

property CITest.is_bound: bool

Whether the test is bound to data and can answer queries.

CITest.data#

property CITest.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

CITest.n#

property CITest.n: int

Rows the oracle answers from (after dropping incomplete rows).

CITest.variables#

property CITest.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

CITest.bind#

def CITest.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

CITest.test#

def CITest.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

CITest.independent#

def CITest.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

CITest.diagnostics#

def CITest.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

CITest.statistic#

def CITest.statistic(*args: Any) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

DSeparationOracle#

class DSeparationOracle(graph: MixedGraph) -> None

The ground-truth oracle: answers queries by m-separation in a known graph.

Quacks like any bound test (p-value 1.0 when separated, 0.0 when connected), so every constraint algorithm can be checked for exactness and run in teaching mode. Needs no data; bind returns self.

Latent and selection nodes are not query-able variables: the oracle exposes the observed margin, which is what a constraint algorithm searching over data would see. Separation is still computed on the full graph, latents included, and every query conditions on the selection nodes, because data sampled under selection only ever shows the distribution given them.

DSeparationOracle.is_bound#

property DSeparationOracle.is_bound: bool

Always True: the oracle answers from its graph and needs no data.

DSeparationOracle.graph#

property DSeparationOracle.graph: MixedGraph

The graph whose d-separations (m-separations) the oracle reports.

DSeparationOracle.bind#

def DSeparationOracle.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

DSeparationOracle.data#

property DSeparationOracle.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

DSeparationOracle.n#

property DSeparationOracle.n: int

Rows the oracle answers from (after dropping incomplete rows).

DSeparationOracle.variables#

property DSeparationOracle.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

DSeparationOracle.test#

def DSeparationOracle.test(
    x: str,
    y: str,
    given: Iterable[str] = (),
) -> CIResult

The cached conditional-independence query, symmetric in x and y.

DSeparationOracle.independent#

def DSeparationOracle.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

DSeparationOracle.diagnostics#

def DSeparationOracle.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

DSeparationOracle.statistic#

def DSeparationOracle.statistic(
    *args: Any,
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

DecodeResult#

class DecodeResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
    decoder: str = '',
    candidates: tuple[MixedGraph, ...] = (),
    probabilities: EdgeProbabilities | None = None,
) -> None

A DAG (or ADMG) decoded from edge probabilities, with how it was decoded.

A DiscoveryResult whose graph is the decoded graph: the single decoded graph for a deterministic decoder, the consensus (medoid) of the candidates for a sampler. method is "decode:<decoder id>" (a Welt run decoded locally keeps "welt:<model>"). Its kind is dag, or admg when the knowledge holds confound pairs: each such pair gets a <-> b on top of the decoded directed edges. When the decoder also chose a directed edge on that pair (a bow), the confounder is written as an explicit latent node L__a__b (in graph.latent, with an edge into each end), since a graph holds one edge per pair; identification reads both forms alike.

diagnostics always labels the graph: decoded is True, native is False, and it records decoder, parameters, seed, samples, rejected_cycles (nominated edges skipped because they would close a directed cycle), order_skipped (skipped only because they would break the order knowledge) and search_expansions (each one per candidate), threshold_source, edge_budget, objective (the factorized pair log score of graph's directed edges on the input's pair states), knowledge_objective (the same score on the knowledge-weighted pair states the decoder read; equal to objective without soft statements), n_edges (directed edges), graph_kind, confounded_pairs, latent_confounders (latent node to pair, bows only), source_kind, source_provenance, knowledge_applied, n_hard_statements, n_soft_statements and note ("not a graph posterior; orientation not certified").

With knowledge it also holds the audit:

• knowledge_audit: one entry per statement (statements(): hard ones one per stored item, soft ones as made), a dict with statement (readable text), kind, confidence (None for hard), probability (the model's chance the whole statement holds, from the input's pair states with pairs taken as independent), min_pair_probability (the smallest chance among the pairs it covers that the statement's part on that pair holds), conflict_pairs (the pairs whose part gets less than 0.1, weakest first), holds (in graph) and holds_without_knowledge (in a decode with the same settings and seed but no knowledge). A confounder, which no edge state expresses, has None probabilities, no conflict pairs and the reason in note. holds reads each statement's own meaning: an order or tier fails when a later variable is an ancestor of an earlier one, through any directed path. The probabilities read direct edges only (an order's part on a pair is "no direct edge against it"), which note says for order and tier entries.
• knowledge_violations and knowledge_violations_without_knowledge: what violations() reports for the hard knowledge as a whole on graph (always empty: the decode enforces it) and on the knowledge-free decode's graph (None with audit=False), leaving confound pairs out: the knowledge-free decode is a DAG and never holds one, and their audit entries already say so. It reads order through the transitive closure, so it catches what no single statement's entry shows, such as c -> a under order(a, b) and order(b, c).
• knowledge_cost: the knowledge-free decode's objective minus this one's. Non-negative for exact_map; the other decoders are not exact on the score, so it can be negative.
• knowledge_conflicts: the audit entries with at least one conflict pair, a pair the model gives less than 0.1 for the statement's part on it (for tiers, an edge against the order above 0.9; for root(x), some P(v -> x) above 0.9). Any such entry also raises one KnowledgeConflictWarning. The whole-statement probability does not decide it: a statement over many pairs can be unlikely as a whole with no pair surprising.

Without knowledge, or with audit=False, these are (), None and ().

Decoder#

class Decoder(
    method: str = 'greedy',
    threshold: float | None = None,
    samples: int = 1,
    seed: int | None = None,
    max_nodes: int = 20,
    audit: bool = True,
) -> None

A decoder and its settings, checked up front and reusable.

Decoder("greedy", threshold=0.3).decode(probabilities) is decode(probabilities, "greedy", threshold=0.3). Pass one to discover() as decoder= to fix every setting before a paid run.

Decoder.stochastic#

property Decoder.stochastic: bool

Whether this decoder samples candidates (and so needs a seed).

Decoder.decode#

def Decoder.decode(
    probabilities: EdgeProbabilities,
    *,
    knowledge: DomainKnowledge | None = None,
) -> DecodeResult

Decode probabilities with these settings; see decode().

DirectLiNGAM#

class DirectLiNGAM(prune: str = 'bic') -> None

DirectLiNGAM over continuous tabular data.

fit(data, knowledge) returns a DiscoveryResult whose graph is a DAG: with independent non-Gaussian noises the whole orientation is identified, so there is no equivalence class to report. The class consumes neither a CI test nor a decomposable score (ingredient is None); the data is its only ingredient.

DirectLiNGAM.fit#

def DirectLiNGAM.fit(
    data: TabularData,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Find the causal order, prune the parents, and return the DAG.

DiscoveryResult#

class DiscoveryResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
) -> None

What a discovery run found, and how.

graph is the headline. test or score names the ingredient the algorithm consumed (one of them is None). diagnostics carries the run counters (tests by conditioning size, cache hits, score evaluations); a welt:<model> run carries its Welt provenance there instead (ids, versions and the native edge-score dataframe, see discover()). sepsets maps an unordered removed pair to the set that separated it, when a constraint method produced one.

Discrete#

class Discrete(data: TabularData | None = None) -> None

The G squared conditional-independence test on contingency tables.

Degrees of freedom are adjusted for sparsity, bnlearn-style: each nonempty stratum of the conditioning configuration contributes (card_x - 1) * (card_y - 1), and configurations with no observed rows contribute nothing. The unadjusted product over the full configuration space would overstate the dof badly on sparse tables and make the test blind. When the dof comes out zero (a single-level endpoint, or no rows), the query is vacuous and the result is (0.0, 1.0, 0.0).

Discrete.statistic#

def Discrete.statistic(
    x: npt.NDArray[np.int64],
    y: npt.NDArray[np.int64],
    z: tuple[npt.NDArray[np.int64], ...],
    cards: tuple[int, int, tuple[int, ...]],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Discrete.is_bound#

property Discrete.is_bound: bool

Whether the test is bound to data and can answer queries.

Discrete.data#

property Discrete.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Discrete.n#

property Discrete.n: int

Rows the oracle answers from (after dropping incomplete rows).

Discrete.variables#

property Discrete.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Discrete.bind#

def Discrete.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Discrete.test#

def Discrete.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Discrete.independent#

def Discrete.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Discrete.diagnostics#

def Discrete.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

DiscreteBIC#

class DiscreteBIC(
    data: TabularData | None = None,
    *,
    penalty: float = 1.0,
) -> None

The multinomial BIC, decomposable and score-equivalent.

Every column must be discrete (categorical, ordinal, or binary): the score runs on integer codes, so binding to data with a continuous column raises TypeError. Parameter counts agree across a Markov equivalence class, so Markov-equivalent DAGs get equal totals, exactly as with the Gaussian BIC.

DiscreteBIC.penalty#

property DiscreteBIC.penalty: float

The multiplier on the complexity term.

DiscreteBIC.local_from_codes#

def DiscreteBIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The multinomial local from integer codes; see the module docstring.

DiscreteBIC.is_bound#

property DiscreteBIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

DiscreteBIC.data#

property DiscreteBIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

DiscreteBIC.n#

property DiscreteBIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

DiscreteBIC.variables#

property DiscreteBIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

DiscreteBIC.bind#

def DiscreteBIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

DiscreteBIC.local#

def DiscreteBIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

DiscreteBIC.local_diff#

def DiscreteBIC.local_diff(
    node: str,
    parents: Iterable[str],
    add: str,
) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

DiscreteBIC.enumerate_locals#

def DiscreteBIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

DiscreteBIC.graph_score#

def DiscreteBIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

DiscreteBIC.diagnostics#

def DiscreteBIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

DiscreteBIC.local_from_arrays#

def DiscreteBIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

EdgeProbabilities#

class EdgeProbabilities(
    variables: tuple[str, ...],
    values: npt.NDArray[np.float64],
    kind: ProbabilityKind = 'independent',
    calibrated: bool = False,
    provenance: Mapping[str, Any] = <factory>,
) -> None

A matrix of directed-edge probabilities over named variables.

values[i, j] is the number for an edge from variables[i] (the source row) to variables[j] (the target column). The diagonal is zero, every entry is finite and in [0, 1], and the stored matrix is a read-only copy, so the caller's array is never shared or changed.

kind says how to read the numbers:

• "independent": one marginal per directed edge, each on its own. Both directions of a pair may be high at once (AVICI and CDFM return these).
• "pair_state": each unordered pair has three exclusive states (no edge, i -> j, j -> i) and the two directed cells are two of them, so P[i, j] + P[j, i] <= 1; the rest is the no-edge mass.
• "frequency": how often each directed edge appeared over resampled runs (see from_bootstrap()). One edge per pair per run, so the pair simplex holds here too.

calibrated is False unless the source states the numbers match observed edge frequencies; ergodic never sets it on its own. provenance is free-form, read-only metadata (model, version, job id, threshold, ...).

None of this is a posterior over graphs. The marginals do not account for acyclicity or for dependence between edges; a decoded graph built from them is a modelling choice, not the model's claim.

EdgeProbabilities.to_pair_states#

def EdgeProbabilities.to_pair_states() -> EdgeProbabilities

The pair-state reading: three exclusive states per unordered pair.

For "independent" marginals p = P[i, j] and q = P[j, i], the two edges are treated as independent coin flips and conditioned on "not both" (a DAG cannot hold a two-cycle): no edge (1 - p)(1 - q), i -> j p(1 - q), j -> i q(1 - p), renormalized. When both are exactly 1 the limit is undecided, so the pair becomes (0, 0.5, 0.5). This keeps each model's direction ranking and never looks at a true graph. It is the conversion Ergodic's research decoders were evaluated with.

"pair_state" values are returned unchanged (the same object). "frequency" values already satisfy the pair simplex (one edge per pair per run), so they are read as pair states as they are: a new object of kind "pair_state" with the same numbers.

Returns: Edge probabilities of kind "pair_state", with the conversion recorded in provenance["pair_state_conversion"].

EdgeProbabilities.existence#

def EdgeProbabilities.existence() -> npt.NDArray[np.float64]

Skeleton probabilities: the chance each pair is adjacent, either way.

Read from the pair-state view (to_pair_states()), so it is P[i, j] + P[j, i] there: for "independent" marginals that is (p(1 - q) + q(1 - p)) / (1 - pq). The matrix is symmetric with a zero diagonal and read-only.

EdgeProbabilities.threshold#

def EdgeProbabilities.threshold(t: float = 0.5) -> MixedGraph

The directed edges whose value is above t, as a diagnostic graph.

Keeps every edge with values[i, j] > t (strictly, as AVICI and CDFM threshold natively). The result is a mixed-kind graph over every variable and may hold directed cycles. A pair whose two directions both pass cannot be drawn as two edges (one edge per pair), so it shows as an undirected X --- Y edge: adjacent, direction unresolved. This is not a DAG and not a decoded graph; use decode() for one.

Returns: A mixed MixedGraph.

EdgeProbabilities.to_frame#

def EdgeProbabilities.to_frame() -> pd.DataFrame

The matrix as a dataframe: index source, columns target.

EdgeProbabilities.from_frame#

classmethod EdgeProbabilities.from_frame(
    frame: pd.DataFrame,
    *,
    kind: ProbabilityKind = 'independent',
    calibrated: bool = False,
    provenance: Mapping[str, Any] | None = None,
) -> EdgeProbabilities

Build from a square dataframe: rows the source, columns the target.

The index gives the variable order; the columns must hold the same names (in any order) and are reordered to match. Labels are read as strings.

EdgeProbabilities.from_bootstrap#

classmethod EdgeProbabilities.from_bootstrap(
    result: BootstrapResult,
    *,
    undirected: Literal['split', 'drop'] = 'split',
) -> EdgeProbabilities

Directed-edge frequencies from a bootstrap() run.

Each canonical glyph frequency lands on its directed cell: "X --> Y" on [X, Y] and "X <-- Y" on [Y, X]. An undirected X --- Y (a CPDAG edge, "one of the two directions") is split half and half under undirected="split" and dropped under "drop". Bidirected and circle-marked glyphs say something other than "this direction" and are always dropped; the dropped mass is recorded in provenance["dropped_frequency"]. Frequencies measure resampling stability, not probability.

Returns: Edge probabilities of kind "frequency".

FeatureProbability#

class FeatureProbability(
    probability: float,
    mc_error: float,
    exact: bool,
    bound: bool = False,
) -> None

A posterior probability with its Monte Carlo standard error.

When every draw agrees (a feature that holds in all draws or in none), the binomial error would be 0 although the value is only estimated. mc_error is then the rule-of-three bound 3 / N for N draws, a one-sided 95% bound on the mass the draws never saw, and bound is True. For MCMC the draws are correlated, so that bound is optimistic.

GRaSP#

class GRaSP(score: Score | str = 'auto', *, depth: int = 3) -> None

Greedy relaxations of the sparsest permutation over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

GRaSP.fit#

def GRaSP.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Walk tucks to convergence and return the class found.

Glasso#

class Glasso(
    lam: float | None = None,
    *,
    gamma: float = 0.5,
    n_lambdas: int = 100,
    tol: float = 0.001,
    max_iter: int = 100,
) -> None

The graphical-lasso screen: an undirected skeleton from a sparse precision.

Like MMPC, this is a cheap screen to run before an expensive search: DomainKnowledge().restrict_to(result.graph) turns the recovered support into a constraint that PC, GES, or AStar then searches inside, with far fewer tests or score calls. The honesty note: a precision zero is a conditional independence under Gaussianity, so the graph is a skeleton source, not causal structure; it says which adjacencies the data supports, never which way they point.

The solve is FHT 2008 block coordinate descent on the correlation matrix (the graph, not the scale, is the deliverable), with the diagonal unpenalized (the scikit-learn convention; the paper's differs, see the module docstring) and the shrunk initialization that keeps every column subproblem positive definite.

Knowledge enters through the penalty, and only its adjacency content binds, because the screen is undirected. A pair forbidden in both directions (an adjacency ban, including pairs outside a restrict_to support) is a hard zero, excluded from the updates entirely, unless it is confounded and not under forbid_adjacency() (a hidden common cause, neither causing the other): that pair is left to the data like any other. A required pair (either direction) is unpenalized, so it survives every lambda the data gives it any conditional association at. Directional knowledge (tiers, one-way forbids, roots, sinks) does not constrain an undirected adjacency and is left to the search that runs inside the support.

fit(data, knowledge) needs a TabularData whose variables are all continuous (a Gaussian model) and returns a DiscoveryResult whose graph holds one tail-tail edge per nonzero precision pair, on the mixed face. The chosen penalty lives in the diagnostics, with the EBIC value when the path ran.

Glasso.fit#

def Glasso.fit(
    data: TabularData,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run the screen and return the undirected skeleton with its provenance.

GrangerResult#

class GrangerResult(
    pvalues: Mapping[tuple[str, str], float],
    significant: frozenset[tuple[str, str]],
    max_lag: int,
    alpha: float,
    n: int,
) -> None

Granger tests for every ordered pair: a predictive table, not a graph.

pvalues maps (cause, effect) to the likelihood-ratio p-value; significant holds the pairs with p <= alpha. n counts the stacked window rows the regressions used.

GraphClaim#

class GraphClaim(
    claim: str,
    kind: str,
    variables: tuple[str, ...],
    probability: float | None,
    mc_error: float | None,
    exact: bool,
    note: str | None = None,
    bound: bool = False,
) -> None

One statement a candidate graph makes, and its posterior probability under the method.

GraphEvaluation#

class GraphEvaluation(
    candidate: MixedGraph,
    candidate_kind: str,
    method: str,
    log_target: float | None,
    best_graph: MixedGraph,
    best_log_target: float,
    best_source: str,
    log_ratio_to_best: float | None,
    credible_level: float | None,
    credible_level_error: float | None,
    mass_above: float | None,
    claims: tuple[GraphClaim, ...],
    class_probability: FeatureProbability | None,
    expected_shd: float | None,
    expected_shd_error: float | None,
    expected_shd_between_draws: float,
    p_acyclic: float | None,
    p_acyclic_error: float | None,
    not_assessable: tuple[str, ...],
    samples: GraphSamples,
    diagnostics: dict[str, Any],
) -> None

A candidate graph read against a posterior over DAGs; see evaluate_graph().

GraphEvaluation.claims_frame#

def GraphEvaluation.claims_frame() -> pd.DataFrame

The claims as a dataframe, one row each, weakest assessable claim first.

bound marks the rows whose mc_error is the rule-of-three bound.

GraphSamples#

class GraphSamples(
    variables: tuple[str, ...],
    method: str,
    target: str,
    adjacency: BoolArray,
    weights: FloatArray,
    log_targets: FloatArray,
    counts: IntArray,
    draws: IntArray,
    seed: int | None,
    knowledge: DomainKnowledge | None,
    probabilities: EdgeProbabilities | None,
    score: Score | None,
    diagnostics: dict[str, Any],
    _target: GraphTarget,
    _edges: FloatArray,
    _edge_errors: FloatArray,
    _best: tuple[BoolArray, float, str],
    _log_z: float | None,
    _conditional: FloatArray | None = None,
) -> None

Distinct DAGs with their posterior weights under one method.

adjacency[k] is a DAG over variables (rows the source), and weights[k] its weight: its exact posterior probability when every DAG was enumerated, otherwise its share of the draws (exact draws from the subset recursion, or MCMC draws after warmup). counts holds the draws per graph and draws the index of each kept draw, one row per chain (empty when enumerated). log_targets[k] is the unnormalized log target.

Build one with graph_posterior(). Every probability it returns is a posterior probability under the method, with a Monte Carlo error unless it is exact; see graph_posterior() for what "the method" is.

GraphSamples.exhaustive#

property GraphSamples.exhaustive: bool

Whether every DAG is listed with its exact probability (enumeration).

GraphSamples.n_draws#

property GraphSamples.n_draws: int

Kept draws over all chains (0 when enumerated).

GraphSamples.log_normalizer#

property GraphSamples.log_normalizer: float | None

log Z of the target when computed exactly, else None.

GraphSamples.exact_probabilities#

property GraphSamples.exact_probabilities: FloatArray | None

Each listed DAG's exact posterior probability, exp(log target - log Z).

None for MCMC, where Z is unknown.

GraphSamples.graph#

def GraphSamples.graph(index: int) -> MixedGraph

The index-th distinct DAG as a graph.

GraphSamples.top#

def GraphSamples.top(k: int = 5) -> list[tuple[MixedGraph, float]]

The k heaviest DAGs and their weights, heaviest first (ties by log target).

GraphSamples.best_graph#

property GraphSamples.best_graph: MixedGraph

The best DAG found: the exact MAP when exact, else the best of the searches.

GraphSamples.best_log_target#

property GraphSamples.best_log_target: float

The best DAG's log target.

GraphSamples.log_target#

def GraphSamples.log_target(graph: MixedGraph) -> float

The unnormalized log target of a DAG over variables; -inf off the support.

For edge probabilities it is the pair log score plus the knowledge's log prior; for a score, the score plus the log prior. Only differences matter: exp(log_target(G1) - log_target(G2)) is the posterior ratio.

GraphSamples.edge_probabilities#

def GraphSamples.edge_probabilities() -> EdgeProbabilities

Posterior edge marginals, as pair states.

values[i, j] is the posterior probability under the method of the edge variables[i] -> variables[j]. A DAG holds at most one direction per pair, so the two directions and "no edge" are exclusive and the result has kind "pair_state" ("frequency" means bootstrap stability, which this is not). The marginals are exact for the exact computation; for MCMC on edge probabilities they are Rao-Blackwellized (each draw contributes the pair's conditional state probabilities given the rest of the graph), else draw frequencies. provenance records the method, the estimator and the largest Monte Carlo error; edge_errors() has them all.

GraphSamples.edge_errors#

def GraphSamples.edge_errors() -> FloatArray

The Monte Carlo standard error of each edge marginal; zeros when exact.

GraphSamples.probability#

def GraphSamples.probability(
    feature: Callable[[MixedGraph], bool],
) -> FeatureProbability

The posterior probability that feature(G) holds, with its Monte Carlo error.

feature is any predicate on a DAG (a MixedGraph over variables), called once per distinct DAG. Exact (mc_error 0) when every DAG was enumerated. For exact draws the error is the binomial one, sqrt(p (1 - p) / draws); for MCMC it uses the feature's effective sample size across chains.

GraphTest#

class GraphTest(
    statistic: float,
    p_value: float,
    dof: int,
    k: int,
    alpha: float,
    test: str,
    kind: str,
    n: int,
    implications: tuple[ImpliedIndependence, ...],
    untestable: tuple[tuple[str, str], ...] = (),
    vacuous: tuple[ImpliedIndependence, ...] = (),
    notes: tuple[str, ...] = (),
) -> None

The result of test_graph(): Fisher's C over the basis set.

p_value is a frequentist p-value for "the data are compatible with the candidate's implied independences". A small value is evidence against the graph. A large one is not evidence for it: every Markov-equivalent graph gets the same answer, and with few rows any graph passes.

GraphTest.rejected#

property GraphTest.rejected: bool

Whether Fisher's combined p-value is at or below alpha.

This is the level-alpha test of the graph. The Holm decisions on the claims say which claims fail; rejecting the graph when either one rejects would spend more than alpha.

GraphTest.n_rejected#

property GraphTest.n_rejected: int

How many claims the Holm procedure rejects at alpha.

GraphTest.to_dict#

def GraphTest.to_dict() -> dict[str, Any]

The result as a JSON-ready dict.

GraphTest.to_frame#

def GraphTest.to_frame() -> pd.DataFrame

The claims as a dataframe, one row each, in the implications order.

HillClimb#

class HillClimb(
    score: Score | str = 'auto',
    *,
    tabu: int = 10,
    max_iter: int = 10000,
    max_parents: int | None = None,
) -> None

Greedy DAG search: add, remove, and reverse moves with a tabu ring.

Steepest ascent from the empty (or knowledge-seeded) graph, returning the best state seen. n_iterations in the diagnostics counts applied moves, so it always equals n_adds + n_removes + n_reverses.

HillClimb.fit#

def HillClimb.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Climb to a high-scoring DAG and return it with its provenance.

ImpliedIndependence#

class ImpliedIndependence(
    x: str,
    y: str,
    given: tuple[str, ...],
    rule: str,
    statistic: float,
    p_value: float,
    p_holm: float,
    rejected: bool,
    dof: float | None = None,
) -> None

One claim of the basis set and its test on the data.

ImpliedIndependence.to_dict#

def ImpliedIndependence.to_dict() -> dict[str, Any]

The claim as a JSON-ready dict.

Linear#

class Linear(data: TabularData | None = None) -> None

Partial correlation through the Fisher z transform.

Residualizes x and y on [1, z], computes the correlation r of the residuals, and reports sqrt(n - dz - 3) * atanh(r) against the standard normal, two-sided. The sign of r survives in the statistic; the p-value does not depend on it.

Linear.statistic#

def Linear.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Linear.is_bound#

property Linear.is_bound: bool

Whether the test is bound to data and can answer queries.

Linear.data#

property Linear.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Linear.n#

property Linear.n: int

Rows the oracle answers from (after dropping incomplete rows).

Linear.variables#

property Linear.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Linear.bind#

def Linear.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Linear.test#

def Linear.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Linear.independent#

def Linear.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Linear.diagnostics#

def Linear.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

Mixed#

class Mixed(data: TabularData | None = None) -> None

Partial canonical correlations with Bartlett's chi-square.

Residualizes the x and y blocks on [1, z], takes the reduced QR of each residual block, and reads the canonical correlations rho_i off the singular values of Qx.T @ Qy (clipped below 1). Bartlett's statistic -(n - dz - (dx + dy + 3) / 2) * sum(log(1 - rho_i^2)) is referred to a chi-square with dx * dy degrees of freedom.

With scalar endpoints (dx == dy == 1) the single canonical correlation is the absolute partial correlation, so Mixed agrees with Linear asymptotically. A zero-width block (a single-level categorical) carries no information; that query is vacuous and returns (0.0, 1.0, 0.0).

Mixed.statistic#

def Mixed.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Mixed.is_bound#

property Mixed.is_bound: bool

Whether the test is bound to data and can answer queries.

Mixed.data#

property Mixed.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Mixed.n#

property Mixed.n: int

Rows the oracle answers from (after dropping incomplete rows).

Mixed.variables#

property Mixed.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Mixed.bind#

def Mixed.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Mixed.test#

def Mixed.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Mixed.independent#

def Mixed.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Mixed.diagnostics#

def Mixed.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

MixedBIC#

class MixedBIC(
    data: TabularData | None = None,
    *,
    penalty: float = 1.0,
) -> None

The degenerate-Gaussian BIC, decomposable and score-equivalent.

Any column kind is a valid node: a continuous node is its single column and a categorical node is its one-hot block (first level dropped), scored jointly. local = -(n/2) log det Sigma_{Y|X} - penalty * p/2 * log n with p = dy dx + dy (dy + 1) / 2 for a dy-column node and dx total parent columns.

MixedBIC.local_from_arrays#

def MixedBIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The same local by direct least squares on (n, dy) y and (n, k) X.

Agrees with the covariance path to numerical precision; the tests hold the two to 1e-8.

MixedBIC.penalty#

property MixedBIC.penalty: float

The multiplier on the complexity term.

MixedBIC.is_bound#

property MixedBIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

MixedBIC.data#

property MixedBIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

MixedBIC.n#

property MixedBIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

MixedBIC.variables#

property MixedBIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

MixedBIC.bind#

def MixedBIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

MixedBIC.local#

def MixedBIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

MixedBIC.local_diff#

def MixedBIC.local_diff(
    node: str,
    parents: Iterable[str],
    add: str,
) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

MixedBIC.enumerate_locals#

def MixedBIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

MixedBIC.graph_score#

def MixedBIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

MixedBIC.diagnostics#

def MixedBIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

MixedBIC.local_from_codes#

def MixedBIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

PosteriorConvergenceWarning#

The structure MCMC chains disagree or mix too slowly for the reported numbers.

Score#

class Score(data: TabularData | None = None) -> None

Base for decomposable scores; see the module docstring.

Subclasses set name, needs, and the property flags, may take options in __init__ (set them before calling super().__init__(data)), and implement the method matching needs: local_from_arrays for "numeric", local_from_codes for "codes". _prepare runs once after binding for sufficient statistics; a native score may also override _local_value to skip array assembly entirely (the covariance path in BIC).

Score.is_bound#

property Score.is_bound: bool

Whether the score is bound to data and can be evaluated.

Score.data#

property Score.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

Score.n#

property Score.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

Score.variables#

property Score.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Score.bind#

def Score.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

Score.local#

def Score.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

Score.local_diff#

def Score.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

Score.enumerate_locals#

def Score.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

Score.graph_score#

def Score.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

Score.diagnostics#

def Score.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

Score.local_from_arrays#

def Score.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

Score.local_from_codes#

def Score.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

TimeSeriesDiscoveryResult#

class TimeSeriesDiscoveryResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
    max_lag: int = 0,
    links: Mapping[tuple[str, str, int], float] = <factory>,
) -> None

A DiscoveryResult over window nodes, with the per-link p-values.

links maps (source, target, lag) to the p-value of the test that kept the link. A lag of one or more is the lagged edge source[t-lag] -> target[t]. A lag of zero is a contemporaneous edge: an oriented edge is keyed (source, target), an undirected edge is keyed once with the pair in sorted order.

TimeSeriesDiscoveryResult.summary_edges#

def TimeSeriesDiscoveryResult.summary_edges() -> frozenset[tuple[str, str]]

Variable-level pairs implied by the links.

(x, y) is present when any lagged link x -> y exists or a contemporaneous edge x -> y is oriented in the window graph. An undirected contemporaneous pair contributes both orders, and a lagged self-link (autocorrelation) contributes (x, x).

adjustment_probability#

def adjustment_probability(
    samples: GraphSamples,
    treatment: str | Iterable[str],
    outcome: str | Iterable[str],
    adjustment: str | Iterable[str],
) -> FeatureProbability

The posterior probability that adjustment is a valid adjustment set.

Each distinct DAG is checked with the complete adjustment criterion (is_valid_adjustment_set()), and the masses of the DAGs where it holds are summed. Exact (mc_error 0) when every DAG was enumerated; otherwise an estimate from the draws (exact draws or MCMC) with its Monte Carlo error. A posterior probability under the method. On observational data these probabilities typically stay well below 1, because only the Markov equivalence class is identified; in the registered B4 run (simulated linear-Gaussian graphs, the true parents of the treatment under the BGe posterior) the medians were 0.37 at n = 200 and 0.48 at n = 1,000.

average_effect#

def average_effect(
    samples: GraphSamples,
    data: TabularData,
    treatment: str,
    outcome: str,
    *,
    rule: str = 'canonical',
    strategy: str | None = None,
    learner: Any = None,
    max_sets: int = 20,
    seed: int = 0,
    n_boot: int = 1000,
    **kwargs: Any,
) -> AveragedEffect

The effect of treatment on outcome, averaged over a posterior over DAGs.

For each distinct DAG: no directed path from the treatment to the outcome means the effect is 0; otherwise rule picks an adjustment set over the data's columns ("canonical", the set identify_effect() uses, or "minimal"). A variable of the posterior that data lacks is treated as latent, so a DAG can have no valid set; its mass is unidentified_mass. The DAGs are grouped by set, and the max_sets heaviest sets are each estimated once with estimate() (strategy and learner as there; the default is the family default, the doubly-robust AIPW, which needs a binary treatment; use strategy="dml" for a continuous one), seeded with seed.

The result mixes the group estimates by posterior mass, renormalized over the mass it covers: the zero-effect mass (at 0, standard error 0) plus the estimated sets. Its variance is the law of total variance, the mass-weighted mean of the squared standard errors (within) plus the mass-weighted variance of the estimates (between), reading each estimate as a normal around its value. Unidentified and truncated mass are left out and reported, never averaged in silently. Masses are posterior probabilities under the method, exact when every DAG was enumerated and otherwise from the draws, each with its Monte Carlo error.

Returns: An AveragedEffect.

bootstrap#

def bootstrap(
    data: TabularData | TimeSeriesData | PanelData,
    *,
    method: str = 'pc',
    n_bootstrap: int = 100,
    seed: int = 0,
    knowledge: DomainKnowledge | None = None,
    block_length: int | None = None,
    **kwargs: Any,
) -> BootstrapResult

Rerun a discovery method on bootstrap resamples and count what survives.

Answers the question a single run cannot: how stable is this edge under resampling? Each replicate draws a resample, runs discover(sample, method=method, knowledge=knowledge, **kwargs), and the result aggregates adjacency and glyph frequencies over the runs that finished. Extra keyword arguments flow to the algorithm unchanged (alpha, max_cond_size, score=..., max_lag=...), so any registered method, iid or temporal, sits behind the same call.

An iid method on a TabularData draws data.n_rows rows with replacement, the classic bootstrap; block_length does not apply there and raises. A temporal method ("pcmci") never resamples rows, which would shred the serial dependence it reads. On a TimeSeriesData it runs a moving-block bootstrap: ceil(T / L) blocks of L consecutive rows (block_length, default max(ceil(sqrt(T)), 2 * (max_lag + 1))) are drawn with replacement and enter the replicate as separate entities of a panel, "block_0", "block_1", ... with fresh time stamps. Blocks-as-entities is the point: gluing blocks back into one series would build lagged windows that straddle two unrelated fragments and let fake adjacencies in, while per-entity pooling never reads across a block seam. On a PanelData it runs the cluster bootstrap instead, drawing n_entities entities whole with replacement (relabeled in draw order so repeated draws stay distinct); block_length does not apply to a panel and raises. One honesty note: block resampling treats the series as stationary at the block scale, so a trend or a regime change inside the sampling window leaks into the frequencies. An iid method on temporal data still refuses, inside the replicate, exactly as discover() does.

The frequencies measure resampling variability, not Bayesian uncertainty; they are not posterior probabilities over graphs (ergodic has no graph posteriors). The run is deterministic: the same data, seed, and settings give the same result. A temporal method keys the frequencies on window-graph glyphs like "X[t-1] --> Y[t]".

A replicate that raises ValueError or RuntimeError (degenerate resamples happen: a constant column, a knowledge conflict on noisy data) counts in n_failed and is skipped; frequencies are over the successful runs only. When every replicate fails, the wrapper raises RuntimeError chaining the last error.

decode#

def decode(
    probabilities: EdgeProbabilities,
    method: str = 'greedy',
    *,
    threshold: float | None = None,
    knowledge: DomainKnowledge | None = None,
    samples: int = 1,
    seed: int | None = None,
    max_nodes: int = 20,
    audit: bool = True,
) -> DecodeResult

Decode edge probabilities into a DAG with a named, versioned decoder.

The input is read as pair states (to_pair_states()) and every decoder works on the edge weight w(i -> j) = log P(i -> j) - log P(no edge). The decoders:

• "greedy" (greedy_v1, the default): for each pair nominate its better direction when w > 0, or, with threshold=t, when the input value of that direction is above t; add nominations by decreasing weight and skip any that would close a cycle. On independent marginals, w > 0 is p > 0.5, AVICI's native rule; threshold lets a model's own cut (CDFM's adaptive one) choose the edges instead.
• "exact_map" (exact_map_v1): the DAG maximizing the factorized pair score, solved exactly by A* over topological orders.
• "expected_count" / "soft_f1" (expected_count_v2, soft_f1_v2): the best K edges compatible with the exact MAP order, K the rounded expected edge count or the size that maximizes a label-free soft-F1 proxy, both over the pair states with the hard knowledge applied (a forbidden direction holds no mass, a required one all of it). They may keep edges a MAP decoder drops. expected_count_v1 and soft_f1_v1 stay available with their 0.2 budget, which also counts the mass of directions the hard knowledge forbids; without knowledge the two versions agree.
• "pair_sampler" / "perturb_and_map" (pair_sampler_v1, perturb_and_map_v1): samples candidate DAGs, from independent pair-state draws with random-order cycle rejection, or from exact MAP under Gumbel-perturbed pair states. graph is their medoid (the candidate with the smallest total pair-state distance to the others, earliest on an exact tie; the research posterior_medoid can pick another tied candidate). Candidates are not posterior samples.

Ties are broken by variable order, so a decoder is deterministic given its inputs (and seed). The input is never changed and the global numpy random state is never touched.

Knowledge constrains the decode. Hard forbidden edges, tiers, roots, sinks and a skeleton restriction remove directed edges (their weight becomes minus infinity). Hard order (order, tiers) is enforced as one time order: the decoded graph has a topological order consistent with every order statement and every edge at once. That is exactly the tier reading, and stricter than violations() for separate order pairs (it only asks that no later variable be an ancestor of an earlier one), so a DAG that satisfies each pair alone but no single order is ruled out. Required edges go in first, whatever their weight, and must not form a cycle. The decoded graph is checked against the hard knowledge afterwards.

Soft statements (any builder's confidence= below 1) are weighed, not enforced. After the input is read as pair states, each pair's directed states are multiplied by exp(knowledge.log_odds(variables)) and the pair is renormalized; every decoder then runs unchanged on those states, so exact_map is exact on the combined score (pair score plus log_prior()). The _v2 cardinality decoders take their edge budget from the same weighted states with the hard constraints applied (a forbidden direction holds no mass, a required one all of it). A soft weight on a direction a hard constraint forbids, or on a pair a required edge settles, is dropped, so it changes nothing for any decoder. Greedy with threshold=t multiplies the input value's odds by the ratio on that direction before comparing it with t; for a statement that weighs both directions of a pair (forbid_adjacency, restrict_to) on pair_state input this is not quite the renormalized pair state the other decoders read, because the other direction's weight is left out.

confound(a, b) adds a <-> b after the directed part is decoded, so the graph and every candidate is an ADMG (see DecodeResult for a pair the decoder also joined with a directed edge).

With knowledge, diagnostics records an audit: per statement, the model's probability that it holds and whether it holds with and without the knowledge, the score the knowledge cost, and the strong conflicts (a pair the statement covers whose part the model gives less than 0.1), which also raise one KnowledgeConflictWarning. The audit decodes a second time without the knowledge (same settings and seed, so a sampler draws its samples candidates again). That can cost more than the decode itself: hard constraints often make an exact order search cheap, and the knowledge-free search has none. Pass audit=False to skip it.

Returns: A DecodeResult: graph a DAG over every variable (an ADMG with confounders), with the decoder, its parameters and the diagnostics recorded.

discover#

def discover(
    data: TabularData | TimeSeriesData | PanelData,
    method: str = 'pc',
    *,
    test: Any = None,
    score: Any = None,
    knowledge: DomainKnowledge | None = None,
    alpha: float = 0.05,
    **kwargs: Any,
) -> DiscoveryResult

Learn causal structure from data with the named method.

test feeds a constraint method (a CITest instance built without data, since discover binds it to data, or an alias like "linear"; a bound test runs through the algorithm class, as in PC(test=GCM(data)).fit()); score feeds a score-based method (a Score instance or an alias like "bic", "bdeu"); both default to "auto", picked from the column schema. Each algorithm declares which ingredient it consumes. Extra keyword arguments pass to the algorithm (max_cond_size, max_parents, node_limit, max_degree, tabu, max_lag).

The iid methods read a TabularData and reject temporal shapes (rows are not exchangeable); method="pcmci" reads a TimeSeriesData or PanelData and returns a window graph over lagged nodes like "X[t-1]".

A few methods consume neither ingredient because the machinery is the algorithm: "lingam" (regression and residual independence) and "glasso" (a sparse precision estimate). They reject test= and score=, and alpha does not apply (their own knobs, like lam or prune, pass through the keyword arguments).

method="welt:<model>" (for example "welt:arrow") runs a hosted Welt model through discover(), one paid Welt job. It needs pip install "ergodic[welt]", takes no test, score or alpha, and passes only client, seed, timeout and decoder through the keyword arguments; anything else raises before any request. decoder= asks for local decoding of the model's edge scores (decode()), which is what lets a non-DAG model such as "welt:avici" return a (decoded) graph; knowledge applies to that decoding and, without decoder=, implies greedy.

evaluate_graph#

def evaluate_graph(
    candidate: MixedGraph,
    source: EdgeProbabilities | GraphSamples,
    knowledge: DomainKnowledge | None = None,
    *,
    method: str = 'auto',
    max_exact: int = 15,
    samples: int = 10000,
    chains: int = 4,
    warmup: int = 500,
    thin: int = 4,
    seed: int = 0,
) -> GraphEvaluation

Read a candidate graph against a method's posterior over DAGs.

With edge probabilities this first builds the posterior with graph_posterior() (knowledge and the sampler settings go there); with GraphSamples, for example a score posterior, it reads that one, and the settings are not used (pass knowledge to graph_posterior(), not here).

The measures are those of DEC-082. For a DAG: the log-ratio to the best graph (the normalizer cancels), the credible level (the posterior mass of DAGs the method ranks strictly above the candidate, ties counted half; a Bayesian tail probability, not a p-value), every pair state as a claim, and the expected structural Hamming distance to a draw against that between two draws. For a CPDAG, PAG or ADMG: the claims that kind makes (see the module docstring for the reading of each mark) and, for a CPDAG, the equivalence class's probability. For independent marginals, P(acyclic) under the raw matrix: exact by the subset recursion up to max_exact variables, else from 20,000 seeded draws.

Pair-state claims are exact when the posterior is (enumeration or the subset recursion). Other claims, and the credible level, come from the draws with a Monte Carlo error unless every DAG was enumerated. All of them are posterior probabilities under the method: as calibrated as the method is.

Returns: A GraphEvaluation.

granger#

def granger(
    data: TimeSeriesData,
    *,
    max_lag: int = 1,
    alpha: float = 0.05,
) -> GrangerResult

Granger-test every ordered pair of variables in data.

For each effect Y, the full regression puts Y[t] on an intercept plus every variable's lags 1..max_lag; the restricted one drops the candidate cause's lags. All variables must be continuous. See the module docstring for what this does and does not mean.

graph_posterior#

def graph_posterior(
    source: EdgeProbabilities | TabularData,
    *,
    score: Score | str | None = None,
    knowledge: DomainKnowledge | None = None,
    method: str = 'auto',
    max_exact: int = 15,
    samples: int = 10000,
    chains: int = 4,
    warmup: int = 500,
    thin: int = 4,
    temperatures: tuple[float, ...] | None = None,
    seed: int = 0,
) -> GraphSamples

A posterior over DAGs under a method, exact when feasible, else by MCMC.

From edge probabilities the target is exp(pair log score(G) + log_prior(G)): the pair states the decoders read (independent marginals conditioned on no reciprocal pair, frequencies as they are, floored at 1e-12) and the knowledge's log_prior(). From data it is exp(score(G) + log_prior(G)), uniform over DAGs times the knowledge. BDeu and BGe are marginal likelihoods; a BIC-family score only approximates one and diagnostics["approximate"] says so.

Hard knowledge restricts the DAGs: banned edges, required edges, and the order statements read ancestrally (a before b: b is not an ancestor of a, as violations() reads them). Soft statements weigh the pair states. confound statements speak about no directed edge, so a DAG posterior ignores them and lists them in diagnostics["ignored_confounders"]; log_prior above is then that of the knowledge without them. A single variable has one DAG, so any method returns it enumerated.

method="auto" picks, in order: every DAG enumerated (up to 4 variables, 5 for a score); the exact subset recursion (a pair-factorized target, edge probabilities whose order knowledge reduces to banned edges, up to max_exact variables), which gives the exact normalizer and edge marginals and samples exact draws; otherwise chains structure MCMC chains. Chain 0 starts from the greedy decode (edge probabilities) or a hill climb (a score), the others from random DAGs; each runs warmup sweeps and keeps one draw every thin sweeps, samples // chains per chain (a sweep proposes one move per movable pair). Each chain is a parallel-tempering ladder, one replica per inverse temperature in temperatures, with draws kept from the untempered one; the tempered replicas carry the chain across valleys between modes that single-pair moves rarely cross. "exact" and "mcmc" force one; "exact" raises when it does not apply.

MCMC reports the rank-normalized split R-hat and bulk effective sample size of the log target and of every directed-edge indicator, and warns with PosteriorConvergenceWarning when the log target's R-hat is above 1.01, its effective sample size below 100 per chain, or an edge indicator with a posterior mean in [0.01, 0.99] has an R-hat above 1.05.

Every probability derived from the result is a posterior probability under the method: as good as the method's own calibration, which the docs report.

Returns: GraphSamples. diagnostics holds method, engine ("enumeration", "subset_recursion" or "structure_mcmc"), target (its kind, variables, whether it factorizes by pair, the ancestral constraints, approximate), seed, samples, log_normalizer (exact only), n_distinct, edge_estimator ("exact", "rao_blackwell" or "frequency"), best_source, ignored_confounders, approximate and note; for MCMC also chains, draws_per_chain, warmup, thin, temperatures, proposals_per_sweep, acceptance, swap_acceptance, rhat_log_target, ess_log_target, max_edge_rhat, min_edge_ess, converged and convergence_problems.

lagged_node#

def lagged_node(name: str, lag: int) -> str

The window-column name for name at lag steps back.

Lag zero reads "X[t]", lag k reads "X[t-k]". The name must carry no whitespace, since the graph layer parses edges on whitespace.

lagged_tabular#

def lagged_tabular(
    data: TimeSeriesData | PanelData,
    max_lag: int,
) -> TabularData

Stack the windows of data into a tabular frame of lagged columns.

Row t holds every variable at lags 0..max_lag; the column for variable v at lag k is named lagged_node(v, k) and holds v at time t - k. Column types carry over: continuous stays continuous, categorical keeps its categories, ordinal keeps its level order.

A time series needs at least max_lag + 8 rows. Panel windows are built within each entity (sorted entity order) and stacked; entities with fewer than max_lag + 2 rows are skipped with a warning, and it is an error when no rows survive. Unequal time steps warn too: lags count rows, not calendar time.

parse_lagged#

def parse_lagged(node: str) -> tuple[str, int]

Invert lagged_node(): "X[t-2]" becomes ("X", 2).

test_graph#

def test_graph(
    candidate: MixedGraph,
    data: TabularData,
    *,
    test: CITest | str | None = None,
    alpha: float = 0.05,
) -> GraphTest

Test the conditional independences a candidate graph implies against data.

Lists one implied independence per non-adjacent pair of observed variables (Shipley's basis set for a DAG; m-separation for an ADMG, MAG or PAG; see the module docstring), tests each, ranks them, adjusts their p-values by Holm's method, and combines them with Fisher's C, C = -2 * sum(log p_i), referred to a chi-square on 2k degrees of freedom. The chi-square reference assumes the k tests are independent, which Shipley argues for a DAG's basis set when the model holds (and only asymptotically for asymptotic tests such as Fisher z and G squared); for ADMGs, MAGs and PAGs it is an approximation. The combined p-value is the level-alpha decision (GraphTest.rejected). The Holm decisions locate which claims fail and hold the family-wise error rate at alpha whatever the dependence; rejecting when either one rejects spends more than alpha.

Only the default partial-correlation test is calibrated: its size was checked against alpha on linear-Gaussian data. Any other test carries a note on the result saying it is not calibrated; a GCM with its default nearest-neighbour learner and G squared with fewer than 1,000 rows are known to over-reject true graphs.

Claims are ranked by p-value, smallest first, with ties (such as p-values that underflow to 0.0) going to the larger |statistic|. Fisher's C is unaffected by the order.

What the p-value means: how surprising the data would be if every implied independence held. It cannot separate Markov-equivalent graphs (they imply the same independences), it never tests an edge that is present (only absent ones), failing to reject does not confirm the graph, and its power grows with the number of rows.

Returns: A GraphTest with Fisher's C, its p-value, and the ranked claims.