sirus 0.1.0

SIRUS (Stable and Interpretable RUle Set): interpretable rule-based models from random forests, scikit-learn compatible

sirus — SIRUS (Stable and Interpretable RUle Set) in Python

A scikit-learn compatible re-implementation of the R package sirus (CRAN mirror: cran/sirus) by Clément Bénard et al.

SIRUS is a regression and binary-classification algorithm based on random forests that takes the form of a short, stable list of rules. It keeps the aggregation principle of random forests but, instead of aggregating predictions, it selects the most frequent nodes of a forest of shallow, quantile-constrained trees. The result combines the simplicity of decision rules with accuracy close to random forests — and, crucially, the selected rules barely change when the data is perturbed.

Proportion of class 1 = 0.333 - Sample size n = 150
if Petal.Length < 1.7 then 1 (n=44) else 0.0566 (n=106)
if Petal.Width < 0.4 then 1 (n=41) else 0.0826 (n=109)
if Sepal.Length >= 5.6 then 0.033 (n=91) else 0.797 (n=59)
...

Validated head-to-head against the R package on identical data: the two implementations select the same rules (bit-identical predictions on the classification benchmarks) — see Benchmarks.

Installation

pip install sirus          # or:  uv add sirus

From source:

pip install git+https://github.com/asubbaswamy/sirus-py

Requires Python ≥ 3.10 with numpy, scipy, scikit-learn, joblib (installed automatically). pandas is optional (DataFrames, categorical variables, rule tables, sirus_cv), as is matplotlib (CV plots): pip install "sirus[pandas,plot]". The library is a single module — if you prefer, drop sirus.py next to your code instead of installing.

Quick start

from sirus import SirusClassifier, SirusRegressor

# ---- classification (binary) ------------------------------------------
clf = SirusClassifier(num_rules=10, random_state=42)
clf.fit(X_train, y_train)              # X: ndarray or DataFrame
clf.print_rules()                      # human-readable rule list
proba = clf.predict_proba(X_test)      # average of the rule outputs
pred  = clf.predict(X_test)            # threshold at 0.5

# ---- regression ---------------------------------------------------------
reg = SirusRegressor(num_rules=10, random_state=42)
reg.fit(X_train, y_train)              # rules combined by non-negative ridge
reg.print_rules()                      # weights, intercept, rules
y_hat = reg.predict(X_test)

# rules as data
df = clf.rules_to_dataframe()          # rule / frequency / outputs / weight
for rule in clf.rules_:                # Rule and Condition dataclasses
    print(rule.frequency, rule.conditions, rule.output_in, rule.output_out)

Both estimators are plain scikit-learn estimators: they support clone, get_params/set_params, cross_val_score, pipelines and grid search. See example.py for worked examples (iris, breast cancer, diabetes).

How it works (same six steps as the R package)

1. Discretization — numerical variables are binned on their empirical q-quantiles (default q=10), variables with ≤ discrete_limit distinct values keep their observed values as split points, and categorical variables (non-numeric DataFrame columns) are target-ordered as in ranger.
2. Forest — a random forest of depth-max_depth trees (default 2) is grown on the binned data, so every split falls on a quantile. By default the number of trees is chosen automatically: batches of num_trees_step trees are grown until the estimated stability of the rule selection reaches 1 - alpha (95%).
3. Path extraction — every tree node defines a path (a hyperrectangle); the frequency of each path across trees is computed.
4. Selection — the paths occurring in a fraction > p0 of the trees are selected, or simply the num_rules most frequent ones when p0=None.
5. Post-treatment — a rule is discarded when its indicator function is a linear combination of the intercept and higher-ranked selected rules (e.g. X1 >= 2 is dropped if X1 < 2 was already selected).
6. Aggregation — each rule outputs the mean training response inside or outside its hyperrectangle. For classification the rule outputs are simply averaged; for regression they are combined by a ridge regression constrained to non-negative coefficients, with the penalty chosen by cross-validation (the glmnet(alpha=0, lower.limits=0) step of the R package).

API mapping from R

R (sirus)	Python
sirus.fit(data, y, type='classif')	SirusClassifier(...).fit(X, y)
sirus.fit(data, y, type='reg')	SirusRegressor(...).fit(X, y)
sirus.predict(m, newdata)	predict_proba(X)[:, 1] / predict(X)
sirus.print(m)	m.print_rules()
num.rule, p0, num.rule.max	num_rules, p0, num_rules_max
q, discrete.limit, max.depth	q, discrete_limit, max_depth
num.trees, num.trees.step, alpha	num_trees, num_trees_step, alpha
mtry (default p/3)	mtry
replace, sample.fraction	replace, sample_fraction
m$rules, m$rules.out, m$proba	m.rules_ (conditions, outputs, frequency)
m$rule.weights, m$rule.glm	m.rule_weights_, m.intercept_, m.lambda_
m$num.trees, m$mean	m.num_trees_, m.mean_
sirus.cv(data, y, nfold, ncv)	sirus_cv(X, y, nfold=10, ncv=10, ...)
cv.grid$p0.pred, $p0.stab	cv.p0_pred, cv.p0_stab
cv.grid$error.grid.p0	cv.error_grid (pandas DataFrame)
sirus.plot.cv(cv.grid)	cv.plot() (matplotlib)

from sirus import sirus_cv, SirusRegressor

cv = sirus_cv(X, y, nfold=10, ncv=10, random_state=0)   # tune p0
model = SirusRegressor(p0=cv.p0_stab).fit(X, y)          # p0_pred for classif
cv.plot()                                                # error/stability path

sirus_cv follows the R implementation exactly: one model with a large rule list is fit per fold, and the whole 500-point p0 path is obtained by truncating its rule list (selection by any p0 is a prefix of the post-treatment output, because the redundancy filter is greedy-sequential). Error is the pooled out-of-fold 1-AUC (classification) or unexplained variance (regression); stability is the average proportion of rules shared by two models fit on distinct folds, compared in quantile-index space. For regression it is slow for the same reason the R version is (one non-negative-ridge cross-validation per rule-list prefix per fold); pass a fixed num_trees (a few thousand) to bound the runtime.

Implementation notes / differences from the R package

• The forest is built with scikit-learn regression trees on the binned data (the R package fits ranger regression trees on the numeric 0/1 output; variance splitting on a binary target is equivalent to Gini, so the split criterion matches). min_samples_split=5 reproduces ranger's default minimal node size.
• Tree paths are counted symbolically as sets of (feature, quantile, side) splits, exactly as defined in the papers and the R package; region-equivalent paths (e.g. a rule and a redundant refinement of it) are merged afterwards by the linear-dependence post-treatment, as in R.
• The post-treatment implements the paper's linear-dependence filter with a rank test on data resampled from the product of the marginal bin distributions — the exact procedure of paths.filter.d in the R sources (used there for depth > 2), augmented with the empirical sample. For depth ≤ 2 the R package uses an equivalent symbolic algorithm.
• The automatic stopping criterion estimates the expected Dice–Sørensen stability of the selection between two independent forests via a Gaussian approximation of the path frequencies, averaged over the range of thresholds that actually determines the fitted model. The trajectory and the resulting number of trees (≈10⁴ on small datasets) match the R behavior (live A/B on LA Ozone: R stops at 11000 trees, this port at 8000); the exact C++ metric may differ marginally.
• The non-negative ridge CV rebuilds glmnet's lambda path (lambda_max = max|Xᵀy|/(0.001·n), 100 log-spaced values, 10-fold CV) and solves each fit as a small exact NNLS problem, so no glmnet dependency is needed.
• Extras not in the R package: arbitrary binary labels (not just 0/1), n_jobs parallel forest growing, rules_to_dataframe().
• The R package's C++ forest is ~6–20× faster to fit than the sklearn-based forest here (~0.1 s vs 0.6–2 s for ~10⁴ trees on small data). Use n_jobs=-1 for repeated fits (e.g. sirus_cv); a one-off fit gains little because of the worker-pool start-up cost.

