simple-scaling-laws

Warning: internal research tooling, largely AI-generated. This package was written primarily by an AI coding agent to support internal experiments in the McDermott Health AI Lab. It is public so it can be installed and cited easily, not because it is a hardened, externally supported library. Read the statistics section before you trust a number from it, check the warnings it emits, and treat its extrapolations as triage signal rather than evidence.

Fit empirical scaling laws to repeated ML training/evaluation records, with conservative uncertainty, and predict performance at model or dataset sizes you have not run yet.

It exists to answer one question for an automated experiment platform: given a handful of models trained at different scales, is this hypothesis worth pursuing at a larger scale? Answering it well requires being honest about how little a small experiment tells you, so the package is opinionated about the statistics and quiet about the knobs.

Install

uv add simple-scaling-laws # or: uv pip install simple-scaling-laws

Quick start

scaling-laws fit runs.parquet --law separable-power --output experiment.slaw

scaling-laws predict experiment.slaw points.parquet --output predictions.parquet

The Python API is the same thing, and is a first-class interface rather than a wrapper:

>>> from simple_scaling_laws import fit
>>> from simple_scaling_laws.simulate import simulate_runs
>>> runs = simulate_runs(
...     {
...         "test_loss__cross_entropy": {"E": 1.0, "A": 2.0, "alpha": 0.30, "B": 1.5, "beta": 0.25},
...         "test_metric__auroc": {"E": 0.92, "A": -0.15, "alpha": 0.40, "B": -0.08, "beta": 0.30},
...     },
...     runs_per_config=2,
...     evaluations_per_run=5,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
... )
>>> model = fit(runs, law="separable-power", n_draws=200, seed=0)
>>> model.targets
('test_loss__cross_entropy', 'test_metric__auroc')

simulate_runs generated that table from a known law, so we can check the fit recovers it:

>>> {k: round(v, 2) for k, v in model.params("test_loss__cross_entropy").items()}
{'E': 1.04, 'A': 1.96, 'alpha': 0.3, 'B': 1.5, 'beta': 0.25}

Prediction returns a dataframe, so plotting is whatever you already use:

>>> points = pl.DataFrame(
...     {
...         "point_id": ["current", "10x-model", "10x-data"],
...         "model_size__n_params": [1e8, 1e9, 1e8],
...         "dataset_size__n_subjects": [1e5, 1e5, 1e6],
...     }
... )
>>> predictions = model.predict(points, quantiles=[0.025, 0.5, 0.975])
>>> for row in predictions.iter_rows(named=True):
...     print(
...         f"{row['point_id']:<10} {row['domain']:<14} "
...         f"{row['test_loss__cross_entropy__median']:.3f} "
...         f"[{row['test_loss__cross_entropy__q025']:.3f}, "
...         f"{row['test_loss__cross_entropy__q975']:.3f}]"
...     )
current    interpolation  2.856 [2.838, 2.878]
10x-model  extrapolation  2.367 [2.323, 2.413]
10x-data   extrapolation  2.487 [2.456, 2.530]

Input format

One table, long over evaluation resamples and wide over metrics. Column roles are inferred from prefixes, so a conventionally named table needs no configuration at all:

Every target column found is fit, each with its own parameters and its own uncertainty. Columns that match nothing are ignored, so extra bookkeeping is harmless.

Three structural rules carry all the statistical weight:

1. One training_run_id per independently trained model. Different seeds or different splits mean different runs.
2. Many rows per run, one per evaluation resample. The spread within a run estimates evaluation noise. Do not pre-aggregate into bootstrap quantiles -- give the package the rows.
3. Reuse test_set_id across models. When the same resample scores several models, the package detects the shared noise instead of mistaking it for a real difference between them.

Non-conventional column names can be mapped without renaming anything:

>>> renamed = runs.rename({"training_run_id": "run", "test_set_id": "bootstrap_id"})
>>> remapped = fit(
...     renamed,
...     columns={"training_run_id": "run", "test_set_id": "bootstrap_id"},
...     n_draws=50,
...     seed=0,
... )
>>> remapped.primary_target
'test_loss__cross_entropy'

or in a small YAML config passed to --config:

law: separable-power

columns:
  training_run_id: run
  test_set_id: bootstrap_id

What it fits

Seven law families. Each is linear in its offset and amplitudes, which are solved exactly, and nonlinear in a handful of other parameters -- exponents, and for the richer laws a shift, a ratio or a break -- which are searched:

>>> from simple_scaling_laws import available_laws
>>> for name, description in available_laws().items():
...     print(f"{name}\n    {description}")
separable-power
    E + sum_i A_i * x_i**-alpha_i (one additive power term per predictor)
multiplicative-power
    E + A * prod_i x_i**-alpha_i (one joint multiplicative power term)
shifted-separable-power
    E + sum_i A_i * (x_i + s_i)**-alpha_i (power terms with a finite small-scale limit)
kaplan-joint-power
    E + A * (N**-rho + R / D)**beta (Kaplan et al. joint model/data law)
skaling-joint-power
    E + A * (N**-alpha + R * D**-beta)**k (Videau et al. Skaling law: coupled model/data terms)
bounded-rational-power
    (E + I * g) / (1 + g), g = sum_i A_i * x_i**-alpha_i (plateau I at small scale, floor E at large)
broken-power
    E + A * x**-c0 * (1 + (x / d)**(1 / f))**(-c1 * f) (one smooth change of slope; single predictor)

With one model-size and one dataset-size column, separable-power is the conventional E + A·N^-α + B·D^-β. Start there. The other families each relax one specific assumption it makes, at the price of more parameters and a larger grid, and each except multiplicative-power (a single joint power term, E + A·N^-α·D^-β) has a section below saying when it is worth that price. A lower in-sample residual is not by itself a reason to prefer one: every richer law here can fit separable data at least as well as the separable law does.

Losses are fit with non-negative amplitudes, which gives the usual "irreducible floor plus decaying terms" reading. Metrics get free-signed amplitudes so an increasing metric such as AUROC can approach its asymptote from below:

>>> auroc = model.params("test_metric__auroc")
>>> auroc["A"] < 0 and auroc["B"] < 0
True

The examples below

The richer laws are shown on a 4x4 grid, and some are compared against data that follow the plain separable law exactly:

>>> grid = {"model_sizes": (1e6, 1e7, 1e8, 1e9), "dataset_sizes": (1e3, 1e4, 1e5, 1e6)}
>>> separable = simulate_runs(
...     {"test_loss__cross_entropy": {"E": 1.0, "A": 2.0, "alpha": 0.3, "B": 1.5, "beta": 0.25}},
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
...     **grid,
... )

There is deliberately no law here in which a larger model does worse when data are scarce. The obvious form -- the separable law minus a joint term C·N^-γ·D^-δ -- was tried and removed. Any C large enough to couple the two resources also makes the loss rise with model size somewhere: always at small enough data, and at every data size when γ < α. So where the curve turns up is set by the formula, not by the data. Videau et al. (2026) reject the form on the same grounds. The two joint laws below couple model size and data while staying monotone in both. For a curve that genuinely turns, use broken-power along one axis.

When the small end is not yet a power law

A power law is a straight line on log-log axes all the way down, and diverges at zero scale. Real curves flatten toward something finite at the small end of an experiment, and those pre-asymptotic points drag a plain power law's exponent toward zero. shifted-separable-power puts a non-negative shift inside each power term, E + A·(N + s_N)^-α + B·(D + s_D)^-β, which gives the term a finite value at zero scale and leaves its large-scale behavior alone. Its single-predictor form is the saturating law a·(C + d)^-b + c of Zhai et al. (2022). It is not the "M3" form of Bansal et al. (2022), α·(D^-1 + C)^p, which shifts the reciprocal of the scale and bounds the large end instead.

The cost of ignoring a shift that is there is not subtle. Fitting the plain law to shifted data:

