{ "cells": [ { "cell_type": "markdown", "id": "permutation_based-0", "metadata": {}, "source": [ "# Permutation-based search (BOSS / GRaSP)\n", "\n", "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Abel-ai-lab/andrey/blob/main/examples/permutation_based.ipynb)\n", "\n", "BOSS and GRaSP search over orderings of the variables. For an ordering, a BIC score chooses each\n", "variable's parents from the variables before it; the search keeps the ordering with the best score.\n", "Both return a CPDAG, like GES, and suit linear-Gaussian data." ] }, { "cell_type": "markdown", "id": "7fb27b941602401d91542211134fc71a", "metadata": {}, "source": [ "## Start\n", "\n", "Choose a CPU runtime and run this cell first. In Colab it installs `andrey-core==0.9.0` from PyPI;\n", "local kernels use the installed environment." ] }, { "cell_type": "code", "execution_count": 1, "id": "acae54e37e7d407bbb7b55eff062a284", "metadata": { "tags": [ "colab-setup" ] }, "outputs": [], "source": [ "import os\n", "import sys\n", "\n", "IN_COLAB = \"google.colab\" in sys.modules\n", "\n", "if IN_COLAB:\n", " os.environ[\"ANDREY_DEVICE\"] = \"cpu\"\n", " !pip install --quiet \"andrey-core[viz]==0.9.0\" pandas" ] }, { "cell_type": "code", "execution_count": 2, "id": "permutation_based-1", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from IPython.display import display\n", "\n", "import andrey\n", "import andrey.metrics as metrics\n", "import andrey.viz\n", "from andrey.data import benchmarks" ] }, { "cell_type": "markdown", "id": "permutation_based-2", "metadata": {}, "source": [ "## Data\n", "\n", "`linear_gauss_er` draws an Erdos-Renyi graph with linear mechanisms and Gaussian noise: here 10\n", "variables and 2,000 samples." ] }, { "cell_type": "code", "execution_count": 3, "id": "permutation_based-3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2000 observations, 10 variables, seed 7\n" ] }, { "data": { "text/html": [ "0123456789" ], "text/plain": [ "'0123456789'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "seed = 7\n", "dataset = benchmarks.scm(\"linear_gauss_er\", d=10).sample(n=2000, seed=seed, scale=\"standardize\")\n", "X, truth = dataset.data, dataset.observed_graph\n", "print(f\"{X.shape[0]} observations, {X.shape[1]} variables, seed {seed}\")\n", "display(andrey.viz.draw(truth))" ] }, { "cell_type": "markdown", "id": "permutation_based-4", "metadata": {}, "source": [ "## Fit BOSS and GRaSP\n", "\n", "BOSS starts from an ordering computed from the data and visits the variables in a random order;\n", "GRaSP starts from a random ordering. `seed` fixes the random choices, so a run repeats exactly;\n", "without `seed`, the `ANDREY_SEED` environment variable sets it (default 0). `metadata[\"score\"]` is\n", "the BIC objective of the final ordering, which the search minimizes. A larger `lambda_value` weighs\n", "the BIC penalty more and gives a sparser graph." ] }, { "cell_type": "code", "execution_count": 4, "id": "permutation_based-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BOSS: cpdag, score -22398.8\n" ] }, { "data": { "text/html": [ "0123456789" ], "text/plain": [ "'0123456789'" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "GRaSP: cpdag, score -22406.2\n" ] }, { "data": { "text/html": [ "0123456789" ], "text/plain": [ "'0123456789'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = {\n", " \"BOSS\": andrey.boss(X, seed=0),\n", " \"GRaSP\": andrey.grasp(X, seed=0),\n", "}\n", "for name, result in results.items():\n", " print(f\"{name}: {result.structure.kind}, score {result.metadata['score']:.1f}\")\n", " display(andrey.viz.draw(result.structure))" ] }, { "cell_type": "markdown", "id": "permutation_based-6", "metadata": {}, "source": [ "## Compare with GES over ten datasets\n", "\n", "GES searches over CPDAGs with a BIC score as well. Structural Hamming distance (SHD) counts errors\n", "against the truth's CPDAG; lower is better. The empty graph's SHD is the number of true edges, so\n", "an estimate at or above it has not improved on returning no edges. One dataset can favor any\n", "method, so the table covers ten seeds and their median." ] }, { "cell_type": "code", "execution_count": 5, "id": "permutation_based-7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
seed0123456789median
Empty graph1515151515151515151515.0
GES211001017711170010.0
BOSS6801011013110149.0
GRaSP60011016000.5
\n", "
" ], "text/plain": [ "seed 0 1 2 3 4 5 6 7 8 9 median\n", "Empty graph 15 15 15 15 15 15 15 15 15 15 15.0\n", "GES 21 10 0 10 17 7 11 17 0 0 10.0\n", "BOSS 6 8 0 10 1 10 13 11 0 14 9.0\n", "GRaSP 6 0 0 1 1 0 1 6 0 0 0.5" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "methods = {\n", " \"GES\": lambda X: andrey.ges(X),\n", " \"BOSS\": lambda X: andrey.boss(X, seed=0),\n", " \"GRaSP\": lambda X: andrey.grasp(X, seed=0),\n", "}\n", "rows = {}\n", "for s in range(10):\n", " draw = benchmarks.scm(\"linear_gauss_er\", d=10).sample(n=2000, seed=s, scale=\"standardize\")\n", " d = draw.observed_graph.n_nodes\n", " empty = andrey.GraphStructure.from_numpy(np.zeros((d, d), dtype=np.uint8), kind=\"cpdag\")\n", " rows[s] = {\"Empty graph\": metrics.score(empty, draw.observed_graph)[\"shd\"]}\n", " for name, fit in methods.items():\n", " rows[s][name] = metrics.score(fit(draw.data), draw.observed_graph)[\"shd\"]\n", "shd = pd.DataFrame(rows)\n", "shd.columns.name = \"seed\"\n", "shd[\"median\"] = shd.median(axis=1)\n", "shd" ] }, { "cell_type": "markdown", "id": "permutation_based-8", "metadata": {}, "source": [ "A seed where GES reaches the empty graph's SHD is one where it has not improved on returning\n", "no edges. The [release cards](https://andrey.abel.ai/docs/benchmarks.html) show the same on\n", "larger graphs: GES stays close to the empty graph's SHD, while BOSS and GRaSP come close to zero.\n", "\n", "## Use another dataset\n", "\n", "Change `seed`, or the `d` and `n` passed to `benchmarks.scm` and `sample`, and rerun. The\n", "benchmarks page also compares both methods with `causal-learn`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }