GIES)#

Download this notebook

XGES, hill climbing, and Tabu search look for the graph with the best BIC score. XGES searches over equivalence classes, as GES does, with deletions first. Hill climbing adds, removes, or reverses one edge of a DAG at a time and stops when no move improves the score; Tabu search keeps going past that point and returns the best DAG it visits. All three return a CPDAG and suit linear-Gaussian data. Exact search finds the best-scoring DAG on a few variables, and GIES also learns from interventions.

Start#

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.9.0 from PyPI; local kernels use the installed environment.

import os
import sys

IN_COLAB = "google.colab" in sys.modules

if IN_COLAB:
    os.environ["ANDREY_DEVICE"] = "cpu"
    !pip install --quiet "andrey-core[viz]==0.9.0" pandas

import numpy as np
import pandas as pd
from IPython.display import display

import andrey
import andrey.metrics as metrics
import andrey.viz
from andrey.data import benchmarks

Data#

linear_gauss_er draws an Erdos-Renyi graph with linear mechanisms and Gaussian noise: here 10 variables and 2,000 samples.

seed = 7
dataset = benchmarks.scm("linear_gauss_er", d=10).sample(n=2000, seed=seed, scale="standardize")
X, truth = dataset.data, dataset.observed_graph
print(f"{X.shape[0]} observations, {X.shape[1]} variables, seed {seed}")
display(andrey.viz.draw(truth))

2000 observations, 10 variables, seed 7

Fit XGES, hill climbing, and Tabu search#

metadata["score"] is the BIC score of the returned graph, which each search minimizes. A larger lambda_value weighs the BIC penalty more and gives a sparser graph. With extended_search=False, XGES skips its extended search for a faster fit, and tabu_length sets how many recently visited DAGs Tabu search may not return to.

results = {
    "XGES": andrey.xges(X),
    "Hill climbing": andrey.hc(X),
    "Tabu search": andrey.tabu(X),
}
for name, result in results.items():
    print(f"{name}: {result.structure.kind}, score {result.metadata['score']:.1f}")
    display(andrey.viz.draw(result.structure))

XGES: cpdag, score -22406.2

Hill climbing: cpdag, score -22356.9

Tabu search: cpdag, score -22368.4

Compare with GES over ten datasets#

Structural Hamming distance (SHD) counts errors against the truth’s CPDAG; lower is better. The empty graph’s SHD is the number of true edges, so an estimate at or above it has not improved on returning no edges. One dataset can favor any method, so the table covers ten seeds and their median.

methods = {
    "GES": lambda X: andrey.ges(X),
    "XGES": lambda X: andrey.xges(X),
    "Hill climbing": lambda X: andrey.hc(X),
    "Tabu search": lambda X: andrey.tabu(X),
}
rows = {}
for s in range(10):
    draw = benchmarks.scm("linear_gauss_er", d=10).sample(n=2000, seed=s, scale="standardize")
    d = draw.observed_graph.n_nodes
    empty = andrey.GraphStructure.from_numpy(np.zeros((d, d), dtype=np.uint8), kind="cpdag")
    rows[s] = {"Empty graph": metrics.score(empty, draw.observed_graph)["shd"]}
    for name, fit in methods.items():
        rows[s][name] = metrics.score(fit(draw.data), draw.observed_graph)["shd"]
shd = pd.DataFrame(rows)
shd.columns.name = "seed"
shd["median"] = shd.median(axis=1)
shd

seed	0	1	2	3	4	5	6	7	8	9	median
Empty graph	15	15	15	15	15	15	15	15	15	15	15.0
GES	21	10	0	10	17	7	11	17	0	0	10.0
XGES	0	8	0	9	1	0	1	6	0	0	0.5
Hill climbing	28	6	0	8	25	16	22	11	4	13	12.0
Tabu search	6	0	0	1	20	9	16	12	0	0	3.5

Tabu search starts with hill climbing’s moves, so its score is never worse than hill climbing’s on the same data. The 0.3.0 release card compares XGES with the xges package and the other two with bnlearn, on larger graphs with ten rows per variable.

Exact search#

exact_search scores every parent set of every variable and returns a DAG with the best BIC score. Its run time at least doubles with each added variable, so keep it to about 20 variables.

exact = andrey.exact_search(X)
print("Exact search:", exact.structure.kind)
display(andrey.viz.draw(exact.structure))

Exact search: dag

Interventions: GIES#

gies also takes the rows’ interventions: True where a row’s value of a variable was set from outside. A chain from observational data stays undirected; here the last 1,000 rows set x1, which orients both of its edges. The columns come in an order drawn from the seed, and the interventions follow them.

CHAIN = [("x0", "x1"), ("x1", "x2")]  # the true edges, cause first


def sample_intervened(
    seed: int, n: int = 3000, shuffle: bool = True
) -> tuple[pd.DataFrame, np.ndarray]:
    """Draw n rows of the chain, the last 1,000 setting x1 from outside, and the interventions."""
    rng = np.random.default_rng(seed)
    x0 = rng.standard_normal(n)
    x1 = 0.8 * x0 + rng.standard_normal(n)
    x1[n - 1000 :] = rng.standard_normal(1000)  # the last 1,000 rows set x1 from outside
    x2 = 0.8 * x1 + rng.standard_normal(n)
    frame = pd.DataFrame({"x0": x0, "x1": x1, "x2": x2})
    frame = frame[rng.permutation(frame.columns)] if shuffle else frame
    intervened = np.zeros((n, 3), dtype=bool)
    intervened[n - 1000 :, frame.columns.get_loc("x1")] = True
    return frame, intervened


gies_out = andrey.gies(*sample_intervened(seed=1))
print(gies_out.structure.oriented_edges())
display(andrey.viz.draw(gies_out.structure))

[('x1', 'x2', 'directed'), ('x0', 'x1', 'directed')]

Use another dataset#

Change seed, or the d and n passed to benchmarks.scm and sample, and rerun.