>>> truth = {"E": 1.0, "A": 2.0, "B": 1.5, "alpha": 0.5, "beta": 0.4, "s_n_params": 0.3, "s_n_subjects": 0.2}
>>> bent = simulate_runs(
...     {"test_loss__cross_entropy": truth},
...     law="shifted-separable-power",
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
...     **grid,
... )
>>> plain = fit(bent, law="separable-power", n_draws=200, seed=0)
>>> shifted = fit(bent, law="shifted-separable-power", n_draws=200, seed=0)
>>> {k: round(v, 2) for k, v in plain.params("test_loss__cross_entropy").items()}
{'E': -15.81, 'A': 12.84, 'alpha': 0.04, 'B': 7.21, 'beta': 0.05}
>>> [note.code for note in plain.warnings]
['lack_of_fit']
>>> {k: round(v, 2) for k, v in shifted.params("test_loss__cross_entropy").items()}
{'E': 1.0, 'A': 1.99, 'B': 1.51, 'alpha': 0.5, 'beta': 0.39, 's_n_params': 0.3, 's_n_subjects': 0.19}

Two decades beyond the grid, where the true loss is 1.10, the plain law predicts a negative one and its interval does not reach the truth; the shifted law's does:

>>> far = {"model_size__n_params": [1e11], "dataset_size__n_subjects": [1e8]}
>>> for name, candidate in (("plain", plain), ("shifted", shifted)):
...     row = candidate.predict(far).row(0, named=True)
...     print(
...         f"{name:<8} {row['test_loss__cross_entropy__median']:6.2f} "
...         f"[{row['test_loss__cross_entropy__q025']:.2f}, {row['test_loss__cross_entropy__q975']:.2f}]"
...     )
plain     -1.33 [-2.18, -0.27]
shifted    1.10 [1.05, 1.15]

lack_of_fit on a separable fit is therefore the cue to try this law. It needs four scales of every predictor rather than three -- with three there are two contrasts from which to learn an amplitude, an exponent and a shift -- and says so with an error-severity too_few_predictor_levels note if it does not get them. On data that are already a clean power law the shifts sit at zero and are reported as parameter_at_bound. Measured on 4x4 and 5x5 designs, its nominal 95% prediction intervals covered the truth 88-98% of the time, inside the grid and a decade beyond it.

When one resource is the bottleneck

kaplan-joint-power is the joint law of Kaplan et al. (2020), E + A·(N^-ρ + R/D)^β, with an irreducible floor E added. Whichever resource is scarcer sets the loss: with ample data it is a power law in model size (exponent ρ·β), with an ample model it is a power law in data (exponent β), and in between, adding more of the abundant resource does almost nothing.

>>> truth = {"E": 1.0, "A": 2.0, "rho": 0.5, "beta": 0.6, "R": 1.0}
>>> limited = simulate_runs(
...     {"test_loss__cross_entropy": truth},
...     law="kaplan-joint-power",
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
...     **grid,
... )
>>> kaplan = fit(limited, law="kaplan-joint-power", n_draws=200, seed=0)
>>> {k: round(v, 2) for k, v in kaplan.params("test_loss__cross_entropy").items()}
{'E': 1.01, 'A': 2.0, 'rho': 0.5, 'beta': 0.6, 'R': 1.0}
>>> sweep = kaplan.predict(
...     {
...         "model_size__n_params": [1e6, 1e7, 1e8, 1e9, 1e10, 1e11] * 2,
...         "dataset_size__n_subjects": [1e3] * 6 + [1e6] * 6,
...     }
... )
>>> losses = [round(v, 2) for v in sweep["test_loss__cross_entropy__median"]]
>>> losses[:6]  # scarce data: a hundred thousand times the model buys almost nothing
[18.53, 17.42, 17.06, 16.94, 16.91, 16.89]
>>> losses[6:]  # ample data: the model is the bottleneck, and scaling it pays
[6.65, 3.86, 2.47, 1.79, 1.47, 1.34]

It is monotone in both resources, so it describes saturation, not a larger model being worse. It is defined for exactly one model-size and one dataset-size column, and both must actually vary.

Measured on 4x4 designs, the Kaplan law's nominal 95% prediction intervals covered the truth 92-100% of the time, inside the grid and a decade beyond it. On a 3x3 grid with two runs per configuration the Kaplan law's covered 93-95%, but only because the bootstrap puts back noise that a fit this close to saturated absorbs (78-85% without it); prefer the larger grid.

When model size and data are coupled, but neither is ever harmful

