pybnlearn

A Python port of bnlearn, Marco Scutari's R package for Bayesian network structure learning, parameter learning and inference.

Status: pre-release. 139 of bnlearn's 160 exported functions are ported and checked against R — see What works and Verified against R. Wheels build in CI for Linux (x86-64 and aarch64, manylinux and musllinux), macOS arm64 and Windows, and each one runs the whole parity suite before it is kept. Nobody outside this repository has used it yet, so expect the first releases to be alphas.

This is a port, not a drop-in replacement. Plotting (graphviz.*, strength.plot, the bn.fit.* lattice plots) and the conversions to igraph, gRain, graphNEL and lm are not here at all — they would be rewrites rather than ports. And direct_lingam reproduces R's causal ordering exactly but picks each node's parents with a score-based search where R uses glmnet's adaptive lasso: measured over eight cases it agrees on six and is less sparse on two, keeping arcs R prunes and never the reverse. Both are set out in full below.

What makes this a port rather than a reimplementation

The numerical core is bnlearn's C code. The sources under src/c/bnlearn/ are byte-identical to the CRAN 5.2.1 tarball; what makes them build without R is a compatibility layer (src/c/compat/) that reimplements the subset of the R C API bnlearn uses — the SEXP object model, attributes, allocVector, error() and about sixty other entry points. Forwarding headers let the vendored sources keep including <R.h> and friends by their usual names, so tracking a new bnlearn release means dropping in a new tarball.

One thing here is not a port: bnlearn has no exact inference over discrete networks of its own, so query() implements the junction tree directly rather than delegating to the gRain package as bnlearn does. Its tests compare numbers against gRain, which is a weaker claim than the rest of the suite makes, and src/pybnlearn/exact.py says why. Gaussian networks are a port: bnlearn turns those into a multivariate normal itself, and src/pybnlearn/mvnorm.py follows it.

R's own standalone maths library is vendored too (src/c/nmath/), along with R's Mersenne-Twister, so distribution functions and random streams agree with R bit for bit. That is what makes it meaningful to assert that results match rather than approximate: the parity suite compares against numbers generated by R itself, at a relative tolerance of 1e-11 — a bound set by measuring the worst disagreement across macOS and Linux (1.05e-12), not chosen for comfort. On one machine the agreement is 1e-12, which is what the fixtures were generated to. Where it cannot be that good, and why, is set out under what "matches R" does not cover.

Because bnlearn is GPL-2 | GPL-3, pybnlearn is distributed under the GPL v3 or later. See NOTICE for attribution.

Install

pip install bnlearn-port

import pybnlearn        # the install name and the import name differ

Installed as bnlearn-port, imported as pybnlearn. PyPI declines names too close to an existing project, and there is already a package called bnlearn there — a separate, MIT-licensed one built on pgmpy, unrelated to this and to Scutari's R package. The rule is doing its job: two Bayesian-network packages a prefix apart would be genuinely confusing. So the distribution took a name that says what it is, and the import name stayed put rather than churn every example. Python does this often enough to be unremarkable — scikit-learn imports as sklearn, Pillow as PIL.

Use

import pandas as pd
import pybnlearn

data = pd.read_csv("learning.test.csv", dtype="category")

# conditional independence testing
pybnlearn.ci_test(data, "A", "B", sx=["C"], test="mi")
# {'statistic': {'mi': 2347.258587894275}, 'parameter': {'df': 12.0}, ...}

# structure learning
net = pybnlearn.hc(data, score="bic")
net.modelstring()          # '[A][C][F][B|A][D|A:C][E|B:F]'
net.arcs                   # [('A', 'B'), ('A', 'D'), ...]
pybnlearn.score(net, data) # -24006.734232498

# constraint-based learning: these return a partially directed graph
pybnlearn.gs(data, test="mi", alpha=0.05).arcs
pybnlearn.iamb(data).arcs
pybnlearn.inter_iamb(data, undirected=True).arcs

# parameter learning
fitted = pybnlearn.fit(net, data)
fitted["B"].probabilities        # the conditional probability table
fitted["B"].as_frame()           # ... in long format

# inference: seeded sampling reproduces R's set.seed() exactly
pybnlearn.set_seed(1)
pybnlearn.cpquery(fitted, {"B": "a"}, {"A": "a"}, method="ls", n=10000)
pybnlearn.rbn(fitted, 100)       # sample from the network

# ... or compute it exactly, no sampling involved
pybnlearn.query(fitted, "B", {"A": "a"}).values

# resampling
pybnlearn.set_seed(1)
pybnlearn.boot_strength(data, algorithm="hc", replicates=200)
pybnlearn.bn_cv(data, "hc", k=10).mean

# constraints
pybnlearn.hc(data, whitelist=[("A", "F")], blacklist=[("A", "B")], maxp=3)

# graph priors: what you believe about the structure before seeing the data.
# A whitelist forces an arc; a prior only leans on it, so on data this
# decisive all three still land on the network above -- the scores differ.
pybnlearn.score(net, data, type="bde", prior="vsp", beta=0.1)
pybnlearn.score(net, data, type="bde", prior="marginal", beta=0.9)
pybnlearn.score(net, data, type="bde", prior="cs", beta=pd.DataFrame(
    {"from": ["A", "B"], "to": ["B", "C"], "prob": [0.9, 0.05]}))
pybnlearn.hc(data, score="bde", prior="vsp", beta=0.1)

# the whole graph at once, in and out
pybnlearn.amat(net)              # a labelled square frame of 0s and 1s
pybnlearn.set_amat(net, matrix)  # ... and back to a network

Two notes on pandas.read_csv and categorical data

Missing values. read_csv treats NA, None, N/A and similar strings as missing by default. bnlearn's own insurance data set has a category literally called "None", so the default reads 14681 of its rows as missing. Pass keep_default_na=False, na_values=[] when the labels are meant to be data. pybnlearn raises rather than letting missing values reach the C core, but it cannot tell that a category was lost before it ever saw the frame.

Level order. read_csv orders categories alphabetically. Most of the library does not care — mutual information and the network scores are invariant to it — but a conditional probability table is indexed by level, so fit() returns its axes in whatever order the categories are in. If the order matters to you, set it explicitly:

data["Species"] = data["Species"].cat.reorder_categories(["Sagrei", "Distichus"])

What works

139 of bnlearn's 160 exported functions are covered, and nothing of substance is left. Fifteen of the remaining twenty-one are deliberately out of scope -- the nine plotting functions and the six conversions to igraph, graphNEL, gRain and lm -- because they would be rewrites rather than ports, with nothing to check against. The other six are R's assignment forms (nodes<-, arcs<-, parents<-, children<-, modelstring<-, alst<-), whose effect is available here as ordinary functions.

Three things are deliberately not reproduced, and each is documented where it lives. direct_lingam gives R's causal ordering exactly, but picks the arcs with a score-based search rather than with glmnet's adaptive lasso, which is not vendored. tests/parity/test_lingam_arcs.py measures how far apart that leaves them -- six of eight cases identical, two differing only by arcs this keeps and R prunes -- and fails if the gap widens, narrows or changes direction. Porting glmnet would not close it cleanly: its lambda path stops when the gain in deviance falls below 1e-5 times the deviance, and that comparison came down to 7.865e-06 against 7.634e-06 on the first data set tried, so the number of lambdas -- and thus which one is selected -- flips with the convergence tolerance. Matching it needs glmnet's arithmetic rather than its method. mi = "gkernel" raises rather than returning a number: bnlearn builds its Gram matrix with a matrix product where it means an elementwise one, leaving a condition number around 1e20, so R's own answer is floating-point noise and cannot be reproduced by anything. And count_graphs("dags-with-r-arcs", r = 0) returns one rather than R's n: R's inner loop runs seq(from = 1, to = 0), which counts downwards, and the extra pass overwrites the count -- unlike the other two this one is checkable, and the counts for a given number of nodes only add up to the total with the corrected value. Exact inference covers discrete and Gaussian networks, but not mixtures of the two.

Three limits on how far "matches R" goes, all found by running the code somewhere new rather than reasoned about in advance.

The last digit belongs to the platform, not to this package. The vendored nmath and Mersenne-Twister make distributions and random streams agree with R exactly, but nmath itself calls log, exp and lgamma from the system maths library, and glibc and Apple's libm differ in the last unit in the last place. Nothing vendored can fix that. Measured across macOS and Linux the worst relative disagreement over the whole suite is 1.05e-12, so the parity tests compare at a relative tolerance of 1e-11 — ten-fold headroom over what was observed, tight enough that a real regression still fails. On a single machine agreement is far better, at the 1e-12 the fixtures were generated to.

BLAS is the suspect people reach for here, and it is the wrong one: building against scipy-openblas32 instead of Accelerate — a genuinely different BLAS on the same machine — changes nothing at all, and the whole suite passes unaltered.

Ill-conditioned data determine fewer digits than they appear to. For a data set built with two columns correlated to 1 − 5e-13, a partial correlation agrees with R only to about four digits — and measured against exact rational arithmetic, neither answer is better: R is wrong by 6.6e-5, this by 2.9e-4, against a forward error bound of 1.5e-3. The data do not determine those digits. test_the_ill_conditioned_disagreement_is_not_a_defect does that arithmetic rather than asserting it.

Some results are decided by ties, and those cannot be compared at all. Hartemink's discretization merges 30 initial bins over 88 observations on the marks data, so it chooses between differences that small — and a one-ulp change to a single observation, which is the same data by any reasonable reading, moves the answer in 4 of 12 trials. There the tests compare the levels and their number, which are reproducible, and not which side of a boundary each observation fell. test_hartemink_on_marks_really_is_decided_by_ties measures the sensitivity, so the exemption cannot quietly widen into covering a real disagreement.

One inconsistency is reproduced rather than smoothed over, because smoothing it over would be the divergence: bdj accepts prior and beta and then ignores them, scoring as if the prior were uniform -- while still being reported as non-decomposable under cs and marginal, which moves the search even though the score stands still. pybnlearn does the same, and the fixtures pin down both halves of it.

Verified against R

pytest runs 9381 checks, 8707 of which compare directly against values produced by R 4.6.1 with bnlearn 5.2.1:

• 318 conditional independence tests across discrete and Gaussian data, each comparing the statistic, the degrees of freedom and the p-value;
• 82 hill-climbing runs across 8 data sets, 13 scores, non-default hyperparameters, whitelists, blacklists and parent limits, each comparing the arc set, the model string and the per-node scores;
• 231 constraint-based runs across pc_stable, gs, iamb, inter_iamb, iamb_fdr, mmpc and si_hiton_pc, 6 data sets, 7 independence tests, several significance levels, constraint sets and both directed and undirected output, comparing the arc set including direction;
• 64 hybrid runs across mmhc and rsmax2, covering every ported restrict/maximize pair and arguments passed through to each phase;
• 110 hpc and h2pc runs across 7 data sets -- including 37-node alarm, where the two superset stages earn their keep and an off-by-one in them would be least visible -- 6 independence tests, conditioning-set limits and constraints;
• 63 tabu searches across 8 data sets, 9 scores, tabu list sizes from 1 to 30, constraints and parent limits — 13 of which R's tabu resolves differently from R's hc, so the tabu-specific paths are actually covered rather than incidentally agreeing with hill climbing;
• 55 fitted networks, comparing every cell of every conditional probability table and every regression coefficient, over three structures per data set and both the maximum likelihood and Bayesian estimators;
• 61 seeded simulation and inference results -- generated data compared observation by observation, and conditional probabilities compared exactly rather than statistically, since both implementations draw the same numbers in the same order from R's Mersenne-Twister;
• 61 bootstrap and cross-validation results, again exact rather than statistical, including R's sample() itself -- everything here rests on drawing the same rows in the same order;
• 48 predictions and prediction-based losses, covering both prediction methods and all five predictive cross-validation losses;
• 42 classifier structures and class posteriors, over four data sets and several class variables, including the tree root that decides how a TAN's feature tree is oriented;
• 40 exact inference results for discrete networks checked against gRain -- marginals, conditionals on up to three variables, joints, and exact prediction;
• 82 Gaussian exact inference results, which unlike the discrete ones are genuine parity fixtures: global means and covariances, the factorisation back into a network, conditional expectations, and exact prediction of every node of every structure -- including a network with a deterministic node, whose covariance matrix is singular and which therefore exercises the diagonal patching R does rather than failing;
• 70 conditional Gaussian results over two mixed data sets: structure learning with all four -cg scores under both hc and tabu, and parameter learning across all three of the estimators mle-cg dispatches to;
• 138 arc-strength and network-averaging results: p-values and score differences for six data sets and thirteen criteria, bootstrap strengths with the inclusion threshold they carry, and the averaged network at four thresholds -- plus thirteen strength vectors of awkward shapes checking that R's optimiser is reproduced rather than approximated, since the threshold it finds is used as a cutoff;
• 632 checks of the node and arc utilities, enumerated rather than chosen: nine graphs -- including one whose undirected part is not chordal and one with a real directed cycle -- crossed with every ordered pair of nodes and all five arc operations;
• 148 interchange-format results: every node of four networks written by R as BIF, DSC and NET and read back here, which is what pins down where each conditional distribution belongs -- NET does not record it at all;
• 24 hand-built networks from custom_fit, compared as parameters and then sampled from and queried, since a network with an axis in the wrong place reads back correctly and samples from a different distribution;
• 524 checks of d-separation, colliders, consistent extensions, structural intervention distance, random graph generation and discretization -- d-separation enumerated over every pair of nodes in six graphs and three kinds of conditioning set, and the generated graphs compared arc for arc rather than statistically, since they come from R's own generator;
• 245 checks of fast.iamb, local structure learning, entropy and Kullback-Leibler divergence, the last two comparing pybnlearn's junction tree against gRain's;
• 40 incomplete-data results: parameters estimated from data with gaps, all three imputation methods, and structural EM under both maximisers -- plus a latent variable, never observed at all, which is the case EM cannot start from an empty network on;
• 183 causal results: interventions on every node of four graph shapes, twin networks, and counterfactuals with and without node merging -- plus the parameterised half, where a node's residual variance becomes a node of its own feeding both copies;
• 18 LiNGAM causal orderings, including a data set with genuinely non-Gaussian noise, which is the assumption the method rests on;
• 139 results for the remaining utilities: seeded perturbations and the hill-climbing restarts built on them, every extension of an equivalence class, exact graph counts, Bayes-factor arc strengths in extended precision, deduplication, and bootstrapping an arbitrary statistic;
• 4727 results from an exhaustive sweep, which exists because the fixtures above are cases chosen by hand and share the blind spots of the hand that chose them. Here the arguments are crossed rather than varied one at a time; the ten data sets are awkward on purpose — two observations per contingency-table cell, eleven levels against two, 98% of the mass on one level, near-determinism, more variables than the sample supports, nearly collinear columns, ten orders of magnitude between columns, Cauchy tails, thirty rows; and 60 combinations have no solution, so the paths that refuse are exercised too. Every score on every data set with its per-node breakdown, every independence test over every ordered pair at three conditioning-set sizes, every constraint algorithm against every test, and ntests() for all eleven searches — the last being the sharpest check there is, since two different search paths can reach the same arcs but not the same amount of work. Writing it found four defects the hand-chosen fixtures could not: nparams refused a fitted network and miscounted Gaussian nodes, connected_components returned a bare string for a one-node component, and hc/tabu accepted whitelists R rejects;
• 22 adjacency-matrix results over ten graphs, including two partially directed ones -- an undirected arc is a symmetric pair of ones, so the matrix cannot tell it from two opposed directed arcs, and the round trip is recorded rather than assumed. The row/column convention and the arc order coming back out are compared too: both are invisible in a same-shaped answer and decide how everything downstream reads;
• 235 results for the non-uniform graph priors: four Castelo & Siebes completions, and nine prior settings crossed with three data sets, five scores and both hc and tabu -- which is what pins down the non-decomposable path, since a prior that touches every node makes a move's effect wider than its cached delta. hc and tabu are expected to disagree with each other under the marginal prior, and do, on the same pair of networks R's do;
• 27 checks of the graph utilities: CPDAG, moral graph and skeleton for six learned networks, shd/hamming/compare over five network pairs, model2network round trips, and chow_liu and aracne on six data sets;
• set.seed(42) reproduces R's uniform and normal streams to 15 digits.

Regenerate the fixtures (needs R with bnlearn installed):

Rscript tools/gen_r_fixtures.R
Rscript tools/gen_r_hc_fixtures.R
Rscript tools/gen_r_constraint_fixtures.R
Rscript tools/gen_r_graph_fixtures.R
Rscript tools/gen_r_tabu_fixtures.R
Rscript tools/gen_r_hybrid_fixtures.R
Rscript tools/gen_r_fit_fixtures.R
Rscript tools/gen_r_levels.R
Rscript tools/gen_r_inference_fixtures.R
Rscript tools/gen_r_resampling_fixtures.R
Rscript tools/gen_r_predict_fixtures.R
Rscript tools/gen_r_classifier_fixtures.R
Rscript tools/gen_r_exact_fixtures.R      # also needs gRain
Rscript tools/gen_r_cg_fixtures.R
Rscript tools/gen_r_mvnorm_fixtures.R
Rscript tools/gen_r_strength_fixtures.R
Rscript tools/gen_r_nodes_fixtures.R      # also needs igraph
Rscript tools/gen_r_custom_fixtures.R
Rscript tools/gen_r_foreign_fixtures.R
Rscript tools/gen_r_analysis_fixtures.R
Rscript tools/gen_r_local_fixtures.R     # also needs gRain
Rscript tools/gen_r_hpc_fixtures.R
Rscript tools/gen_r_missing_fixtures.R   # also needs gRain
Rscript tools/gen_r_causal_fixtures.R
Rscript tools/gen_r_lingam_fixtures.R
Rscript tools/gen_r_misc_fixtures.R      # also needs gmp and Rmpfr
Rscript tools/gen_r_priors_fixtures.R
Rscript tools/gen_r_amat_fixtures.R
Rscript tools/gen_r_sweep_fixtures.R
Rscript tools/gen_r_rejection_fixtures.R
Rscript tools/gen_r_lingam_arcs_fixtures.R  # also needs glmnet

Performance

Roughly two to three times slower than R on hc — 0.64s vs 0.32s for the 37-node alarm data set — because the network is marshalled into R objects afresh on each iteration of the search. Memory is flat across repeated runs.

Threads and processes

One call runs at a time per process. bnlearn's C keeps its state in process-wide statics — the arena results are allocated in, the interned symbol table, the random number generator, the jmp_buf error() unwinds to — so every entry point takes a lock. Calling from several threads is safe; it does not go faster.

This is bnlearn's own model rather than a limitation added here. R cannot be called from two threads at all, the vendored C contains no OpenMP and no pthreads, and bnlearn's cluster argument takes a parallel::makeCluster() cluster — which is separate R processes, each with its own memory.

So for real parallelism, use processes, exactly as R does:

from concurrent.futures import ProcessPoolExecutor

with ProcessPoolExecutor() as pool:                  # ~ parallel::makeCluster
    networks = list(pool.map(pybnlearn.hc, folds))

Each process gets its own arena and takes no lock. Seed inside the worker: the generator is per-process, so set_seed in the parent does not reach it.

Building from source

Needs a C compiler and a BLAS/LAPACK. No Fortran: R's dqrdc2 and dqrsl are translated to C in src/c/linpack/, and tools/check_linpack.sh checks the translation against the Fortran original bit for bit.

pip install -e . --no-build-isolation

To build and exercise the C core without Python in the loop — useful when bisecting a numerical disagreement with R:

tools/build_core.sh

Licence

GPL-3.0-or-later. See LICENSE and NOTICE.