Running in Colab? Install first: %pip install beamfeat (before the PyPI release: %pip install git+https://github.com/LD-Shell/beamFeat).

Reading a fit¶

The three notebooks before this one plant an answer and check that beamfeat finds it. Real fits do not arrive that way. They arrive as an equation you did not write, containing columns you did not expect, with fdr_controlled_ reporting True and no obvious way to tell a triumph from a coincidence.

This notebook is about that moment. We take a problem with a known physical answer, hide it among noise, and then get it wrong on purpose — because the ways it goes wrong are more instructive than the way it goes right, and every one of them has a fix.

In [1]:

import re
import warnings

import numpy as np
import pandas as pd

# beamfeat's warnings, formatted without this machine's file paths
warnings.filterwarnings("default")
warnings.formatwarning = lambda message, category, *_: f"{category.__name__}: {message}\n"

from beamfeat import BeamFeatRegressor

rng = np.random.default_rng(99)
n = 300

import re import warnings import numpy as np import pandas as pd # beamfeat's warnings, formatted without this machine's file paths warnings.filterwarnings("default") warnings.formatwarning = lambda message, category, *_: f"{category.__name__}: {message}\n" from beamfeat import BeamFeatRegressor rng = np.random.default_rng(99) n = 300

The problem¶

Hydraulic power of a centrifugal pump, P = rho * g * Q * H / eta: density, gravity, flow rate, head, efficiency. Five real variables, buried among eight columns of plausible-looking nonsense — the ambient noise any plant historian is full of.

Note that g is a constant. Keep an eye on it.

In [2]:

columns = {
    "rho": rng.uniform(800, 1200, n),   # density, kg/m^3
    "g": np.full(n, 9.81),              # gravity, m/s^2 — constant
    "Q": rng.uniform(0.01, 0.5, n),     # flow rate, m^3/s
    "H": rng.uniform(10, 100, n),       # head, m
    "eta": rng.uniform(0.6, 0.9, n),    # efficiency, dimensionless
}
power = columns["rho"] * columns["g"] * columns["Q"] * columns["H"] / columns["eta"]
y = power + rng.normal(0, 0.02 * power.std(), n)

for index in range(8):
    columns[f"noise_{index}"] = rng.uniform(0, 50, n)

X = pd.DataFrame(columns)
print(f"{X.shape[0]} rows, {X.shape[1]} columns, {X.shape[1] - 5} of them irrelevant")

columns = { "rho": rng.uniform(800, 1200, n), # density, kg/m^3 "g": np.full(n, 9.81), # gravity, m/s^2 — constant "Q": rng.uniform(0.01, 0.5, n), # flow rate, m^3/s "H": rng.uniform(10, 100, n), # head, m "eta": rng.uniform(0.6, 0.9, n), # efficiency, dimensionless } power = columns["rho"] * columns["g"] * columns["Q"] * columns["H"] / columns["eta"] y = power + rng.normal(0, 0.02 * power.std(), n) for index in range(8): columns[f"noise_{index}"] = rng.uniform(0, 50, n) X = pd.DataFrame(columns) print(f"{X.shape[0]} rows, {X.shape[1]} columns, {X.shape[1] - 5} of them irrelevant")

300 rows, 13 columns, 8 of them irrelevant

Start at the defaults¶

Always. The defaults are max_depth=2, beam_width=50.

In [3]:

model = BeamFeatRegressor(random_state=99).fit(X, y)

print(model.equation())
print(f"R^2: {model.score(X, y):.4f}")

model = BeamFeatRegressor(random_state=99).fit(X, y) print(model.equation()) print(f"R^2: {model.score(X, y):.4f}")

UserWarning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are meant to carry data.

y = 9.799*((H / eta) * (Q * rho)) + 279.9   [1 of 50 certified terms; parsimony=None prints all 50]
R^2: 0.9996

(H / eta) * (Q * rho) is rho * Q * H / eta — the physics, exactly, with none of the eight noise columns in it.

The suffix is the equation being honest about its own scope: screening certified fifty formulas, most of them pairing the physics core with a noise column, and the parsimony step kept the one that carries the model. Those fifty are the set the guarantee covers. parsimony=None prints them all, and the section below on the marginal null explains why the noisy ones are there and why certifying them is not an error.

Where did gravity go?¶

It is in the coefficient. g never varies, so no expression containing it can be distinguished from the same expression without it: rho * g is just 9.81 * rho, perfectly collinear with rho, and the redundancy pass keeps one representative of the pair. The scaling then lands in the linear model, which is why the fitted coefficient reads 9.80 rather than 1. Under parsimony=None the same scaling is shared across fifty collinear terms and the leading coefficient reads lower, which is why a constant is easiest to read off the compact fit.

This generalises. A multiplicative constant is absorbed into the coefficient rather than recovered as a symbol, so its absence from the formula is not a miss. Check the coefficient before concluding a factor was lost.

Searching harder makes it worse¶

Five variables feels like it should need a deep search. Let us give it one: depth 3, and a wider beam. verbose=1 shows what each stage did.

In [4]:

deep = BeamFeatRegressor(
    random_state=99, max_depth=3, beam_width=60, verbose=1
).fit(X, y)

print()
print(deep.equation())
print(f"R^2: {deep.score(X, y):.4f}")
print(f"fdr_controlled_: {deep.fdr_controlled_}")

deep = BeamFeatRegressor( random_state=99, max_depth=3, beam_width=60, verbose=1 ).fit(X, y) print() print(deep.equation()) print(f"R^2: {deep.score(X, y):.4f}") print(f"fdr_controlled_: {deep.fdr_controlled_}")

[beamfeat] fit: 300 rows x 13 columns, split 150 search / 150 selection

UserWarning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are meant to carry data.

[beamfeat] search: depth 3, beam 60 -> 42413 proposed, 60 candidates (4.33s)
[beamfeat] screen: 60 candidates, permutation (BY) at q=0.1 -> 60 certified
[beamfeat] parsimony: 60 certified -> 1 term (|S|/|S'| = 60.00)
[beamfeat] result set: 1 term in the fitted equation (the guarantee covers the certified set, not this subset)
[beamfeat] result: y = 9.743*(((H / eta) * (Q * rho)) - ((H + noise_0) / (Q - eta))) - 483.6   [1 of 60 certified terms; parsimony=None prints all 60]  [FDR controlled over the screened set]

y = 9.743*(((H / eta) * (Q * rho)) - ((H + noise_0) / (Q - eta))) - 483.6   [1 of 60 certified terms; parsimony=None prints all 60]
R^2: 0.9996
fdr_controlled_: True

The physics is still in there — (H / eta) * (Q * rho) — but a junk term containing noise_0 has been bolted onto it, and fdr_controlled_ is still True. The R^2 is unchanged to four decimals. We paid three times the runtime for an equation nobody can read.

Why the flag is still True¶

This is the single most important thing to understand about the guarantee, and it is not a defect.

The null being tested is marginal: is this expression, as a whole, independent of the target? An expression that combines the true core with a noise column is emphatically not independent of the target — the core sees to that — so its null is false and selecting it is not an error under the criterion. The guarantee certifies that a selected formula is not pure noise. It does not certify that every column inside it earns its place.

So fdr_controlled_ = True means what it says, and it is not the property you wanted here. Noise hitchhikes on a strong signal.

Depth is not the number of variables¶

The instinct that sent us to max_depth=3 was wrong arithmetic. Depth counts composition levels, not variables, and a binary operator combines any two surviving expressions — not a variable onto a chain.

So rho * Q * H / eta, four variables, is depth two:

depth	built
1	(Q * H) and (eta / rho)
2	(Q * H) / (eta / rho)

The default reached it on the first try. The extra level bought no accuracy and gave the search room to decorate a finished answer with noise. Raise max_depth when a target you expect is not being recovered, never speculatively: a deeper search enlarges the candidate pool, every threshold scales as 1/m, and each extra candidate costs power on the ones that matter.

Units are the fix¶

Statistics cannot tell that noise_0 does not belong — as we just saw, it is genuinely associated with the target once attached to the core. Physics can. You cannot add a dimensionless number to a quantity in kg·m/s, and beamfeat rejects such expressions at construction, before they are ever scored.

In [5]:

units = {
    "rho": "kg/m**3", "g": "m/s**2", "Q": "m**3/s", "H": "m",
    "eta": "dimensionless",
    **{f"noise_{index}": "dimensionless" for index in range(8)},
}

gated = BeamFeatRegressor(
    random_state=99, max_depth=3, beam_width=60, units=units
).fit(X, y)

print(gated.equation())
print(f"R^2: {gated.score(X, y):.4f}")


# Formulas that add or subtract a noise column outright, across the whole
# certified set rather than the pruned equation.
ADDED_NOISE = re.compile(r"(noise_\d+ [+-]|[+-] noise_\d+)")


def added_noise(model):
    full = BeamFeatRegressor(**{**model.get_params(), "parsimony": None}).fit(X, y)
    return [f for f in full.formulas() if ADDED_NOISE.search(f)], full.n_features_out_


for label, m in (("no units", deep), ("units", gated)):
    hits, total = added_noise(m)
    print(f"certified formulas adding a noise column, {label:8s}: {len(hits)} of {total}")

units = { "rho": "kg/m**3", "g": "m/s**2", "Q": "m**3/s", "H": "m", "eta": "dimensionless", **{f"noise_{index}": "dimensionless" for index in range(8)}, } gated = BeamFeatRegressor( random_state=99, max_depth=3, beam_width=60, units=units ).fit(X, y) print(gated.equation()) print(f"R^2: {gated.score(X, y):.4f}") # Formulas that add or subtract a noise column outright, across the whole # certified set rather than the pruned equation. ADDED_NOISE = re.compile(r"(noise_\d+ [+-]|[+-] noise_\d+)") def added_noise(model): full = BeamFeatRegressor(**{**model.get_params(), "parsimony": None}).fit(X, y) return [f for f in full.formulas() if ADDED_NOISE.search(f)], full.n_features_out_ for label, m in (("no units", deep), ("units", gated)): hits, total = added_noise(m) print(f"certified formulas adding a noise column, {label:8s}: {len(hits)} of {total}")

UserWarning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are meant to carry data.

y = 9.799*((H / eta) * (Q * rho)) + 279.9   [1 of 60 certified terms; parsimony=None prints all 60]
R^2: 0.9996
[beamfeat] fit: 300 rows x 13 columns, split 150 search / 150 selection

UserWarning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are meant to carry data.

[beamfeat] search: depth 3, beam 60 -> 42413 proposed, 60 candidates (3.69s)
[beamfeat] screen: 60 candidates, permutation (BY) at q=0.1 -> 60 certified
[beamfeat] parsimony: off, keeping all 60 certified terms
[beamfeat] result set: 60 terms in the fitted equation, the screened set entire; the guarantee covers it
[beamfeat] result: y = 7.803*(((H / eta) * (Q * rho)) - ((H + noise_0) / (Q - eta))) - 0.7173*(((H + noise_0) * (Q * rho)) + ((Q * rho) * (H - noise_2))) - 0.5565*(((H + noise_7) * (Q * rho)) + ((Q * rho) * (H - noise_3))) + 0.4962*(((H + noise_6) * (Q * rho)) + ((Q * rho) * (H - noise_5))) + 0.5378*(((Q * rho) * (H - noise_2)) + ((Q * rho) * (H - noise_3))) + 0.484*(((H + noise_7) * (Q * rho)) / eta) + 0.1134*(log(H) * ((H + noise_1) * (Q * rho))) + 595.7*(((H / eta) * (H * Q)) / (H + noise_3)) + 2006*(((H * Q) / (H + noise_4)) * (log(Q) * (H * Q))) + 0.2827*(((H + noise_4) * (Q * rho)) + ((Q * rho) * (H - noise_6))) + 0.3142*(((H + noise_0) * (Q * rho)) / eta) + 0.2031*(((H + noise_2) * (Q * rho)) + ((H + noise_5) * (Q * rho))) + 0.3047*(((Q * rho) * (H - noise_1)) + ((Q * rho) * (H - noise_2))) - 479.4*((H + noise_2) * ((H * Q) / (H + noise_0))) + 0.2198*(((H + noise_3) * (Q * rho)) + ((Q * rho) * (H - noise_1))) + 0.8161*(((H * Q) * (Q * rho)) + ((H * Q) * (H - noise_7))) + 0.2482*(((Q * rho) * (H - noise_3)) + ((Q * rho) * (H - noise_4))) - 6049*(((H * Q) / (H + noise_6)) * sqrt((H * Q))) - 1.083*(((H + noise_4) / (Q - eta)) * (log(Q) * (H * Q))) + 57.37*(((H * eta) / (Q - eta)) - (H * Q)) - 2.837*((H / rho) * ((H + noise_4) * (Q * rho))) + 283.4*(((H / eta) * (H * Q)) / (H + noise_0)) + 276.3*((H + noise_5) * ((H * Q) / (H + noise_6))) - 4.504*(((H + noise_2) / (Q - eta)) * sqrt((H * Q))) + 4.392*((log(H) * (Q * rho)) - ((H + noise_3) / (Q - eta))) + 1692*(((H + noise_2) / (rho / Q)) * (H / eta)) - 3.198e-06*(((H + noise_3) * (Q * rho)))^2 - 6.063*(((H * Q) / (H + noise_6)) * ((H / eta) * (H * Q))) + 3.887*(((H / eta) * (H * Q)) - ((H * Q) * (H - noise_3))) + 3.17e-06*(((H / eta) * (Q * rho)))^2 - 0.001444*((H + noise_6) * ((H / eta) * (Q * rho))) - 0.0001463*(rho * ((H + noise_6) * (Q * rho))) + 99.94*(((H * Q) / (H + noise_6)) * ((H * eta) / (Q - eta))) - 1.507e-06*(((H + noise_0) * (Q * rho)) * (H * rho)) + 231.1*(((H / rho) + (Q)^2) * (H + noise_5)) - 4094*(((H * Q) / (H + noise_0)) * log(H)) - 0.003781*(((H + noise_5) * (Q * rho)) * (log(Q) * (H * Q))) - 0.1996*(((H * Q) / (H + noise_7)) * (H * rho)) + 143*((H + noise_3) * ((H * Q) / (H + noise_6))) - 0.0935*(((Q * rho) * (H - noise_1)) + ((Q * rho) * (H - noise_3))) - 1.846e-06*(((H + noise_6) * (Q * rho)) * ((H + noise_7) * (Q * rho))) - 4.525*((Q * rho) + ((Q * noise_5) * (H - noise_5))) - 952.9*(sqrt((H * Q)) - eta) - 0.0002273*(((H + noise_1) / (Q - eta)) * ((H + noise_4) * (Q * rho))) - 0.0002061*(((H + noise_0) / (Q - eta)) * ((H + noise_3) * (Q * rho))) + 2.854e+04*(((H / rho) / log(Q)) - ((H + noise_6) / (rho / Q))) - 1.313e-06*(((H + noise_4) * (Q * rho)) * ((H + noise_5) * (Q * rho))) - 1.662*(((Q * noise_6) * (H - noise_6)) + ((H * Q))^2) + 0.04865*(((Q * rho) * (H - noise_1)) + ((Q * rho) * (H - noise_5))) + 2.329*((Q * rho) + ((Q * noise_3) * (H - noise_3))) + 120.2*(((Q)^2 * (H - noise_2)) + sqrt((H * Q))) + 4.57*(((H * Q) / (Q - eta)) / ((H / rho) + (Q)^2)) + 0.1039*(((H * Q) / (H + noise_6)) * ((H + noise_5) * (Q * rho))) + 42.18*(((H + noise_6) * (Q * rho)) / (H + rho)) - 0.3636*(((H / eta) * (H * Q)) - (noise_4 / Q)) + 0.000298*((H + noise_7) * ((H + noise_2) * (Q * rho))) + 1.004*((Q * rho) + ((Q * noise_6) * (H - noise_6))) - 1.346*sqrt(((H + noise_3) * (Q * rho))) + 1.764e-05*((H + noise_7) * ((H + noise_3) * (Q * rho))) - 0.001119*(((H + noise_4) * (Q * rho)) / eta) + 664.5  [FDR controlled over the printed equation]
certified formulas adding a noise column, no units: 55 of 60

UserWarning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are meant to carry data.

certified formulas adding a noise column, units   : 0 of 60

The same over-deep search now returns the clean core. Identical R^2, no parasites — and the effect is not just in the pruned equation: across the whole certified set, not one formula adds a noise column to a dimensioned quantity, against almost all of them before. That is the general lesson — where a statistical criterion cannot separate two candidates, a structural constraint often can, and it costs nothing because it applies before any numerical work.

Two cautions about what units check.

They enforce that an expression is internally coherent, not that it matches the target's dimension. rho * Q * H / eta is kg·m/s while power is watts, and nothing objects, because g was absorbed. Units reject nonsense; they will not tell you a factor is missing.

And a dimensionless column is legal almost everywhere. Q / noise_7 is still m³/s, so noise can ride in multiplicatively even under a full labelling, and some of the certified set still mentions one. What units removed is the whole class of violations — adding a bare number to a physical quantity — not every appearance of an irrelevant column.

Cover every column¶

Labelling the five real measurements and leaving the noise columns blank looks like the careful thing to do. It is not.

In [6]:

partial = {
    "rho": "kg/m**3", "g": "m/s**2", "Q": "m**3/s", "H": "m",
    "eta": "dimensionless",
}

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    half_gated = BeamFeatRegressor(
        random_state=99, max_depth=3, beam_width=60, units=partial
    ).fit(X, y)

for entry in caught:
    if "units cover" in str(entry.message):
        print(str(entry.message)[:200])
        break

print()
print(half_gated.equation())

partial = { "rho": "kg/m**3", "g": "m/s**2", "Q": "m**3/s", "H": "m", "eta": "dimensionless", } with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always") half_gated = BeamFeatRegressor( random_state=99, max_depth=3, beam_width=60, units=partial ).fit(X, y) for entry in caught: if "units cover" in str(entry.message): print(str(entry.message)[:200]) break print() print(half_gated.equation())

beamfeat: units cover 5 of 13 columns; 8 are unlabelled (noise_0, noise_1, noise_2, noise_3, noise_4, ...) and are treated as dimensionally unconstrained, so expressions combining them with labelled c

y = 9.799*(noise_0 + ((H / eta) * (Q * rho))) + 16.98   [1 of 60 certified terms; parsimony=None prints all 60]

noise_0 is back, added directly to a kg·m/s quantity — the exact violation the units pass exists to catch. An unlabelled column is dimensionally unconstrained, so it combines freely with everything, and the gate stops binding on precisely the columns you did not vouch for. Units are all or nothing; give the genuinely unitless columns "dimensionless".

The equation is not the certified set¶

equation() prints the parsimony subset, and its suffix says so. The guarantee covers the larger screened set, which selection_report_ holds in full, with an exact p- and q-value per candidate.

In [7]:

report = deep.selection_report_
screened = [row for row in report if row["screened"]]
printed = [row for row in report if row["kept"]]

print(f"candidates screened: {len(report)}")
print(f"certified (guarantee applies here): {len(screened)}")
print(f"terms in equation(): {len(printed)}")
print(f"|S|/|S'|: {deep.fdp_inflation_:.2f}")
print(f"fdr_controlled_: {deep.fdr_controlled_}   fdr_scope_: {deep.fdr_scope_!r}")
print()
for row in sorted(screened, key=lambda r: r["p_value"])[:3]:
    formula = row["formula"]
    shown = formula if len(formula) <= 62 else formula[:59] + "..."
    print(f"  {shown:62} p={row['p_value']:.2e} q={row['q_value']:.2e}")

report = deep.selection_report_ screened = [row for row in report if row["screened"]] printed = [row for row in report if row["kept"]] print(f"candidates screened: {len(report)}") print(f"certified (guarantee applies here): {len(screened)}") print(f"terms in equation(): {len(printed)}") print(f"|S|/|S'|: {deep.fdp_inflation_:.2f}") print(f"fdr_controlled_: {deep.fdr_controlled_} fdr_scope_: {deep.fdr_scope_!r}") print() for row in sorted(screened, key=lambda r: r["p_value"])[:3]: formula = row["formula"] shown = formula if len(formula) <= 62 else formula[:59] + "..." print(f" {shown:62} p={row['p_value']:.2e} q={row['q_value']:.2e}")

candidates screened: 60
certified (guarantee applies here): 60
terms in equation(): 1
|S|/|S'|: 60.00
fdr_controlled_: True   fdr_scope_: 'screened set'

  (((H / eta) * (Q * rho)) - ((H + noise_0) / (Q - eta)))        p=1.78e-04 q=8.33e-04
  (((H + noise_2) * (Q * rho)) + ((H + noise_5) * (Q * rho)))    p=1.78e-04 q=8.33e-04
  (((H * Q) * (Q * rho)) + ((H * Q) * (H - noise_7)))            p=1.78e-04 q=8.33e-04

The gap matters. Parsimony picks a subset by greedy forward selection on the same rows, and a data-dependent subset of an FDR-controlled set does not inherit the guarantee: pruning cannot add false selections, but it can raise the false discovery proportion, because the denominator shrinks faster than the numerator, and fdp_inflation_ reports by how much. This is why fdr_controlled_ and fdr_scope_ are two attributes rather than one: the flag is True and the scope is the screened set, not the printed equation. What that costs in fit is small and has been measured — about a thousandth of held-out R^2 across 308 paired fits, in benchmarks/PARSIMONY_COST.md. Report the screened set when you need the guarantee; use the equation when you need something to read.

When nothing is selected¶

Two very different situations produce an empty selection, and they have different fixes. Here is the honest one — a target that is pure noise.

In [8]:

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    empty = BeamFeatRegressor(random_state=0).fit(X, rng.normal(0, 1, n))

print(f"features constructed: {empty.n_features_out_}")
print(f"fdr_controlled_: {empty.fdr_controlled_}")
for entry in caught:
    print(f"  warning: {str(entry.message)[:130]}")

with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always") empty = BeamFeatRegressor(random_state=0).fit(X, rng.normal(0, 1, n)) print(f"features constructed: {empty.n_features_out_}") print(f"fdr_controlled_: {empty.fdr_controlled_}") for entry in caught: print(f" warning: {str(entry.message)[:130]}")

features constructed: 0
fdr_controlled_: False
  warning: beamfeat: 1 of 13 columns are constant (g) and cannot be selected. They are ignored and cost nothing, but check whether they are m
  warning: beamfeat: no feature passed FDR-controlled selection at target FDR 0.1 on the held-out split (on_no_discoveries='empty'). Returnin

An intercept-only model and a visible warning, which is the correct answer: there was nothing to find, and the constructor declined to invent it. fdr_controlled_ reads False here rather than True — no feature was returned, so nothing carries the guarantee, and the flag describes what came back rather than how hard the procedure tried.

The other case looks identical from the outside but is not. If the candidate pool is large enough that the multiplicity threshold falls below the smallest p-value the permutation budget can produce, the correction cannot reject anything, whatever the data say. The selector detects this and says so in its warning, naming the correction and what to change. So read the warning before concluding there is no signal: an empty selection with a budget warning means the test could not fire, not that the data are silent.

A checklist¶

• Start at the defaults. Raise beam_width when an expected target is missed; raise max_depth only when the target genuinely needs the depth, and count composition levels, not variables.
• An equation that grew noise terms is a sign you searched too deep, not that you need to search deeper.
• fdr_controlled_ = True says the formula is not noise. It does not say every column in it belongs.
• Supply units when you have them, for every column.
• A missing constant is usually in the coefficient.
• Quote selection_report_ when you need the guarantee; quote equation() when you need a sentence.
• Read the warnings. They distinguish states that look identical from outside.

What to read next¶

• 02: Search and scoring — beam width, depth and the scoring strategies, with the trace that shows where the cost went.
• 03: Selection and units — the calibration behind the guarantee, why the default is not knockoffs, and the dimensional analysis used above.