skaling-joint-power is the Skaling law of Videau, Youbi-Idrissi, Lopez-Paz & Ahuja (2026), E + A·(N^-α + R·D^-β)^k: the two power terms of the separable law, summed and raised to one coupling exponent k. At k = 1 it is exactly separable-power; with the inner data exponent fixed at one it is kaplan-joint-power. For 0 < k < 1 each resource makes the other more valuable -- the coupling the authors measured on large language-model grids, with k from about 0.3 to 0.9 -- and the loss still falls in both. With ample data it is a power law in model size with exponent α·k, with an ample model a power law in data with exponent β·k.

>>> truth = {"E": 1.0, "A": 2.0, "alpha": 0.5, "beta": 0.4, "k": 0.5, "R": 1.0}
>>> joint = simulate_runs(
...     {"test_loss__cross_entropy": truth},
...     law="skaling-joint-power",
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
...     **grid,
... )
>>> skaling = fit(joint, law="skaling-joint-power", n_draws=200, seed=0)
>>> {k: round(v, 2) for k, v in skaling.params("test_loss__cross_entropy").items()}
{'E': 1.03, 'A': 1.98, 'alpha': 0.5, 'beta': 0.4, 'k': 0.5, 'R': 0.99}

On the separable data from the section above, which have no coupling at all, the fit puts k at one and says nothing, because the separable law it then reduces to is the right answer. Its A and R·A are then that law's two amplitudes:

>>> skaling_on_plain = fit(separable, law="skaling-joint-power", n_draws=200, seed=0)
>>> round(skaling_on_plain.params("test_loss__cross_entropy")["k"], 2)
1.0
>>> [note.code for note in skaling_on_plain.warnings]
[]

It is defined for exactly one model-size and one dataset-size column. The paper fit 117 to 404 runs per grid; on simulated 3x3 grids with two runs per configuration k still came back within about 0.01 of the truth and nominal 95% prediction intervals covered 88-100%.

When the experiment reaches down to chance