Benchmark against the R package

Two benchmarks live in benchmarks/ (results from July 2026, R sirus 0.3.3).

Live A/B on identical data (live_ab_*)

Both implementations are fit on identical train/test CSVs with identical settings (10 rules, 10000 trees, q=10, depth 2, same mtry). Summary of benchmarks/output/live_ab_report.md:

dataset (task)	shared rules (Dice)	max |Δfreq|	prediction agreement (test)
breast cancer (classif)	10/10 (1.00)	0.007	bit-identical (max Δ 5e-14), AUC 0.9891 both
categorical (classif)	10/10 (1.00)	0.010	bit-identical (max Δ 5e-14), AUC 0.7384 both
diabetes (regression)	9/10 (0.90)	0.011	r = 0.992, R² 0.331 (R) vs 0.337 (py)
LA ozone (regression)	9/10 (0.90)	0.009	r = 0.995, R² 0.667 (R) vs 0.683 (py)

Each regression disagreement is a frequency near-tie at rank 10 (ozone 0.048 vs 0.043, diabetes 0.051 vs 0.052) — which rule lands the last slot is seed noise. Classification predictions are bit-identical because the selected rules, their supports, and their outputs match exactly; regression predictions differ slightly through the cross-validated ridge penalty (different fold RNG). sirus.cv / sirus_cv agree too: on ozone p0_stab is 0.030 for both, and the error-vs-rules paths differ by ~0.01 (within CV noise at ncv=2).

Reproduce with R (install.packages("sirus")) available:

python benchmarks/live_ab_prepare.py     # shared CSVs + manifest
Rscript benchmarks/live_ab_run_r.R --cv  # R side  -> benchmarks/output/r_*.json
python benchmarks/live_ab_run_py.py --cv # Python  -> benchmarks/output/py_*.json
python benchmarks/live_ab_compare.py     # report  -> live_ab_report.md

Published rule list (benchmark_vs_r.py)

Comparison to Table 1 of the AISTATS 2021 paper (LA Ozone, produced by the R implementation with ~9000 trees): 9 of 11 published rules are recovered identically at num_rules=11 (11/11 at 13; the last two are frequency near-ties at the cutoff), rule outputs match exactly, frequencies within ±0.012, intercept −7.3 vs −7.8. Note the paper prints display-rounded thresholds (temp < 65 for the actual decile 65.4): the live R package run on the same CSV produces exactly the thresholds this port produces.

Development

git clone https://github.com/asubbaswamy/sirus-py && cd sirus-py
conda create -n sirus python=3.12 -y          # interpreter from conda ...
uv pip install --python "$(conda run -n sirus which python)" -e ".[dev]"
conda run -n sirus python -m pytest -q        # 17 tests

(Any interpreter works — plain uv venv + uv pip install -e ".[dev]" too.) CI runs the test suite on Linux/macOS/Windows, Python 3.10–3.13, plus the examples and the published-table benchmark.

Files

• sirus.py — the library (SirusClassifier, SirusRegressor, Rule, Condition, sirus_cv, SirusCVResult).
• example.py — worked examples on iris, breast cancer and diabetes.
• test_sirus.py — pytest suite.
• benchmarks/ — benchmark_vs_r.py (published-table comparison, no R needed) and the live_ab_* scripts (R-vs-Python A/B on identical data); LAozone.csv is the LA Ozone dataset from the ESL website.

References

• Bénard C., Biau G., Da Veiga S., Scornet E. (2021a). SIRUS: Stable and Interpretable RUle Set for classification. Electronic Journal of Statistics, 15:427–505.
• Bénard C., Biau G., Da Veiga S., Scornet E. (2021b). Interpretable Random Forests via Rule Extraction. AISTATS 2021, PMLR 130:937–945.
• Breiman L. (2001). Random forests. Machine Learning, 45:5–32.

License

MIT. This is an independent re-implementation of the published algorithm with the R sources read only as a reference — no code was taken from the GPL-3 R/C++ implementation. If you redistribute this package together with the R package, check license compatibility yourself.