Every law above grows without limit as scale shrinks. A real loss does not: a model trained on almost nothing sits at chance and stays there. bounded-rational-power has a plateau I at the small end as well as a floor E at the large one, and moves between them as a weighted average, (E + I·g) / (1 + g), whose weight g = A·N^-α + B·D^-β is an ordinary separable sum of power terms. It is the M4 estimator of Alabdulmohsin, Neyshabur & Zhai (2022) with its shape exponent (the paper's α, in (f - E) / (I - f)^α = g) fixed at one -- the one case that can be written down in closed form, and so the one whose two limits the fitter can solve exactly rather than search for.

>>> truth = {"E": 1.0, "I": 4.0, "A": 1.0, "B": 0.8, "alpha": 0.5, "beta": 0.4}
>>> saturating = simulate_runs(
...     {"test_loss__cross_entropy": truth},
...     law="bounded-rational-power",
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
...     **grid,
... )
>>> bounded = fit(saturating, law="bounded-rational-power", n_draws=200, seed=0)
>>> {k: round(v, 2) for k, v in bounded.params("test_loss__cross_entropy").items()}
{'E': 1.03, 'I': 4.0, 'A': 1.01, 'B': 0.78, 'alpha': 0.49, 'beta': 0.41}

What it buys is a sane answer below the grid. Three decades under the smallest configuration the true loss is 3.99, just short of the plateau. The separable law, fit to the same runs, has nothing to stop it:

>>> tiny = {"model_size__n_params": [1e3], "dataset_size__n_subjects": [1e0]}
>>> for name, candidate in (("bounded", bounded), ("separable", fit(saturating, n_draws=200, seed=0))):
...     row = candidate.predict(tiny).row(0, named=True)
...     print(
...         f"{name:<10} {row['test_loss__cross_entropy__median']:5.2f} "
...         f"[{row['test_loss__cross_entropy__q025']:.2f}, {row['test_loss__cross_entropy__q975']:.2f}]"
...     )
bounded     3.99 [3.95, 4.03]
separable   8.14 [5.06, 144.22]

The plateau is only identified if the experiment reaches far enough down to see the curve bend toward it. On data that are still a clean power law I runs off to its bound while the weights shrink to compensate, and the fit is reported as ill_conditioned_fit: use separable-power. Both limits are estimated, never assumed; for a metric confined to [0, 1] they are limits on the logit scale the metric is fit on. Its intervals covered the truth 82-98% of the time on 4x4 designs and 97-98% on a 3x3 grid with two runs per configuration -- there only because the bootstrap puts back noise the fit absorbed (75-82% without it), so this law too wants the larger grid.

On the plain separable data above, the note names the parameters that trade off against each other with almost no change in the fitted curve -- read it as the instruction to drop back to separable-power:

>>> unbent = fit(separable, law="bounded-rational-power", n_draws=200, seed=0)
>>> next(note for note in unbent.warnings if note.code == "ill_conditioned_fit").details["parameters"]
['I', 'A', 'B']

When the slope changes partway

Every law above has one slope on log-log axes. broken-power has two, joined smoothly at a break scale d: E + A·x^-c0·(1 + (x/d)^(1/f))^(-c1·f), the one-break Broken Neural Scaling Law of Caballero et al. (2022). Below the break the excess loss falls as x^-c0, above it as x^-(c0+c1), and c1 may have either sign -- so the curve can steepen, flatten, or, when c1 < -c0, turn back up. It is a law for a sweep along one predictor:

>>> truth = {"E": 1.0, "A": 2.0, "c0": 0.3, "c1": -0.6, "d": 3.0, "f": 0.8}
>>> sizes = tuple(10.0**k for k in range(6, 14))
>>> sweep = simulate_runs(
...     {"test_loss__cross_entropy": truth},
...     model_sizes=sizes,
...     dataset_sizes=(),
...     law="broken-power",
...     runs_per_config=2,
...     evaluations_per_run=4,
...     run_sd=0.01,
...     eval_sd=0.02,
...     seed=0,
... )
>>> broken = fit(sweep, law="broken-power", n_draws=200, seed=0)
>>> {k: round(v, 2) for k, v in broken.params("test_loss__cross_entropy").items()}
{'E': 0.99, 'A': 2.0, 'c0': 0.3, 'c1': -0.6, 'd': 2.98, 'f': 0.8}
>>> at_sizes = {"model_size__n_params": list(sizes)}
>>> [round(v, 1) for v in broken.predict(at_sizes)["test_loss__cross_entropy__median"]]
[23.4, 12.2, 6.6, 3.9, 3.0, 4.0, 6.8, 12.6]

A single power law cannot turn, so on the same runs it settles for a floor through the middle of them, and says the fit is poor:

>>> single = fit(sweep, n_draws=200, seed=0)
>>> [round(v, 1) for v in single.predict(at_sizes)["test_loss__cross_entropy__median"]]
[23.3, 11.8, 8.6, 7.3, 6.8, 6.2, 6.2, 6.2]
>>> "lack_of_fit" in [note.code for note in single.warnings]
True

Three things to know before using it:

• It needs at least six scales of its predictor -- it has six parameters on one axis -- and raises too_few_predictor_levels otherwise.
• The break can only be found inside the range that was run. A break estimated at either end of it is reported as parameter_at_bound, and means no break was observed. An extrapolation simply continues the upper slope; it cannot anticipate another change of regime beyond the grid.
• Its intervals rely on knowing the noise. With six parameters on eight scales the fit passes almost exactly through some of them, and the bootstrap puts back the noise the fit absorbed there from the separately estimated run and evaluation variances. With that, nominal 95% prediction intervals covered 90-97% in simulation; without it they covered 70-88%. Replicate runs at a few scales make that noise estimate direct rather than inferred from the residuals.

Bounded metrics

AUROC lives in [0, 1]; a fitted asymptote of 4.9 is not a statement about a classifier. Metrics known to be confined to the unit interval are therefore fit on their logit and mapped back for prediction, which makes the bound hold by construction rather than by clipping afterwards. It also makes the constant-variance assumption far more plausible, since an AUROC near 0.99 has much less room to move than one near 0.7.

The common metrics are recognized by name, with no configuration:

>>> from simple_scaling_laws.metrics import describe
>>> for name in ("auroc", "auprc", "accuracy", "brier", "ece", "cross_entropy"):
...     kind = describe(name, is_loss=name == "cross_entropy")
...     print(f"{name:<14} {kind.support:<5} {kind.direction} is better")
auroc          unit  higher is better
auprc          unit  higher is better
accuracy       unit  higher is better
brier          unit  lower is better
ece            unit  lower is better
cross_entropy  real  lower is better

Direction matters because the role prefix does not carry it: test_metric__auroc improves upward and test_metric__brier improves downward, and only the direction turns a difference between two systems into a win rate. Anything the table does not recognize is left untransformed, with its direction marked as assumed rather than known. Describe such a target explicitly:

>>> described = fit(
...     runs,
...     targets={"test_metric__auroc": {"support": "unit", "direction": "higher"}},
...     n_draws=50,
...     seed=0,
... )
>>> described.target_kind("test_metric__auroc").support
'unit'

The same goes in a config file under a targets: key. Note that a transformed target's parameters are on the logit scale — its E is an asymptote in logit units — while everything predict returns is on the metric's own scale, and cannot leave it:

>>> auroc = described.predict(
...     {
...         "model_size__n_params": [1e6, 1e20],
...         "dataset_size__n_subjects": [1e3, 1e20],
...     }
... )
>>> lo, hi = auroc["test_metric__auroc__q025"], auroc["test_metric__auroc__q975"]
>>> bool((lo >= 0).all() and (hi <= 1).all())
True

Comparing systems

Fitting a law per system says how each one scales. The question a triage platform actually asks is comparative — of these arms, which is worth scaling up? — and two overlapping intervals are perfectly compatible with one arm being reliably better. What you need is the distribution of the difference:

scaling-laws compare points.parquet \
	--model baseline=baseline.slaw \
	--model tuned=tuned.slaw \
	--reference baseline

>>> from simple_scaling_laws import compare
>>> def arm(alpha, seed):
...     return simulate_runs(
...         {"test_loss__cross_entropy": {"E": 1.0, "A": 2.0, "alpha": alpha, "B": 1.5, "beta": 0.25}},
...         runs_per_config=2,
...         evaluations_per_run=4,
...         run_sd=0.01,
...         eval_sd=0.02,
...         seed=seed,
...     )
>>> arms = {
...     "baseline": fit(arm(0.30, 0), n_draws=300, seed=0),
...     "tuned": fit(arm(0.34, 1), n_draws=300, seed=0),
...     "risky": fit(arm(0.28, 2), n_draws=300, seed=0),
... }
>>> table = compare(
...     arms,
...     {
...         "point_id": ["now", "10x"],
...         "model_size__n_params": [1e8, 1e9],
...         "dataset_size__n_subjects": [1e5, 1e5],
...     },
...     reference="baseline",
... )
>>> for row in table.iter_rows(named=True):
...     print(
...         f"{row['point_id']:<4} {row['system']:<9} "
...         f"{row['value__median']:.3f}  best {row['p_best']:>4.0%}  "
...         f"better than baseline {row['p_better_than_reference']:>4.0%}"
...     )
10x  baseline  2.352  best   0%  better than baseline   0%
10x  risky     2.382  best   0%  better than baseline  29%
10x  tuned     2.244  best 100%  better than baseline 100%
now  baseline  2.848  best   0%  better than baseline   0%
now  risky     2.891  best   0%  better than baseline   3%
now  tuned     2.747  best 100%  better than baseline 100%

Any number of arms works, and the output is long — one row per point per arm — so ten arms do not become a hundred columns. p_best is the probability that arm is the best of the set at that point; with a reference, each arm also gets the distribution of its difference from the reference (difference__median, difference__q025, …) and the probability it beats it. Direction is handled for you: a loss is better when smaller, AUROC when larger.

Why the draws have to be paired

Arms run at the same scales share a configuration-level deviation — a split that is hard for one is hard for the other, and where the law misfits it bends both curves the same way. That shared part cancels in the difference, so the difference is pinned down far more sharply than either level is.

The wild bootstrap makes this nearly free. It perturbs residuals with a multiplier drawn once per configuration; if every arm gets the same multiplier sequence, the perturbation of the difference is driven by the difference's own residual. That is right whether or not the arms share splits: when their residuals are correlated the difference is sharp, and when they are independent the arithmetic recovers the independent answer on its own. This is common random numbers, the standard device for comparing simulations, and it is why a paired t-test beats an unpaired one.

It matters a great deal. Measured against the true sampling spread of the difference, over 120 simulated worlds:

Differencing two independently-seeded fits is 4x to 11x too wide once the arms share structure — wide enough that no real difference is ever detectable. Fits pair automatically when they were run with the same seed on the same configuration grid, which is the normal case for arms of one experiment; compare verifies it, reports paired, and warns loudly when it cannot.

The statistics, briefly

Evaluation rows are not training evidence. Each run is reduced to one observation -- its mean -- plus a count of how many evaluations that mean averages. Scoring one model on a hundred bootstrap resamples makes that model's mean more precise; it does not make it a hundred models.

Two variance components. Spread within a run estimates finite-evaluation-set noise; spread between independently trained models at the same scale estimates training-run stochasticity. The run-level variance is pooled across configurations rather than estimated separately from the two or three repeats you have at each. Without replicate runs it falls back to the scatter around the fitted curve, which also contains the law's misspecification -- biased upward, which is the safe direction.

Shared test sets are recognized. When models share test_set_id values, the correlation between their evaluation residuals is estimated and reported, and only the part of evaluation noise that is not shared is treated as averaging away when runs are compared. This matters more than it sounds: with nine tenths of the evaluation noise coming from the resamples themselves, subtracting the marginal within-run variance drives the estimated run-to-run variance to exactly zero in about half of simulated experiments -- which would make a new-run prediction no wider than a mean one. Subtracting only the unshared part recovers the true value.

Uncertainty is a wild cluster bootstrap over configurations. Residuals around the fitted curve are corrected and re-signed with Webb's six-point weights, drawn once per configuration so that a scale's shared lack of fit is preserved rather than averaged away. This was chosen because it is assumption-free about the error distribution, stays valid with the four-to-a-dozen configurations a real scaling experiment has, keeps the (N, D) design fixed instead of resampling it into degeneracy, and is conservative: residuals carry the law's own misspecification into the intervals. Draws -- not quantiles -- are stored, so any summary can be computed later.

The intervals are calibrated against simulation, not asserted. Resampling residuals naively produces intervals that are badly too narrow at these sample sizes: generating from a known law and counting how often the nominal 95% intervals contained it gave 74% at nine configurations and 87% at sixteen. Two corrections fix it -- a degrees-of-freedom factor for the parameters the fit consumed at the level the bootstrap resamples at, and a Student's-t rather than normal reference for the interval endpoints -- and take measured coverage to 91-95% across the designs this package targets. A third correction matters when a law has nearly as many parameters as configurations: the fit then passes almost exactly through some of them, their residuals understate the noise, and each configuration's pseudo-residual is floored at the share of its (separately estimated) noise that the fit absorbed. It never shrinks a configuration's pseudo-residual, so it pushes intervals wider, not narrower. tests/test_coverage.py re-measures this on every CI run, so it cannot silently regress.

That calibration was measured on separable-power. Every other law is held to the same test (tests/test_law_families.py), and its section above reports what it measured.

Correlations are computed across runs. Loss/metric correlations use run-level means, with intervals from resampling configurations. The naive evaluation-row correlation is reported alongside, with its row count, purely so the inflated sample size is visible -- its estimate can land close to the run-level one while claiming five times the evidence:

>>> correlations = model.metric_correlations()
>>> for row in correlations.iter_rows(named=True):
...     print(
...         f"{row['target']}: run-level r={row['pearson']:.3f} over {row['n_runs']} runs "
...         f"(evaluation-row r={row['evaluation_row_pearson']:.3f} over "
...         f"{row['n_evaluation_rows']} rows)"
...     )
test_metric__auroc: run-level r=-0.979 over 18 runs (evaluation-row r=-0.971 over 90 rows)

Prediction semantics

mean is the default because it is the scaling law itself. new-run adds training-run stochasticity, resampled from the observed run-to-run deviations rather than assumed Gaussian when there are at least ten of them (otherwise drawn as Gaussian with the estimated run SD):

>>> def width(frame, target="test_loss__cross_entropy"):
...     return (frame[f"{target}__q975"] - frame[f"{target}__q025"]).to_list()
>>> mean_widths = width(model.predict(points))
>>> new_run_widths = width(model.predict(points, kind="new-run"))
>>> all(n > m for n, m in zip(new_run_widths, mean_widths, strict=True))
True

Evaluation-set noise is deliberately not included: a scaling prediction is about the model, not about which bootstrap resample you happen to score it on.

Interpolation and extrapolation

Both go through predict. The difference is reported, not enforced -- extrapolating is the entire point of fitting a scaling law. Every prediction carries a domain label and an extrapolation_distance, measured in units of each predictor's observed log range, so one full observed range beyond the edge is a distance of 1.0:

>>> for row in predictions.iter_rows(named=True):
...     print(f"{row['point_id']:<10} {row['domain']:<14} {row['extrapolation_distance']:.2f}")
current    interpolation  0.00
10x-model  extrapolation  0.50
10x-data   extrapolation  0.50

An extrapolated interval reflects parameter uncertainty only. It cannot tell you that the power law stops holding two orders of magnitude out, so a UserWarning accompanies any extrapolated request.

The artifact

model.save(path) writes a self-contained directory that loads without the original dataframe:

>>> artifact = model.save(Path(tempfile.mkdtemp()) / "experiment.slaw")
>>> print_directory(artifact)
├── diagnostics.json
├── draws.npz
├── fits.json
└── manifest.json

• manifest.json -- law, column roles, observed predictor ranges, design counts, provenance.
• fits.json -- point estimates (normalized and raw-scale), variance components, optimizer and goodness-of-fit diagnostics.
• draws.npz -- the uncertainty draws themselves, keyed "<target>/<parameter>", plus run_sd and the observed run-level deviations.
• diagnostics.json -- correlations, cross-target comparisons, fit quality and every warning.

A reloaded model predicts identically, to the bit:

>>> from simple_scaling_laws import load
>>> load(artifact).predict(points).equals(model.predict(points))
True

Warnings

The package returns a fit whenever one is mathematically possible and records its objections instead of raising them, so an automated caller can decide how much to trust a prediction:

>>> thin = simulate_runs(
...     {"test_loss__ce": {"E": 1.0, "A": 2.0, "alpha": 0.3, "B": 1.5, "beta": 0.25}},
...     model_sizes=(1e6, 1e8),
...     dataset_sizes=(1e3, 1e5),
...     runs_per_config=1,
...     evaluations_per_run=4,
...     run_sd=0.02,
...     eval_sd=0.03,
...     seed=0,
... )
>>> for note in fit(thin, n_draws=50, seed=0).warnings:
...     print(note.code, "-", note.severity)
underdetermined - error
too_few_configurations - warning
too_few_predictor_levels - error
too_few_predictor_levels - error
single_run_per_configuration - warning
parametric_uncertainty - warning

Detected conditions include: too few distinct configurations, fewer configurations than parameters, a predictor measured at fewer than three distinct scales, a predictor held fixed, collinear predictors, a single run per configuration, no repeated evaluations, optimizer non-convergence, parameters pinned to a bound, weakly identified exponents, parameters that trade off against each other without changing the curve, lack of fit, a constant target, no uncertainty information at all, and extrapolation at predict time.

The three-scales rule deserves its own note, because it is the one failure widening cannot rescue. Restricted to a single predictor a power-law term is E + A·x^-α: three free parameters, so two scales can be fit exactly by infinitely many different curves. (A law that puts more on a predictor needs more scales of it, and the note says how many: four for shifted-separable-power, six for broken-power.) A design with only two dataset sizes yields parameters that look plausible, intervals that look reasonable, and simulated coverage of 0-38%. The package raises an error-severity too_few_predictor_levels note in 100% of such fits -- it cannot give you the right answer, but it will not pretend to.

Programmatic inspection

>>> {k: round(v, 3) for k, v in model.exponents("test_loss__cross_entropy").items()}
{'alpha': 0.305, 'beta': 0.25}
>>> low, high = model.conf_int("test_loss__cross_entropy")["alpha"]
>>> low < 0.3 < high
True
>>> model.observed_domain["model_size__n_params"]["max"]
100000000.0
>>> model.draws["test_loss__cross_entropy"].params.shape
(200, 5)

model.fits, model.draws, model.diagnostics, model.observed_domain, model.predictors and model.warnings are all exposed. There is no plotting here by design: anything you would want to draw can be built from model.predict().

Development

See CONTRIBUTORS.md. In short:

uv sync
uv run pre-commit install
uv run pytest -v

License

MIT. See LICENSE.