Continuous Learning

Model-free geo response-surface learning: CCD designs, Thompson allocation, funding-line and ENBS stopping, Laplace knowledge-gradient design scoring.

Continuous sequential learning — a model-free geo response-surface bandit.

A self-contained Bayesian sequential-experimentation loop that allocates continuous budget across channels by repeatedly (1) fitting a response surface directly from designed experiment data, (2) choosing the most valuable next experiment, and (3) stopping when further testing no longer pays — without requiring a pre-fit MMM. The experiment’s designed cross-sectional variation identifies the surface, so the priors inform but the data dominates.

This complements (does not replace) the model-anchored planning layer in mmm_framework.planning, which sits on top of a fitted BayesianMMM. Here there is no observational model: the loop learns from experiments alone.

Layers

• surface — the differentiable JAX response surface (single source of truth, shared by the likelihood, the DGP, and the allocator).

• model — the NumPyro generative model + priors + the Posterior container and fit().

• design — central-composite geo cells (central_composite()) and geo assignment (assign_geos()).

• dgp — a synthetic TrueWorld with causal ground truth, the recovery harness simulate_panel(), and the fantasy engine.

• planner — the allocator (allocate_under_sample(), thompson_wave()), the funding line (marginal_roas()), the stopping rule (expected_regret(), enbs()), and decision-aware EVSI (knowledge_gradient()).

• loop — LearningState (carry the posterior across waves) and run_closed_loop() (the end-to-end demo / closure test).

The Hill activation matches SaturationType.HILL (slope = alpha, sat_half = kappa), so a continuous-learning posterior is directly comparable to a BayesianMMM Hill fit on the same channel.

See technical-docs/continuous-learning.md and assets/continous_learning.md.

mmm_framework.continuous_learning.activation(spend, kappa, alpha)

Beta-stripped Hill activation f_c(s_c) in [0, 1), elementwise.

Parameters:

• spend – scaled spend, shape (K,) (non-negative; a shutoff cell is 0).

• kappa – half-saturation per channel, shape (K,) (positive).

• alpha – Hill shape per channel, shape (K,) (positive).

Returns:

The activation fraction, shape (K,).

mmm_framework.continuous_learning.logistic(spend, lam)

Exponential-saturation activation f(s) = 1 - exp(-lam * s) in [0, 1).

Smooth, strictly increasing, strictly concave (no inflection), f(0)=0. One shape parameter per channel; the half-saturation point is ln(2)/lam. Matches the framework’s SaturationType.LOGISTIC.

mmm_framework.continuous_learning.hill_mixture(spend, kappa1, alpha1, kappa2, alpha2, w)

A weighted sum of two Hill curves: w f(κ1,α1) + (1-w) f(κ2,α2).

Still smooth, monotone and saturating (f(0)=0, -> a value in [0,1)), but far more flexible than a single Hill — it can express a two-phase / shoulder shape (a low-κ, high-α component that switches on early plus a high-κ component that keeps rising) that a single Hill or a logistic curve cannot represent. Used as a true DGP to study model misspecification: fit it with a single Hill (mild) or a logistic (severe) and watch what breaks.

mmm_framework.continuous_learning.monotone_spline(spend, w1, w2, w3, w4, w5, w6, w7, w8, w9)

Shape-agnostic monotone-spline activation, values in [0, 1).

f(s) = sum_j w_j I_j(s / MSPLINE_S_MAX) / sum_j w_j with positive weights w1..w9 (one per I-spline basis function; only the ratios matter). Smooth, strictly monotone on [0, MSPLINE_S_MAX), f(0) = 0, fully saturated at/after MSPLINE_S_MAX. Early weights buy an early, steep rise; late weights buy a slow, late one — mixtures express two-phase and plateau shapes no single Hill or logistic can.

mmm_framework.continuous_learning.monotone_spline_basis(spend)

The MSPLINE_J monotone I-spline basis functions at spend.

Returns shape spend.shape + (MSPLINE_J,); column j is I_{j+1}(spend / MSPLINE_S_MAX), rising monotonically 0 -> 1 (earlier columns rise earlier). Useful for plotting and for understanding what the fitted weights mean.

mmm_framework.continuous_learning.incremental(spend, beta, kappa, alpha, gamma)

Incremental response to one scaled spend vector.

Parameters:

• spend – scaled spend, shape (K,).

• beta – channel ceilings, shape (K,) (non-negative).

• kappa – half-saturations, shape (K,).

• alpha – Hill shapes, shape (K,).

• gamma – upper-triangular interaction matrix, shape (K, K) (zeros on and below the diagonal).

Returns:

Scalar incremental response.

mmm_framework.continuous_learning.response_curve(spend_matrix, beta, kappa, alpha, gamma)

Incremental response for every row of a (N, K) scaled-spend panel.

A thin numpy-friendly wrapper over incremental_batch returning a plain jax.Array of shape (N,).

mmm_framework.continuous_learning.surface_value(spend, beta, gamma, act_fn, shape)

Incremental response for one spend vector under an arbitrary activation.

act_fn(spend, *shape) -> f computes the per-channel activations; shape is a tuple of (K,) arrays in the activation’s parameter order. The interaction block is unchanged, so any activation drops straight in.

mmm_framework.continuous_learning.surface_over_rows(spend_matrix, beta, gamma, act_fn, shape)

surface_value() for every row of a (N, K) panel -> (N,).

mmm_framework.continuous_learning.fit(data, *, channels, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', beta_scale=1.0, gamma_scale=0.8, num_warmup=500, num_samples=500, num_chains=2, seed=0, progress_bar=False, spend_ref=None, prior_scaling='auto', discount_half_life=None, spline_prior='iid')

Fit the response surface to a geo-week panel and/or summary readouts.

Parameters:

• data (dict[str, Any]) – the data contract — {"spend": (N, K), "geo_idx": (N,), "n_geo": int, "y": (N,)} (scaled spend, natural-unit y), optionally with "summaries": list[dict] (see the module header), an optional "period_idx": (N,) (0-based int; required when time_effect="national"), or summaries-only ({"summaries": [...], "n_geo": 0} / an empty (0, K) panel). At least one of panel rows / summaries required.

• channels (list[str]) – channel names, length K.

• pairs (list[tuple[int, int]] | None) – interaction pairs to model (default: all upper-triangular pairs).

• pair_signs (dict[tuple[int, int], str] | None) – sign-informed prior family per pair (default: all "weak"). Keys are normalized to (min, max) orientation; conflicting duplicate entries raise.

• likelihood (str) – panel observation family — "normal" (default, byte-identical to the original graph), "studentt" for heavy-tailed robust continuous KPIs (a learned tail df nu joins sigma; outlier geo-weeks stop dragging the surface), or "negbinomial" for count KPIs (y must be non-negative integer counts; the summary block stays Gaussian either way). Link caveat for the planner: with "negbinomial" the marginal readouts (marginal_roas, expected_regret, free-mode allocation) differentiate the LATENT surface R, while the observable count mean is softplus(mu) whose derivative is sigmoid(mu) — so the count-scale marginal response is sigmoid(mu) * dR/ds. Since sigmoid(mu) ≈ 1 whenever mu >> 1 (typical counts), the readouts are asymptotically exact; near mu ≈ 0 (per-row means of a few counts) they OVERSTATE the count-scale marginal response by up to 1/sigmoid(mu). fit() warns loudly when the observed count mean is below ~20. Fixed-budget Thompson argmax is unaffected (softplus is monotone).

• time_effect (str) – "none" (default) or "national" — a zero-centered hierarchical per-period national shock tau_t added to the panel mean (a non-empty panel requires data["period_idx"]; see model()). A summaries-only fit needs NO period identity: there are no panel rows for tau_t to index and a national per-period constant cancels in the lift difference, so the summary path is unchanged (no tau sites are sampled).

• gamma_scale (float) – interaction prior scale — the most consequential knob; audit it with prior_sensitivity() (guide §8.2).

• spend_ref (ndarray | None) – per-channel reference constant used to scale spend, carried on the returned Posterior for dollar mapping (see to_dollars()).

• prior_scaling (str) – "auto" (default) derives the intercept/noise/effect prior scales from the evidence so natural-unit KPIs at any magnitude fit sanely; "unit" reproduces the original O(1) priors exactly (see _resolve_prior_scaling()).

• discount_half_life (float | None) – information-decay half-life h in WEEKS (None — the default — is the historical static fit, every row weighing the same forever). When set, each panel row’s likelihood is discounted by its age: the observation scale is inflated to sigma * exp(0.5 * lambda * age) with lambda = ln2 / h — the math page’s Eq. 22 decay clock moved INTO the likelihood (West–Harrison discounting; the count family power-scales the log-likelihood by exp(-lambda * age) instead, since it has no noise scale). Row ages come from data["row_age"] ((N,) weeks, 0 = newest — LearningState maintains this automatically) or, failing that, are derived from data["period_idx"] as max - period; a discounted panel fit with neither raises. Summaries may carry an "age_weeks" key — their SE is inflated by the same clock. The effect: effective sample size SATURATES instead of growing without bound, so the posterior keeps an honest variance floor when media behaviour drifts (a static fit’s bands shrink like 1/sqrt(rows) while the truth moves — narrow and wrong). Use estimate_half_life() to measure h from the programme’s own wave-to-wave drift; a measured lambda ~ 0 recovers full pooling.

• spline_prior (str) – monotone-spline weight prior — "iid" (default, the historical independent LogNormal(0, 1) weights, byte-identical) or "pspline": a first-order random walk on the log-weights with a learned per-channel smoothness w_tau (P-spline / adaptive-smoothing construction). tau -> 0 collapses to the neutral near-linear ramp, so the EFFECTIVE basis dimension grows only as the designed cells earn it — the adaptive answer to the fixed-basis approximation-bias floor (math page Eq. 24) without hand-tuning the knot count. Only valid with activation="monotone_spline"; the public w1..wJ sites are unchanged (the acquisition layer and planner need no changes).

Return type:

Posterior

Returns:

A Posterior with merged-chain numpy samples, R-hat/ESS diagnostics on the key parameters, and diagnostics["evidence"] = {"n_rows", "n_summaries"}.

class mmm_framework.continuous_learning.Posterior(samples, channels, pairs, pair_signs=<factory>, activation='hill', likelihood='normal', time_effect='none', spend_ref=None, diagnostics=<factory>)

Bases: object

Posterior draws plus the metadata the planner needs.

samples holds plain numpy arrays keyed by site name (beta, kappa, alpha, A, sigma_a, a_geo, sigma, and one gamma_<ci>_<cj> per pair; the intercept/noise sites are absent for a summaries-only fit). A likelihood="negbinomial" fit carries phi (the NB concentration) instead of sigma; a likelihood="studentt" fit carries nu (the tail df) alongside sigma; a time_effect="national" fit adds sigma_tau and tau ((draws, n_period)). spend_ref is the per-channel reference constant used to scale spend — convert with to_dollars() / to_scaled().

samples: dict[str, ndarray]

channels: list[str]

pairs: list[tuple[int, int]]

pair_signs: dict[tuple[int, int], str]

activation: str = 'hill'

likelihood: str = 'normal'

time_effect: str = 'none'

spend_ref: ndarray | None = None

diagnostics: dict[str, Any]

property n_channels: int

property n_draws: int

gamma_matrix(d)

Assemble the (K, K) upper-triangular gamma matrix for draw d.

Return type:

ndarray

property shape_names: tuple[str, ...]

The activation’s shape-parameter site names, in order.

draw_params(d)

Per-draw params for the surface functions, activation-agnostic.

Returns {beta, gamma, shape, act_fn, activation} where shape is a tuple of (K,) arrays in the activation’s parameter order and act_fn is its JAX activation. (For the Hill default the shape is (kappa, alpha).)

Return type:

dict[str, Any]

gamma_summary()

Per-pair posterior mean and 5/95 percentiles of the synergy.

Return type:

dict[str, dict[str, float]]

__init__(samples, channels, pairs, pair_signs=<factory>, activation='hill', likelihood='normal', time_effect='none', spend_ref=None, diagnostics=<factory>)

mmm_framework.continuous_learning.default_pairs(n_channels)

All upper-triangular channel pairs (i, j) with i < j.

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.pair_name(channels, pair)

Posterior site name for a pair’s interaction, gamma_<ci>_<cj>.

Return type:

str

mmm_framework.continuous_learning.demote_channel(channels, name, *, pairs=None, base=None)

Demote a non-randomizable channel (a walled garden) to main-effect-only.

Every interaction that touches name becomes "zero" (prior-dominated); its main effect is still identified. Flag gamma(., name) as the least-trustworthy parameter in any decision (guide §5.4).

Return type:

dict[tuple[int, int], str]

mmm_framework.continuous_learning.probe_pairs_excluding(channels, name, *, pairs=None)

Drop the off-axis design cells for a demoted channel’s pairs.

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.refit_fn_from_data(base_data, *, channels, pairs=None, pair_signs=None, activation='hill', likelihood='normal', beta_scale=1.0, gamma_scale=0.8, num_warmup=200, num_samples=200, num_chains=1, seed=7, prior_scaling='auto', time_effect='none')

Build a refit_fn(extra_spend, extra_geo_idx, extra_y) -> Posterior.

The knowledge-gradient acquisition (planner.knowledge_gradient()) fantasizes wave outcomes and refits; this closure appends the fantasy rows to base_data and runs a short NUTS chain. This is the expensive path — see guide §9.1 for the Laplace-update replacement in production.

Requires a panel base dataset: fantasy observations are geo-week rows, so a summaries-only data dict has nothing to append them to.

likelihood must match the posterior the knowledge gradient is scoring: the fantasy y values are drawn from that same observation family (Gaussian, Student-t, or NB counts), so refitting them under a different family would score the design against a model nobody will run.

time_effect other than "none" raises NotImplementedError: the knowledge-gradient fantasy rows carry no period identity, so a national tau_t cannot be refit on the augmented panel. A period_idx present in base_data is DROPPED from the refit base for the same reason (the refit models no time effect, so the index is inert either way).

Return type:

Callable[[ndarray, ndarray, ndarray], Posterior]

mmm_framework.continuous_learning.central_composite(center, delta, probe_pairs)

Build the CCD cells (rows = cells, columns = channels).

Parameters:

• center (ndarray) – the operating allocation, shape (K,) (scaled spend).

• delta (float) – multiplicative trust-region step in (0, 1] (e.g. 0.6).

• probe_pairs (list[tuple[int, int]]) – channel pairs to identify the cross-partial for.

Return type:

ndarray

Returns:

A (n_cells, K) array of non-negative scaled allocations, with n_cells = 1 + 2K + 2 * len(probe_pairs) + K.

mmm_framework.continuous_learning.assign_geos(design, n_geo, rng, *, n_holdout=0, center=None, baseline=None)

Assign CCD cells to geos — shuffled round-robin, or stratified/blocked.

Without a baseline (the default, byte-identical to the historical behavior) cells are tiled round-robin and shuffled: balanced cell counts, zero covariate awareness, and the first n_holdout geos (positional) become holdouts.

With a baseline (a per-geo covariate such as the pre-period KPI level, matched-market style) the assignment is a classic blocked randomization: geos are sorted by baseline, walked in blocks of n_cells, and each block receives a random permutation of the cell indices (the ragged tail gets a random subset without replacement). Every cell’s geos are then spread evenly across the baseline distribution, so between-cell baseline means are nearly equal. Holdouts are carved FIRST and are stratum-aware too: exactly n_holdout evenly spaced positions in baseline-sorted order, so the status-quo counterfactual spans the baseline range AND honors the requested count (a strided pick would silently under-deliver whenever n_holdout does not divide n_geo evenly).

Parameters:

• design (ndarray) – CCD cells, shape (n_cells, K).

• n_geo (int) – number of geos.

• rng (Generator) – a numpy random generator (caller owns the seed).

• n_holdout (int) – hold n_holdout geos at center for the test window (a status-quo counterfactual). Requires center.

• center (ndarray | None) – the status-quo allocation for holdout geos, shape (K,).

• baseline (Union[ndarray, Sequence[float], None]) – optional per-geo covariate, length n_geo, positionally aligned with the geo indices (the caller resolves geo ids). None keeps the legacy shuffled round-robin path.

Return type:

tuple[ndarray, ndarray]

Returns:

(geo_alloc, cell_idx) where geo_alloc is (n_geo, K) (each geo’s test allocation) and cell_idx is (n_geo,) (the design-row index, or -1 for a holdout geo).

class mmm_framework.continuous_learning.TrueWorld(beta, gamma_pairs, channels, kappa=None, alpha=None, activation='hill', shape=<factory>, pairs=<factory>, a_level=4.0, sigma_a=1.0, phi_true=10.0, nu_true=4.0)

Bases: object

A known response surface for recovery / closure testing.

gamma_pairs is aligned to pairs (one synergy value per pair); the matrix form is assembled on demand. A/sigma_a set the geo-intercept distribution. The activation is pluggable: pass Hill kappa/alpha (the default and back-compatible path), or set activation + shape for another family (e.g. activation="logistic", shape={"lam": …}). phi_true is the NegativeBinomial concentration used when a simulation asks for noise_family="negbinomial"; nu_true is the Student-t tail df used by noise_family="studentt" (each ignored otherwise).

beta: ndarray

gamma_pairs: ndarray

channels: list[str]

kappa: ndarray | None = None

alpha: ndarray | None = None

activation: str = 'hill'

shape: dict[str, ndarray]

pairs: list[tuple[int, int]]

a_level: float = 4.0

sigma_a: float = 1.0

phi_true: float = 10.0

nu_true: float = 4.0

property n_channels: int

act_fn()

The JAX activation for this world’s family.

shape_tuple()

Shape-parameter arrays in the activation’s parameter order.

Return type:

tuple

gamma_matrix()

Return type:

ndarray

response_mean(spend_matrix)

Mean incremental response for a (N, K) spend panel (no noise).

Return type:

ndarray

answer_key()

Return type:

dict[str, object]

__init__(beta, gamma_pairs, channels, kappa=None, alpha=None, activation='hill', shape=<factory>, pairs=<factory>, a_level=4.0, sigma_a=1.0, phi_true=10.0, nu_true=4.0)

mmm_framework.continuous_learning.make_world(seed=0, channels=None)

A reproducible, non-trivial 4-channel world for tests and the demo.

Channels ["Chatter", "Pulse", "Orbit", "Vibe"] (generic social-network surfaces) with a mix of strong and weak main effects, distinct half- saturations and Hill shapes, and a synergy structure that includes a genuine cannibalization (negative gamma) and a complementarity (positive gamma) so recovery has to get signs right.

Return type:

TrueWorld

mmm_framework.continuous_learning.make_world_logistic(seed=0, channels=None)

A known logistic (exponential-saturation) world, sibling of make_world().

Same channels and synergy structure, but each channel saturates as f(s) = 1 - exp(-lam * s) — a concave curve with no S-shape — rather than a Hill. Demonstrates that the whole loop (fit, plan, stop, animate) is activation-agnostic. The saturation rates lam are chosen so half- saturation ln(2)/lam lands in the same O(1) scaled range as the Hill kappa above, keeping the two worlds comparable in scale.

Return type:

TrueWorld

mmm_framework.continuous_learning.make_world_hill_mixture(seed=0, channels=None)

A known weighted-sum-of-two-Hills world (a misspecification stress test).

Each channel’s true response is w·Hill(κ1,α1) + (1-w)·Hill(κ2,α2) with a low-κ, high-α early component (a soft activation threshold) plus a high-κ later component — a two-phase / shoulder shape a single Hill can only average over and a logistic (concave, no inflection) cannot represent at all. Fit it with activation="hill" (mild misspecification) or "logistic" (severe) to see what a wrong response family does to the recovered curve, the funding line, and the recommended allocation.

Return type:

TrueWorld

mmm_framework.continuous_learning.simulate_panel(world, center, *, n_geo=80, t_pre=6, t_test=10, delta=0.6, probe_pairs=None, noise=0.6, noise_family='normal', tau_scale=0.0, n_holdout=0, a_geo=None, adstock_alpha=None, adstock_l_max=8, stratify=False, seed=0)

Simulate a CCD wave against a known world.

Returns the data contract (spend, geo_idx, n_geo, y, period_idx) plus design, geo_alloc, cell_idx, a_geo, tau_true and answer_key. Pass an existing a_geo (and the same seed offset) to simulate a later wave over the same geos with their baselines intact.

period_idx is always provided (np.repeat(arange(t_pre + t_test), n_geo) — the row order is week-major, pre block then test block); a fit only uses it when it opts into time_effect="national".

noise_family="negbinomial" draws count outcomes (gamma–Poisson with mean softplus(mu) and concentration world.phi_true); the default "normal" path is byte-identical to the old Gaussian draw. tau_scale > 0 adds true national per-period shocks tau_t ~ Normal(0, tau_scale) to mu (0.0 leaves y byte-identical and draws nothing).

adstock_alpha adds carryover: the response is driven by the geometric-adstocked spend series (within each geo over weeks), while the returned spend stays raw. Fitting on the raw panel is then biased; the preprocess.adstock_prepass recovers it (guide §9.4).

stratify=True (opt-in — the default keeps the historical rng stream byte-identical) resolves a_geo BEFORE the geo assignment and blocks the randomization on it (assign_geos(..., baseline=a_geo)), matching the production practice of stratifying on the pre-period KPI level.

Return type:

dict[str, object]

mmm_framework.continuous_learning.simulate_wave(world, design, a_geo, *, t_test=10, center=None, n_holdout=0, noise=0.6, noise_family='normal', tau_scale=0.0, stratify=False, seed=0)

Simulate a single test window for a recentered design over the SAME geos.

Used by the closed loop to generate wave t > 0 outcomes: the geos and their baselines (a_geo) persist, only the test allocations change. period_idx (wave-LOCAL, np.repeat(arange(t_test), n_geo)) is always provided; ingest() applies the cross-wave offset. noise_family / tau_scale behave as in simulate_panel() (defaults byte-identical to the old draw). stratify=True (opt-in) blocks the geo→cell randomization on the known per-geo baselines a_geo.

Return type:

dict[str, object]

mmm_framework.continuous_learning.drift_world(world, *, rate=0.08, drift_beta=True, drift_shape=True, drift_gamma=False, seed=0)

A slightly different TrueWorld — media behaviour drifting.

Applies one step of multiplicative log-normal jitter, x -> x * exp(Normal(0, rate)), to the channel ceilings beta and the activation shape params (bounded params are clipped back to their family’s support; drift_gamma=True also jitters the synergy magnitudes, sign-preserving). Applying it once per wave makes the response surface a geometric random walk — the changing-media- behaviour regime the information discount (fit(discount_half_life=...)) is built for. Channels, pairs and the geo-intercept hypers are untouched, so pair it with simulate_wave() over a fixed a_geo: the baselines stay pinned while the RESPONSE moves.

Return type:

TrueWorld

mmm_framework.continuous_learning.allocate_under_sample(params, B, value, *, x0=None, n_starts=3, seed=0, mode='fixed', cap=None, group_budgets=None)

Optimal allocation + profit under one posterior draw’s params.

group_budgets fixes sub-budgets over channel groups (arms): [(indices, B_g), ...] adds sum(s[indices]) == B_g per group.

Return type:

tuple[ndarray, float]

mmm_framework.continuous_learning.thompson_wave(post, B, value, *, q=300, seed=0, mode='fixed', cap=None, n_starts=3, group_budgets=None)

A posterior over the optimal split — solve the allocation for q draws.

Returns (allocs, profits, draws): allocs is (q, K) (the mean is the recommended allocation, the spread is the exploration signal), profits is (q,) (each draw’s optimum under itself), and draws are the sampled draw indices.

Return type:

tuple[ndarray, ndarray, ndarray]

mmm_framework.continuous_learning.recommend_allocation(post, B, value, **kwargs)

The recommended allocation — the Thompson posterior mean split.

Return type:

ndarray

mmm_framework.continuous_learning.marginal_roas(post, alloc, value, *, q=300, seed=1)

Posterior of value * dR/ds_c at alloc — the funding line.

Returns (mean_mroas, prob_above_line, draws_mroas). A channel is funded where prob_above_line > 0.5.

Return type:

tuple[ndarray, ndarray, ndarray]

mmm_framework.continuous_learning.expected_regret(post, B, value, *, q=300, seed=2, mode='fixed', cap=None, n_starts=3, group_budgets=None)

Expected regret of acting on the consensus allocation (guide §7.4).

regret_d = profit(best-for-draw-d under d) - profit(consensus under d) >= 0 — the profit still on the table from posterior uncertainty. Each draw’s “best” is max(cold multistart, a solve warm-started FROM the consensus, profit at the consensus) (review fix F3): the warm-started second pass catches per-draw optima the cold multistart missed near the consensus, so the regret is less biased low and the ENBS stop fires less prematurely. Returns (E[regret], consensus_alloc, alloc_sd, optimal_profit_sd).

Return type:

tuple[float, ndarray, ndarray, float]

mmm_framework.continuous_learning.plan_from_posterior(post, B, value, *, q=300, seed=0, mode='fixed', cap=None, n_starts=3, group_budgets=None)

One Thompson pass producing every decision readout coherently (fix F4).

Runs thompson_wave() once; the recommendation is the draw-mean (rescaled to the simplex in fixed mode); the funding line and the warm-started regret pass reuse the SAME draw indices, so the funded set, the consensus, and E[regret] all describe one allocation.

Return type:

PlanResult

class mmm_framework.continuous_learning.PlanResult(channels, B, value, mode, recommendation, consensus, allocs, profits, draws, mroas_mean, prob_above_line, mroas_draws, e_regret, alloc_sd, profit_sd)

Bases: object

Every per-wave decision readout, derived from ONE Thompson sample.

Historically recommend/funding/regret each re-sampled with a different seed (planner seeds 0/1/2), so the funded set and the reported regret referred to different consensus vectors and paid 2-3x the SLSQP cost. Here all readouts share the same q draws: recommendation is the Thompson mean (rescaled onto the budget simplex in fixed mode), consensus is that same vector, the funding line is evaluated at it, and the regret pass is warm-started from it.

channels: list[str]

B: float

value: float

mode: str

recommendation: ndarray

consensus: ndarray

allocs: ndarray

profits: ndarray

draws: ndarray

mroas_mean: ndarray

prob_above_line: ndarray

mroas_draws: ndarray

e_regret: float

alloc_sd: ndarray

profit_sd: float

to_dict()

A small JSON-safe snapshot (summary stats, not the draw matrices).

Return type:

dict[str, Any]

__init__(channels, B, value, mode, recommendation, consensus, allocs, profits, draws, mroas_mean, prob_above_line, mroas_draws, e_regret, alloc_sd, profit_sd)

mmm_framework.continuous_learning.posterior_optimal_allocation(post, B, value, *, q=200, seed=3, mode='fixed', cap=None, n_starts=4, group_budgets=None)

argmax_a E_post[profit(a)] and its value — optimize the mean surface.

Return type:

tuple[ndarray, float]

mmm_framework.continuous_learning.knowledge_gradient(post, candidate_design, refit_fn, B, value, *, n_fantasy=10, t_test=10, n_geo=None, noise=None, mode='fixed', cap=None, q=120, seed=4)

One-step-lookahead EVSI for a candidate test design (guide §7.5).

Return type:

float

``KG(d) = E_y[ max_a profit(a | posterior updated with fantasised y) ]

• max_a profit(a | current posterior)``

For each fantasy: draw params, simulate the candidate design’s outcomes, refit via refit_fn(extra_spend, extra_geo_idx, extra_y) -> Posterior, and re-optimize. Expensive — refit_fn runs a (short) NUTS chain per fantasy; in production swap it for a Laplace/conjugate update (guide §9.1).

noise=None (default) fantasizes with the posterior mean of sigma — the fitted observation noise — rather than a hard-coded constant. Fantasy outcomes are drawn from the posterior’s own observation family (_fantasy_outcomes()): Gaussian, Student-t (posterior-mean nu), or NB counts (posterior-mean phi; noise ignored). refit_fn must refit under that same family — build it with model.refit_fn_from_data()’s matching likelihood.

Requires a panel-fitted posterior: fantasies simulate geo-week outcomes off the geo intercepts (a_geo or the A/sigma_a hypers), which a summaries-only fit never samples.

mmm_framework.continuous_learning.response_grid(post, spend_matrix, draws)

Incremental response for every (spend row, draw) pair -> (G, D).

Activation-agnostic (reads post.activation): the visualizations use this to build the posterior mean/uncertainty surfaces for any activation family, so the same animation code renders a Hill world or a logistic one.

Return type:

ndarray

mmm_framework.continuous_learning.enbs(e_regret, *, margin, population, wave_cost)

Expected net benefit of sampling = E[regret] * margin * population - cost.

E[regret] is profit-per-affected-unit on the response-curve scale; margin and population convert it to a total dollar value of resolving the uncertainty, which a wave must beat to be worth running.

Return type:

float

mmm_framework.continuous_learning.should_stop(e_regret, *, margin, population, wave_cost)

(stop, enbs) — stop when no wave’s expected value clears its cost.

Return type:

tuple[bool, float]

mmm_framework.continuous_learning.adstock_panel(spend, n_geo, t_pre, t_test, *, alpha, l_max=8, normalize=True)

Geometric-adstock each channel’s spend series within each geo.

Reuses transforms.adstock (normalize=True keeps the kernel a weighted moving average so total magnitude stays in the coefficient). Returns a spend array of the same shape, in the same row order.

Return type:

ndarray

mmm_framework.continuous_learning.adstock_prepass(data, t_pre, t_test, *, alpha, l_max=8)

Return a copy of the data contract with the spend adstocked (guide §9.4).

Fit the surface on this instead of the raw panel when the response has carryover; the decay alpha is a hyperparameter (sweep it, or set it from a known channel half-life).

Return type:

dict[str, Any]

mmm_framework.continuous_learning.cuped_covariate(data, t_pre)

Per-geo, mean-centered pre-period KPI — the CUPED control covariate.

Return type:

ndarray

mmm_framework.continuous_learning.cuped_adjust(data, t_pre)

CUPED-adjust the outcome and report the variance reduction.

y_adj = y - theta * x_pre[geo] with theta = Cov(y_test, x_pre) / Var(x_pre) (the regression of the geo’s test-period mean on its pre-period covariate). Returns (adjusted_data, info) where info carries theta, rho (the correlation CUPED exploits) and var_reduction (1 - rho^2, the fraction of geo-level outcome variance removed).

Incompatible with likelihood="negbinomial": the adjustment mutates y to non-integer (and possibly negative) values, which the count likelihood’s validation rejects. For count KPIs fit the raw counts — the geo intercept plus the national time effect (time_effect="national") covers most of what CUPED buys here.

Return type:

tuple[dict[str, Any], dict[str, float]]

mmm_framework.continuous_learning.laplace_knowledge_gradient(post, design_cells, B, value, *, sigma=None, n_geo=80, t_test=10, n_outcomes=64, mode='fixed', cap=None, n_starts=3, cell_weights=None, seed=0)

Decision-aware EVSI of a candidate design, with no refit (guide §9.1).

Gaussian-linear surrogate for mmm_framework.continuous_learning.planner.knowledge_gradient(): the pre-posterior spread of the updated mean is V = Sigma - Sigma_post; fantasy means eta_m ~ N(mu, V) are mapped back to valid surface parameters, re-optimized and averaged against the current best. Orders designs the same way as the NUTS KG at a fraction of the cost.

Works for any registered activation and fitted observation family (the Fisher weights come from observation_unit_info()); sigma=None derives the Gaussian-family noise scale from the posterior.

Each fantasy’s re-optimization is warm-started at the base (mu) allocation, in addition to the usual uniform/random multi-starts. Fantasy means are draws local to mu (V is a reduction in spread from the prior), so the true optimum is almost always near the base one; without the anchor, a non-concave surface (gamma can be negative) evaluated with only a handful of random restarts can converge to a different local optimum depending on the platform’s BLAS/LAPACK build (SLSQP’s floating- point path through the same starting point is not bitwise-portable), occasionally flipping a close two-candidate comparison. Anchoring removes that basin-hopping risk at no extra cost (it replaces one of the existing random starts, per _starts()).

Return type:

float

mmm_framework.continuous_learning.design_eig(post, design_cells, *, sigma=None, target='all', n_geo=80, t_test=10, cell_weights=None)

Pure information gain of a design (D-optimal target="all" or D_s-optimal target="gamma" over the synergy sub-block).

Works for any registered activation and any fitted observation family; sigma=None derives the noise scale (Gaussian families) from the posterior. sigma is ignored for a count posterior — the GLM weights come from phi and the baseline instead.

Return type:

float

mmm_framework.continuous_learning.design_information(design_cells, theta_bar, *, sigma=None, k=None, pairs=None, tmap=None, unit_info=None, cell_weights=None, n_geo=80, t_test=10, residualize=True)

Fisher information Lambda (dim×dim) a design carries about eta.

Lambda = sum_c w_c u_c (g_c - g_bar)(g_c - g_bar)^T where g_c is the unconstrained-space parameter gradient at cell c, u_c the observation family’s per-cell unit information (unit_info; defaults to the homoskedastic Gaussian 1/sigma^2), and g_bar the weighted mean across cells (profiling the geo intercept). residualize=False skips the centering (treats the baseline as known).

theta_bar is the unconstrained eta vector from theta_moments(). Pass the matching tmap; when omitted, a default Hill map with identity gamma transforms is built from (k, pairs) — only correct for a Hill posterior with default (“weak”/”zero”) pair signs.

Return type:

ndarray

mmm_framework.continuous_learning.gaussian_eig(sigma0, lam, *, idx=None)

EIG (nats) of a parameter block under a Gaussian-linear update.

0.5 * (logdet Sigma0_SS - logdet Sigma_post_SS) for the block idx (all parameters if None) — the entropy the design removes. Both log-dets run through the _logdet_psd() symmetrize-and-clip guard: Sigma_post is PSD in exact arithmetic (Lambda >= 0) but the finite-precision double inversion can leave a marginal sub-block — equivalently the Schur complement Lambda_{S|rest} — with tiny negative eigenvalues.

Return type:

float

mmm_framework.continuous_learning.observation_unit_info(post, design_cells, tmap, mu_eta, *, sigma=None)

Per-cell Fisher information of ONE observation, for the fitted family. :rtype: ndarray

• "normal": 1 / sigma^2 (homoskedastic — constant across cells).

• "studentt": (nu + 1) / ((nu + 3) sigma^2) with the posterior-mean tail df — the classic heavy-tail efficiency discount (< the Gaussian info at any finite nu, recovering it as nu -> inf).

• "negbinomial": the GLM weight through the model’s softplus link, sigmoid(eta_c)^2 / (m_c + m_c^2 / phi) with eta_c = baseline + R(s_c; mu), m_c = softplus(eta_c) and the posterior-mean phi. The baseline is the posterior-mean A (or the grand mean of a_geo when the hyper is absent). sigma is ignored — a count family has no Gaussian noise scale.

sigma=None reads the posterior-mean sigma site for the Gaussian families; a missing required site raises ValueError (summaries-only fits never sample the observation sites, and no fixed guess is meaningful across KPI scales).

mmm_framework.continuous_learning.surrogate_validity(post, *, tmap=None, tail_prob=0.25, khat_threshold=0.7, skew_threshold=0.5, kurt_threshold=1.0)

Diagnose whether N(mu, Sigma) is a valid stand-in for the posterior.

The Laplace acquisitions replace the carried NUTS posterior with its Gaussian moment-match in the unconstrained eta space (Eq. 12 of the math page). This runs entirely off the existing draws — no densities, no refit — and checks the two ways that stand-in fails: :rtype: dict[str, Any]

• Shape: per-parameter skewness and excess kurtosis of the unconstrained draws (a Gaussian has both ~0; the transforms already absorb the generic positivity skew, so what remains is real).

• Tails / ridges: the generalized-Pareto shape khat of the Mahalanobis-distance exceedances above their 1 - tail_prob quantile. Under a genuinely Gaussian posterior the squared distances are chi-square (exponential-class tail, khat ~ 0); a skewed, heavy-tailed or ridge-shaped posterior (the negative-gamma cannibalisation ridge is the expected offender) inflates khat. The alarm bar is the same 0.7 PSIS uses (Vehtari et al. 2024), estimated with the same Zhang–Stephens fit.

Returns a JSON-safe dict: ok (False when khat clears the threshold, any parameter is flagged, or there are too few draws to tell), khat, params_flagged, max_abs_skew, max_excess_kurtosis, and the per_param detail. When ok is False, spot-check or replace the wave’s Laplace scores with the NUTS-refit knowledge_gradient() (the correction path in the literature is a skewness-aware expansion — Rue, Martino & Chopin 2009).

mmm_framework.continuous_learning.theta_map(post)

Build the ThetaMap for a posterior’s activation + pair signs.

Return type:

ThetaMap

mmm_framework.continuous_learning.theta_moments(post, *, ridge=1e-06, tmap=None)

Gaussian (mu, Sigma) matched to the posterior — in eta space.

The match happens in the unconstrained reparameterization of theta_map(): positive parameters in log space, bounded ones through a scaled logit, sign-constrained synergies through a sign-aware log. Any activation with a SHAPE_TRANSFORMS entry is supported; an unknown family raises NotImplementedError (the planner’s decision readouts stay activation-agnostic).

Return type:

tuple[ndarray, ndarray]

mmm_framework.continuous_learning.bocd(x, *, hazard=0.04, mu0=None, kappa0=1.0, alpha0=3.0, beta0=None, threshold=0.5)

Bayesian online changepoint detection (Adams & MacKay 2007).

Normal observations with unknown mean and variance: the conjugate Normal–Inverse-Gamma model, whose posterior predictive is a Student-t — so the detector needs no known noise scale. hazard is the constant per-step changepoint prior P(r_t = 0) (1/hazard = expected run length; the default expects a run of ~25 periods, conservative for a 10-week wave). mu0/beta0 default empirically: the series median, and a noise variance from the MAD of the first differences (robust to the very level shift being tested for — a full-series spread estimate would fold the shift into the noise prior and mask it, while differences see only the step at the break). An empirical-Bayes convenience that is appropriate for the retrospective wave check; pass explicit priors for genuinely online use.

Returns a BocdResult; changepoints flags every step t >= 2 with P(a sustained restart happened one step ago) > threshold — see BocdResult for why the flag is P(r_t = 2), not P(r_t = 0). A flag at t means the break onset was at t - 1.

Return type:

BocdResult

class mmm_framework.continuous_learning.BocdResult(cp_prob, run_length_map, changepoints, threshold)

Bases: object

Run-length posterior summary from bocd().

cp_prob[t] is the posterior probability that a sustained restart happened one step ago — P(r_t = 2 | x_{1:t}), the hypothesis that the current regime contains exactly the last two observations. Two naive alternatives fail: P(r_t = 0) is identically the hazard under a constant-hazard BOCD (the changepoint and growth branches share the same predictive weights, so their mass split is H : 1-H regardless of the data), and P(r_t = 1) fires on any single-point outlier (the fresh run’s wide prior predictive absorbs a spike that the tight long-run model rejects). Requiring the restarted run to explain two consecutive points is the minimum evidence for a regime rather than a blip — a spike reverts at the next step and the long run recovers the mass. The price is a one-period detection lag: a flag at t means the break onset was at t - 1. run_length_map[t] is the MAP run length after step t; changepoints are the steps (t >= 2) whose cp_prob cleared the threshold.

cp_prob: ndarray

run_length_map: ndarray

changepoints: list[int]

threshold: float

__init__(cp_prob, run_length_map, changepoints, threshold)

mmm_framework.continuous_learning.wave_stationarity_check(wave_data, *, hazard=0.04, threshold=0.5, min_geos_per_cell=1)

Detect a mid-wave regime change in a designed wave’s contrasts.

For every design cell, build the per-period treatment-minus-control series (mean outcome of the cell’s geos minus the control’s — holdout geos cell_idx == -1 when the wave has them, else the center cell, design row 0) and run bocd() on it. Common shocks cancel in the contrast; a shock that hits treatment differently from control — the kind that biases the readout — appears as a level shift and fires the flag.

wave_data needs the simulate/ingest contract keys spend, geo_idx, y, period_idx and the per-geo cell_idx. Raises ValueError when a required key is missing (the guard cannot run blind). Returns a StationarityReport; when flagged, either censor_periods() the affected periods before ingesting or extend / repeat the wave.

Return type:

StationarityReport

class mmm_framework.continuous_learning.StationarityReport(flagged, break_periods, cell_breaks, control, n_test_periods, n_series, threshold, max_cp_prob, note='', cp_prob=<factory>)

Bases: object

Result of wave_stationarity_check() (JSON-safe via to_dict()).

flagged: bool

break_periods: list[int]

cell_breaks: dict[int, list[int]]

control: str

n_test_periods: int

n_series: int

threshold: float

max_cp_prob: float

note: str = ''

cp_prob: dict[int, list[float]]

to_dict()

Return type:

dict[str, Any]

__init__(flagged, break_periods, cell_breaks, control, n_test_periods, n_series, threshold, max_cp_prob, note='', cp_prob=<factory>)

mmm_framework.continuous_learning.censor_periods(wave_data, periods)

Drop every row of the given (global) periods from a wave dict.

Returns a shallow copy with spend/geo_idx/y/period_idx filtered row-wise (and any other key left untouched) — feed the result to LearningState.ingest() in place of the contaminated wave. Censoring from a detected break onward is the conservative choice when the regime did not revert (a level shift, not a spike); censoring just the flagged periods suits a transient.

Return type:

dict[str, Any]

mmm_framework.continuous_learning.to_scaled(dollars, spend_ref)

Convert dollar spend to scaled units: dollars / spend_ref.

Parameters:

• dollars (ndarray) – spend in dollars — a (K,) allocation or an (N, K) panel (any leading shape; the channel axis is the last one).

• spend_ref (ndarray) – dollars per scaled unit, shape (K,), strictly positive.

Return type:

ndarray

Returns:

The scaled spend, same shape as dollars (float64).

mmm_framework.continuous_learning.to_dollars(scaled, spend_ref)

Convert scaled spend back to dollars: scaled * spend_ref.

Parameters:

• scaled (ndarray) – scaled spend — a (K,) allocation or an (N, K) panel.

• spend_ref (ndarray) – dollars per scaled unit, shape (K,), strictly positive.

Return type:

ndarray

Returns:

The dollar spend, same shape as scaled (float64).

mmm_framework.continuous_learning.experiments_to_summaries(experiments, *, channels, spend_ref, center_scaled, period_days=7.0)

Convert registry experiment dicts into summary observations.

Parameters:

• experiments (list[dict[str, Any]]) – rows from sessions.list_experiments (plain dicts with status/channel/estimand/value/se plus the design and readout JSON snapshots).

• channels (list[str]) – the program’s channel (or flattened arm) names, length K.

• spend_ref (ndarray) – dollars per scaled unit per channel, shape (K,).

• center_scaled (ndarray) – the program’s baseline allocation in scaled units, shape (K,) — the counterfactual every readout is measured against.

• period_days (float) – registry-date cadence (7 = weekly).

Return type:

tuple[list[dict[str, Any]], list[dict[str, Any]]]

Returns:

(summaries, skipped) — summaries are ready for fit(data={"summaries": ...}) (each also carries experiment_id / channel provenance); skipped items are {"id", "reason"}.

mmm_framework.continuous_learning.posterior_to_payload(post)

A JSON-safe dict capturing a Posterior completely.

Sample arrays become (nested) lists of float64 — Python’s JSON float representation round-trips doubles exactly, so posterior_from_payload() reconstructs bit-identical samples.

Return type:

dict[str, Any]

mmm_framework.continuous_learning.posterior_from_payload(d)

Inverse of posterior_to_payload().

Return type:

Posterior

mmm_framework.continuous_learning.state_to_npz(state, path)

Persist a LearningState to one compressed .npz file.

Arrays (panel, posterior samples, summary spend vectors) are stored natively (exact float64); config/history/pair_signs/geo_ids ride an embedded JSON string under the meta key. Returns the written path.

Return type:

str

mmm_framework.continuous_learning.state_from_npz(path)

Reconstruct a LearningState written by state_to_npz().

Return type:

LearningState

class mmm_framework.continuous_learning.ArmSpec(channels, parents, groups)

Bases: object

The flattened arm layout of a (possibly partially) split channel list.

channels are the flattened arm names (surface dimensions, in order); parents[i] is arm i’s parent channel (== the name itself when the channel is unsplit); groups maps every parent to its arm indices.

channels: list[str]

parents: list[str]

groups: dict[str, list[int]]

property n_arms: int

split_parents()

Parents that were actually split into more than one arm.

Return type:

list[str]

__init__(channels, parents, groups)

mmm_framework.continuous_learning.arm_shares(post, spec, parent, spend_ref, *, breakout_name_map, mode='zero_out', draws=500, rng=None)

Export a parent’s within-channel share composition from a CL posterior.

The outbound bridge to the breakout-weighted MMM’s share calibration (ShareMeasurement in mmm_framework.calibration.likelihood): per posterior draw, compute each of the parent’s arms’ incremental response at the reference spend spend_ref, normalize to a within-parent share simplex, and summarize the draws as shares + the empirical covariance of the K-1 additive log-ratios (ALR, last arm as reference).

Location/covariance consistency: the exported shares are the inverse-ALR (softmax) of the mean log-ratios over the surviving draws – NOT the arithmetic mean of the share draws – so the consumer’s MvNormal location and its observed z_hat = ALR(shares) are exactly mean(z), on the same draws as the exported covariance.

Non-positive responses ("zero_out" mode): a draw where any arm’s zero-out response is non-positive (strong cannibalization at spend_ref) has NO well-defined share composition, so such draws are excluded from both the shares and the ALR covariance rather than floored into the statistics (a floored draw would inject log(eps) outliers that inflate the covariance by orders of magnitude and silently down-weight the share evidence). Any exclusion warns; more than 20% excluded raises (the zero-out shares are ill-defined – use mode="main_effect" or a different spend_ref); fewer than 10 surviving draws raises (too few to estimate the covariance).

Estimand caveat (read before use): the MMM’s breakout_share_<C> is an effectiveness share through one shared response curve at the panel’s spend mix, while an arm share here is a ratio of PER-ARM response curves at spend_ref (with interactions in "zero_out" mode). The two coincide only near the panel’s observed sub-stream spend mix and under the breakout model’s own shared-curve assumption – so spend_ref should approximate the panel’s operating point, NOT a reallocated/optimal spend.

Parameters:

• post (Posterior) – A fitted continuous-learning Posterior over the arms surface (each arm a full dimension).

• spec (ArmSpec) – The ArmSpec describing the flattened arm layout.

• parent (str) – The split parent channel to export (must have >= 2 arms).

• spend_ref (ndarray) – Scaled reference spend vector over ALL arms, shape (spec.n_arms,) – the operating point the shares are read at.

• breakout_name_map (dict[str, str]) – REQUIRED mapping from each arm’s sub-name (the part after ARM_SEP, e.g. "Brand" from "Search │ Brand") to the MFF breakout column name (e.g. "Search_Brand"). No fuzzy fallback: a missing key raises. The returned breakouts are the mapped names in the ARM (group) order; the MMM consumer requires them to match its model order exactly, so build the map (and the arm order) against BreakoutWeightedParams.breakout_groups.

• mode (str) – "zero_out" (default) – arm i’s response is R(s_ref) - R(s_ref with arm i zeroed), capturing interactions gamma; "main_effect" – beta_i * act_fn(s_ref)_i, main-effect only (strictly positive since beta is HalfNormal).

• draws (int) – Maximum posterior draws to use; the posterior is subsampled without replacement when it has more.

• rng (Any) – Seed or numpy.random.Generator for the subsample.

Returns:

{channel, breakouts, shares, log_ratio_cov, distribution, source}. channel is set to parent; overwrite it if the MMM’s virtual channel name differs. source records {"mode", "spend_ref", "n_draws", "n_excluded"} as provenance for the double-counting guard (n_draws counts the SURVIVING draws the summary is computed on).

Return type:

A ShareMeasurement-shaped payload

Raises:

ValueError – If more than 20% of the used draws have a non-positive arm response in "zero_out" mode, or if fewer than 10 draws survive the exclusion.

mmm_framework.continuous_learning.expand_arms(channels, arms)

Flatten channels with arms (parent -> sub-names) into an ArmSpec.

Channels absent from arms stay single arms named after themselves; split channels expand in place, preserving the channel order. Raises on an unknown parent, an empty/duplicated sub-list, or a flattened-name clash.

Return type:

ArmSpec

mmm_framework.continuous_learning.within_parent_pairs(spec)

All sibling pairs (i, j), i < j, sharing a parent (substitutes).

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.cross_parent_pairs(spec, pairs_of_parents=None)

Arm pairs spanning two different parents.

pairs_of_parents restricts to specific (unordered) parent pairs — probe only the decision-pivotal cross-channel synergies rather than paying the full off-axis cell cost for every arm combination.

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.default_arm_pair_signs(spec, *, within='neg', base=None)

The default sign map for an arms surface.

Within-parent siblings get within (default "neg" — shared-audience substitution: creatives/keywords cannibalize each other), cross-parent pairs get "weak", and base entries override everything.

Return type:

dict[tuple[int, int], str]

class mmm_framework.continuous_learning.LearningState(channels, center, B, value, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', mode='fixed', cap=None, beta_scale=1.0, gamma_scale=0.8, spend_ref=None, discount_half_life=None, spline_prior='iid', data=None, posterior=None, history=<factory>, summaries=<factory>, geo_ids=None, stationarity_reports=<factory>)

Bases: object

Accumulated experiment data + the current posterior + decision readouts.

spend_ref (dollars per scaled unit, per channel) is carried for the dollar boundary — convert with to_dollars() / to_scaled(); everything inside the state (center, panel spend, summaries) is scaled units. geo_ids (optional, from the data dict’s "geo_ids" key) pins the geo identity across waves — the misspecification study showed the loop diverges if geo baselines are silently re-drawn under the same geo_idx. likelihood / time_effect mirror model.fit()’s knobs (defaults reproduce the old behavior exactly) and are threaded on every refit; a time_effect != "none" program requires and accumulates a global period_idx across waves (see ingest()).

channels: list[str]

center: ndarray

B: float

value: float

pairs: list[tuple[int, int]] | None = None

pair_signs: dict[tuple[int, int], str] | None = None

activation: str = 'hill'

likelihood: str = 'normal'

time_effect: str = 'none'

mode: str = 'fixed'

cap: float | None = None

beta_scale: float = 1.0

gamma_scale: float = 0.8

spend_ref: ndarray | None = None

discount_half_life: float | None = None

spline_prior: str = 'iid'

data: dict[str, Any] | None = None

posterior: Posterior | None = None

history: list[WaveRecord]

summaries: list[dict[str, Any]]

geo_ids: list[str] | None = None

stationarity_reports: list[Any]

ingest(wave_data, *, check_stationarity=False)

Append a wave’s rows to the accumulated panel (same geos throughout).

wave_data may carry an optional "geo_ids" key (list[str] of length n_geo): the first wave pins the program’s geo identities and later waves must match them exactly (order included).

check_stationarity=True (opt-in; requires the wave to carry period_idx and cell_idx) runs the within-wave stationarity guard (wave_stationarity_check()) before ingesting: the report is appended to self.stationarity_reports and a flagged mid-wave regime change raises a UserWarning — the wave is still ingested, so the caller decides whether to censor_periods() and re-ingest, or extend/repeat the wave.

When the program models a national time effect (time_effect != "none"), every wave MUST carry a "period_idx" (wave-local, 0-based); it is accumulated with a per-wave offset (shifted by the accumulated maximum + 1) so two waves’ shocks never alias onto one global tau_t. Programs with time_effect="none" ignore any period_idx a wave carries (the accumulated data dict is unchanged).

Return type:

None

ingest_summaries(items)

Append summary observations (historical lift readouts, no panel).

Each item follows the summary schema in mmm_framework.continuous_learning.model (spend_test/ spend_base in SCALED units, total lift ± se in natural KPI units, scale = geo-periods aggregated). Validated eagerly.

Return type:

None

fit(**fit_kwargs)

Refit the surface on ALL accumulated evidence (carries the posterior).

The state’s discount_half_life / spline_prior are threaded to every refit (both overridable per call via fit_kwargs). When discounting is active, per-row ages are derived from the accumulated row_week column that ingest() maintains (age = newest week − row week, in weeks).

Return type:

Posterior

recommend(**kwargs)

Return type:

ndarray

funding(alloc=None, **kwargs)

regret(**kwargs)

plan(**kwargs)

All per-wave readouts from ONE Thompson pass (see fix F4).

Delegates to plan_from_posterior() with this program’s budget/value/mode/cap; pass q/seed/ group_budgets etc. through as keyword arguments.

Return type:

PlanResult

next_design(delta, *, center=None, probe_pairs=None)

Return type:

ndarray

recenter(alloc, *, min_frac=0.05)

Move the trust-region center to alloc, floored away from zero.

The CCD is multiplicative (review fix F5): recentering a channel onto ~0 zeroes every axial/off-axis/shutoff cell for it — no variation, so beta is unidentified and the spend floor’s zero gradient means the channel can never be resurrected by the loop. Each channel is therefore clamped to >= min_frac * (B / K) in the center’s scaled units (min_frac=0 restores the old unfloored behavior). Warns whenever a clamp fires.

Return type:

None

__init__(channels, center, B, value, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', mode='fixed', cap=None, beta_scale=1.0, gamma_scale=0.8, spend_ref=None, discount_half_life=None, spline_prior='iid', data=None, posterior=None, history=<factory>, summaries=<factory>, geo_ids=None, stationarity_reports=<factory>)

class mmm_framework.continuous_learning.WaveRecord(wave, n_rows, e_regret, enbs, stop, recommendation, funded, mroas_mean, prob_above_line, profit_gap, profit_gap_rel, max_rhat, n_summaries=0, kg_used=False, chosen_delta=None)

Bases: object

A JSON-safe snapshot of one wave’s decision state.

kg_used / chosen_delta describe the designed wave whose data this record’s fit ingested: whether the Laplace knowledge-gradient selected the design (see select_next_design()) and which trust-region delta it chose. Both default to the fixed-design values so historical records (and the byte-stable default loop) are unchanged.

wave: int

n_rows: int

e_regret: float

enbs: float

stop: bool

recommendation: list[float]

funded: list[bool]

mroas_mean: list[float]

prob_above_line: list[float]

profit_gap: float

profit_gap_rel: float

max_rhat: float | None

n_summaries: int = 0

kg_used: bool = False

chosen_delta: float | None = None

__init__(wave, n_rows, e_regret, enbs, stop, recommendation, funded, mroas_mean, prob_above_line, profit_gap, profit_gap_rel, max_rhat, n_summaries=0, kg_used=False, chosen_delta=None)

mmm_framework.continuous_learning.run_closed_loop(world, *, center, B, value, n_geo=80, t_pre=6, t_test=10, delta=0.6, noise=0.6, n_holdout=0, mode='fixed', cap=None, gamma_scale=0.8, pair_signs=None, margin=1.0, population=1.0, wave_cost=0.05, max_waves=4, recenter=True, planner_q=120, fit_kwargs=None, use_laplace_kg=False, candidate_deltas=(0.3, 0.6, 0.9), kg_n_outcomes=32, stratify_geos=False, stop_patience=1, discount_half_life=None, spline_prior='iid', seed=0)

Run the full loop against a known world and return the decision trace.

Each pass: refit on all data, compute the recommendation + funding line + expected regret, evaluate the ENBS stopping rule, then (unless stopping) recenter on the recommendation, run a fresh designed wave over the SAME geos, and ingest it. Returns the per-wave WaveRecord history plus the final recommendation and the truth-optimal target for comparison.

use_laplace_kg=True (opt-in; the default path is byte-identical to the historical loop) scores candidate_deltas with select_next_design() each wave and runs the EVSI-best design instead of the fixed-delta CCD; the choice is recorded on the NEXT wave’s WaveRecord (kg_used / chosen_delta). stratify_geos=True (opt-in) blocks the geo→cell randomization on the true per-geo baselines a_geo (see assign_geos()).

stop_patience (default 1 — byte-identical to the historical loop) requires that many consecutive ENBS <= 0 evaluations before the loop stops. Evaluating ENBS every wave is NOT frequentist optional stopping — a Bayesian expected-value rule is stopping-rule-invariant under the likelihood principle (Edwards, Lindman & Savage 1963; Berger & Wolpert 1988) — but that invariance assumes a correctly specified model, and the misspecification study shows intervals can be narrow and wrong. stop_patience=2 is the cheap belt-and-suspenders guard: one optimistically-low regret read (from a temporarily overconfident, possibly misspecified posterior) cannot end the programme on its own. Each WaveRecord’s stop stays the raw per-wave verdict.

Return type:

dict[str, Any]

mmm_framework.continuous_learning.select_next_design(post, center, pairs, B, value, *, mode='fixed', cap=None, candidate_deltas=(0.3, 0.6, 0.9), candidate_probe_sets=None, n_geo=80, t_test=10, n_outcomes=32, seed=0, fallback_delta=0.6, cost_fn=None)

Score candidate CCDs with the Laplace knowledge-gradient; pick the best.

Candidates are central_composite(center, delta, probe_set) for every delta in candidate_deltas (clipped to (0, 1]) × every probe set in candidate_probe_sets (default: just pairs). Each candidate is scored with laplace_knowledge_gradient() — decision-aware EVSI, no MCMC — using the SAME seed for every candidate (common random numbers, so the Monte-Carlo argmax does not flap). The Fisher weights come from the fitted observation family (observation_unit_info()): Gaussian and Student-t posteriors need the sigma site (the noise scale), NegBinomial posteriors need phi plus the A/a_geo baseline. There is deliberately NO numeric fallback for missing sites: a summaries-only fit (whose beta/gamma live on the KPI’s natural scale under prior_scaling="auto") has no meaningful observation scale to score with, so a hard-coded guess would make every candidate’s Fisher information explode identically and the reported per-candidate scores carry no information while claiming kg_used=True.

cost_fn (optional) makes the selection cost-aware: each candidate is ranked by KG / cost_fn(cells) — decision value per dollar, the knowledge-gradient analogue of expected-improvement-per-unit-cost (EIpu; Lee, Perrone, Archambeau & Seeger 2020). Cell counts grow ~3 per extra arm plus 2 per probed pair, so a cost-blind KG systematically over-favors large designs; dividing by the design’s own cost is the greedy myopic version of the ENBS cost-benefit trade (planner.enbs), which keeps the acquisition and the stopping rule coherent. cost_fn receives the candidate’s cell matrix and must return a positive finite cost (e.g. lambda cells: fixed + per_cell * len(cells)). The raw score and the ranking score_per_cost are both recorded; cost_fn=None is byte-identical to the historical ranking.

Degrades gracefully: an activation without a transform spec, an unknown likelihood family, a posterior missing its observation-noise sites (summaries-only fit), a non-finite candidate score, a non-positive or non-finite candidate cost, a near-singular moment-matched covariance, or any scoring ValueError falls back to the fixed design central_composite(center, fallback_delta, pairs).

Returns:

on success {"kg_used": True, "chosen_delta", "chosen_probe_pairs", "kg_scores": [{"delta", "probe_pairs", "score"[, "cost", "score_per_cost"]}, ...], "sigma", "cost_aware", "surrogate"} (sigma is None for a count posterior, which has no Gaussian noise scale); on fallback {"kg_used": False, "reason", "surrogate"}. surrogate is the compact surrogate_validity() report (ok/khat/params_flagged…) — when ok is False the Gaussian surrogate is suspect for this posterior and the wave’s scores should be spot-checked against (or replaced by) the NUTS-refit knowledge_gradient().

Return type:

(cells, meta) — the chosen design array and a JSON-safe meta dict

mmm_framework.continuous_learning.world_optimal_allocation(world, B, value, *, mode='fixed', cap=None, n_starts=8, seed=0)

The truth-optimal allocation + profit (the recovery/closure target).

Return type:

tuple[ndarray, float]

mmm_framework.continuous_learning.due_for_retest(sigma_post, weeks_elapsed, channel, sigma_exp, *, half_life_overrides=None)

Whether a channel’s evidence has decayed enough to warrant a re-test.

Thin reuse of mmm_framework.planning.eig.reexperiment_due() so the continuous-learning loop’s re-test trigger uses the framework’s existing information-decay model (sigma_eff^2(t) = sigma_post^2 exp(lambda t), lambda = ln2 / half_life) rather than a private one. Returns (due, current_eig_nats).

Return type:

tuple[bool, float]

mmm_framework.continuous_learning.adstock_half_life(adstock_alpha)

Half-life (periods) of a geometric adstock with retention alpha.

alpha**h = 1/2  =>  h = ln(1/2) / ln(alpha). The advertising literature finds carryover half-lives are strongly category-specific (industry FMCG figures cluster near 2–2.5 weeks against academic estimates of 7–12 — Fry, Broadbent & Dixon 2000), which is the empirical argument against one global decay constant: use each channel’s measured carryover as the anchor for its own evidence half-life, via due_for_retest(..., half_life_overrides={channel: h}). The carryover half-life is a lower anchor — a finding decays at least as fast as the mechanism that produced it drifts, and audience/creative drift usually puts the evidence half-life at a multiple of it; when the programme has accumulated waves, prefer the measured estimate_half_life().

alpha <= 0 (no carryover) returns 0.0; alpha >= 1 (no decay) returns inf.

Return type:

float

mmm_framework.continuous_learning.estimate_half_life(gaps_weeks, shifts, sigmas)

Estimate the information-decay half-life from realized posterior drift.

The decay clock (math page Eq. 22) inflates a finding’s variance as sigma_eff^2(t) = sigma_post^2 * exp(lambda * t) with lambda = ln2 / h — but h need not be guessed: the loop observes how fast its own beliefs drift. Between two fits of the same quantity separated by a gap t_j (weeks), the decay model says the realized mean shift is a draw with variance sigma_j^2 * (exp(lambda t_j) - 1) on top of the posterior spread, so z_j = (shift_j / sigma_j)^2 has expectation exp(lambda t_j) - 1. The method-of-moments estimate solves mean_j[z_j] = mean_j[exp(lambda t_j) - 1] for lambda (monotone — a bisection). This is the adaptive-discounting idea of the dynamic-linear-model literature: West & Harrison (1997, §6.3) model each component as losing a constant information fraction delta per period (delta = exp(-lambda) here, returned as "discount"), with separate discounts per component — i.e. estimate one h per channel, not one global constant.

Parameters:

• gaps_weeks (ndarray) – per-observation elapsed weeks between the two fits.

• shifts (ndarray) – realized posterior-mean shifts of the tracked quantity (e.g. a channel’s marginal ROAS) over each gap.

• sigmas (ndarray) – the posterior sd of the quantity at the earlier fit.

Return type:

dict[str, float]

Returns:

{"half_life": h, "lambda": lam, "discount": exp(-lam), "n": J} — half_life = inf when no excess drift is observed (the findings are not decaying measurably). Requires at least one observation with positive gap and sigma.

Surface

The response surface — one differentiable JAX function, used everywhere.

This module is the single source of truth for how scaled spend drives incremental outcome in the continuous-learning loop. The very same incremental() is evaluated inside

• the NumPyro likelihood (mmm_framework.continuous_learning.model),

• the synthetic data-generating process (dgp), and

• the allocator / acquisition planner (planner, which differentiates it with jax.grad()).

Because all three call the identical function, the optimizer can never disagree with the fitted surface and the recovery test exercises the exact map the planner exploits.

Functional form (per geo-week, for a scaled spend vector s of length K):

f_c(s_c)       = s_c^alpha_c / (kappa_c^alpha_c + s_c^alpha_c)     # Hill in [0, 1)
incremental(s) = sum_c beta_c f_c(s_c)
                 + sum_{c<c'} gamma_{cc'} f_c(s_c) f_c'(s_c')

The Hill activation matches the framework’s SaturationType.HILL (x^slope / (x^slope + sat_half^slope)) with slope = alpha and sat_half = kappa — so a continuous-learning posterior is directly comparable to a BayesianMMM Hill fit on the same channel.

The interaction block uses an upper-triangular gamma matrix (gamma[i, j] for i < j); gamma carries the sign of the synergy (positive = complementarity, negative = cannibalization). The cross-partial is d^2 incremental / ds_c ds_c' = gamma_{cc'} f_c'(s_c) f_c''(... ) — see the guide for the algebra; the planner never needs it by hand because it auto-differentiates incremental() directly.

mmm_framework.continuous_learning.surface.activation(spend, kappa, alpha)

Beta-stripped Hill activation f_c(s_c) in [0, 1), elementwise.

Parameters:

• spend – scaled spend, shape (K,) (non-negative; a shutoff cell is 0).

• kappa – half-saturation per channel, shape (K,) (positive).

• alpha – Hill shape per channel, shape (K,) (positive).

Returns:

The activation fraction, shape (K,).

mmm_framework.continuous_learning.surface.incremental(spend, beta, kappa, alpha, gamma)

Incremental response to one scaled spend vector.

Parameters:

• spend – scaled spend, shape (K,).

• beta – channel ceilings, shape (K,) (non-negative).

• kappa – half-saturations, shape (K,).

• alpha – Hill shapes, shape (K,).

• gamma – upper-triangular interaction matrix, shape (K, K) (zeros on and below the diagonal).

Returns:

Scalar incremental response.

mmm_framework.continuous_learning.surface.incremental_batch(spend, beta, kappa, alpha, gamma)

Incremental response to one scaled spend vector.

Parameters:

• spend – scaled spend, shape (K,).

• beta – channel ceilings, shape (K,) (non-negative).

• kappa – half-saturations, shape (K,).

• alpha – Hill shapes, shape (K,).

• gamma – upper-triangular interaction matrix, shape (K, K) (zeros on and below the diagonal).

Returns:

Scalar incremental response.

mmm_framework.continuous_learning.surface.incremental_jit(spend, beta, kappa, alpha, gamma)

Incremental response to one scaled spend vector.

Parameters:

• spend – scaled spend, shape (K,).

• beta – channel ceilings, shape (K,) (non-negative).

• kappa – half-saturations, shape (K,).

• alpha – Hill shapes, shape (K,).

• gamma – upper-triangular interaction matrix, shape (K, K) (zeros on and below the diagonal).

Returns:

Scalar incremental response.

mmm_framework.continuous_learning.surface.response_curve(spend_matrix, beta, kappa, alpha, gamma)

Incremental response for every row of a (N, K) scaled-spend panel.

A thin numpy-friendly wrapper over incremental_batch returning a plain jax.Array of shape (N,).

mmm_framework.continuous_learning.surface.logistic(spend, lam)

Exponential-saturation activation f(s) = 1 - exp(-lam * s) in [0, 1).

Smooth, strictly increasing, strictly concave (no inflection), f(0)=0. One shape parameter per channel; the half-saturation point is ln(2)/lam. Matches the framework’s SaturationType.LOGISTIC.

mmm_framework.continuous_learning.surface.hill_mixture(spend, kappa1, alpha1, kappa2, alpha2, w)

A weighted sum of two Hill curves: w f(κ1,α1) + (1-w) f(κ2,α2).

Still smooth, monotone and saturating (f(0)=0, -> a value in [0,1)), but far more flexible than a single Hill — it can express a two-phase / shoulder shape (a low-κ, high-α component that switches on early plus a high-κ component that keeps rising) that a single Hill or a logistic curve cannot represent. Used as a true DGP to study model misspecification: fit it with a single Hill (mild) or a logistic (severe) and watch what breaks.

mmm_framework.continuous_learning.surface.MSPLINE_J = 9

Number of monotone I-spline basis functions (= weight parameters w1..wJ).

mmm_framework.continuous_learning.surface.monotone_spline_basis(spend)

The MSPLINE_J monotone I-spline basis functions at spend.

Returns shape spend.shape + (MSPLINE_J,); column j is I_{j+1}(spend / MSPLINE_S_MAX), rising monotonically 0 -> 1 (earlier columns rise earlier). Useful for plotting and for understanding what the fitted weights mean.

mmm_framework.continuous_learning.surface.monotone_spline(spend, w1, w2, w3, w4, w5, w6, w7, w8, w9)

Shape-agnostic monotone-spline activation, values in [0, 1).

f(s) = sum_j w_j I_j(s / MSPLINE_S_MAX) / sum_j w_j with positive weights w1..w9 (one per I-spline basis function; only the ratios matter). Smooth, strictly monotone on [0, MSPLINE_S_MAX), f(0) = 0, fully saturated at/after MSPLINE_S_MAX. Early weights buy an early, steep rise; late weights buy a slow, late one — mixtures express two-phase and plateau shapes no single Hill or logistic can.

mmm_framework.continuous_learning.surface.surface_value(spend, beta, gamma, act_fn, shape)

Incremental response for one spend vector under an arbitrary activation.

act_fn(spend, *shape) -> f computes the per-channel activations; shape is a tuple of (K,) arrays in the activation’s parameter order. The interaction block is unchanged, so any activation drops straight in.

mmm_framework.continuous_learning.surface.surface_over_rows(spend_matrix, beta, gamma, act_fn, shape)

surface_value() for every row of a (N, K) panel -> (N,).

Model

Layer 1 — the model-free response model (NumPyro).

A lightweight Bayesian response surface fit directly from designed experiment data — no observational time series, no pre-fit MMM. The whole point of the continuous-learning loop is that the experiment’s designed cross-sectional variation identifies the surface, so the priors inform but the data dominates.

Priors (guide §4.2; y_loc/y_scale default to 0/1):

beta_c   ~ HalfNormal(beta_scale * y_scale)        # ceilings on the y scale
kappa_c  ~ LogNormal(0, 0.5)                       # O(1) half-saturation
alpha_c  ~ TruncatedNormal(2, 0.5, [0.5, 5])       # Hill shape
gamma    ~ sign-informed per pair (see PAIR_SIGNS) # synergy (× y_scale)
A        ~ Normal(y_loc, 5 * y_scale)              # geo-intercept hypers
sigma_a  ~ HalfNormal(2 * y_scale)
a_geo    ~ Normal(A, sigma_a)                      # baseline random intercept
sigma    ~ HalfNormal(y_scale)                     # observation noise

The data contract mandates y in NATURAL units (never normalized), so the prior scales cannot be hard-coded O(1): fit() derives y_loc/y_scale from the evidence (prior_scaling="auto", the default — y_loc = mean(y) and y_scale = `` the nearest DECADE of ``std(y), or of the |lift|/scale magnitudes for summaries-only fits) so a revenue-scale KPI (1e5+) gets priors on its own order of magnitude while O(1) data keeps the original priors exactly (y_scale = 1). prior_scaling="unit" forces the original O(1) priors regardless of the data.

Likelihood:

y ~ Normal(a_geo[geo_idx] + incremental(spend), sigma)

The geo intercept a_geo is pinned by the pre-period (where every geo shares the status-quo allocation), which breaks the within-geo collinearity between baseline and incremental response — the CUPED-style identification requirement of guide §3.2.

Two opt-in extensions (the defaults reproduce the graph above byte-identically — no new sample sites, same rng stream):

• likelihood="studentt" — heavy-tailed robust observation for continuous KPIs with occasional wild geo-weeks (data glitches, local shocks). Same location/scale structure as the Gaussian plus a learned tail nu ~ Gamma(2, 0.1) (the Juárez–Steel robustness prior); y ~ StudentT(nu, mu, sigma). Everything downstream that reads only the surface parameters is unchanged; the sigma site persists (it is the t SCALE — the sd is sigma * sqrt(nu/(nu-2))).

• likelihood="negbinomial" — count KPIs. The panel observation becomes y ~ NegativeBinomial2(softplus(mu), phi) with concentration phi ~ LogNormal(log 30, 2) replacing the sigma site; y must be non-negative integer counts (_validate_panel enforces it — note preprocess.cuped_adjust() mutates y to non-integer/possibly-negative values and is therefore incompatible). The summary block stays Gaussian for BOTH likelihoods: a historical lift ± se readout is already a normal-approximation aggregate, so a count observation model has nothing to add there. Link caveat: the planner’s marginal readouts differentiate the LATENT surface R while the observable count mean is softplus(mu) (derivative sigmoid(mu)) — asymptotically exact for mu >> 1, an overstatement of the count-scale marginal near mu ≈ 0 (see fit()).

• time_effect="national" — a zero-centered hierarchical national per-period shock tau_t ~ Normal(0, sigma_tau), sigma_tau ~ HalfNormal(y_scale), added to every panel row of period t (requires data["period_idx"]). Zero-centered partial pooling, NOT fixed effects: a free tau level is exactly collinear with the intercept hyper A. Summaries need no tau — a national per-period constant cancels in the lift difference exactly as the geo intercept does, so the evidence/summary path is unchanged.

Summary observations (past experiments, no panel required). A historical lift test is a difference measurement, so the geo intercept cancels and no pre-period is needed — the atomic evidence unit is:

{"spend_test": (K,), "spend_base": (K,), "lift": float, "se": float,
 "scale": float}   # scale = n_units * n_periods the total lift aggregates

with likelihood lift ~ Normal(scale * (R(spend_test) - R(spend_base)), se). fit() accepts data["summaries"] alongside the panel, or instead of it ({"summaries": [...], "n_geo": 0}) — a team with historical readouts and no panel can still fit the surface. Same structural-stationarity caveat as the MMM’s off-panel calibration: the response curve is assumed stable across the periods the summaries span. A handful of summaries constrains only a few contrasts of the surface, so κ/α stay prior-dominated — trust the funded set, not the curve shape.

mmm_framework.continuous_learning.model.VALID_LIKELIHOODS = ('normal', 'negbinomial', 'studentt')

Observation families for the panel likelihood (“normal” is the default and reproduces the original graph byte-identically).

mmm_framework.continuous_learning.model.VALID_TIME_EFFECTS = ('none', 'national')

National per-period time-effect modes (“none” default = no tau sites).

mmm_framework.continuous_learning.model.default_pairs(n_channels)

All upper-triangular channel pairs (i, j) with i < j.

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.model.normalize_pair_signs(pair_signs)

Normalize pair-sign keys to (min, max) orientation and merge.

The model’s pairs are always upper-triangular (default_pairs()) and the prior lookup is pair_signs.get((i, j)) — a reversed key like (1, 0) would otherwise validate but be silently ignored (the user’s explicit prior replaced by "weak"). Duplicate entries that agree are merged; conflicting duplicates ((0, 1) vs (1, 0) with different signs) raise.

Return type:

dict[tuple[int, int], str]

mmm_framework.continuous_learning.model.pair_name(channels, pair)

Posterior site name for a pair’s interaction, gamma_<ci>_<cj>.

Return type:

str

mmm_framework.continuous_learning.model.demote_channel(channels, name, *, pairs=None, base=None)

Demote a non-randomizable channel (a walled garden) to main-effect-only.

Every interaction that touches name becomes "zero" (prior-dominated); its main effect is still identified. Flag gamma(., name) as the least-trustworthy parameter in any decision (guide §5.4).

Return type:

dict[tuple[int, int], str]

mmm_framework.continuous_learning.model.probe_pairs_excluding(channels, name, *, pairs=None)

Drop the off-axis design cells for a demoted channel’s pairs.

Return type:

list[tuple[int, int]]

mmm_framework.continuous_learning.model.model(spend, geo_idx, n_geo, y=None, *, channels, pairs, pair_signs, beta_scale=1.0, gamma_scale=0.8, activation='hill', likelihood='normal', time_effect='none', period_idx=None, n_period=0, y_loc=0.0, y_scale=1.0, noise_mult=None, spline_prior='iid', summary_test=None, summary_base=None, summary_scale=None, summary_lift=None, summary_se=None)

NumPyro generative model (guide §4.2/4.3), for any pluggable activation.

noise_mult (optional, (N,), all >= 1) is the per-row information-discount multiplier exp(0.5 * lambda * age) computed by fit(): for the Gaussian families the observation scale becomes sigma * noise_mult (old rows are noisier — exactly the Eq. 22 sigma_eff semantics moved into the likelihood), and for the count family (which has no noise scale) the same discount is applied as a power likelihood, log p scaled by 1 / noise_mult^2 = exp(-lambda * age) (the two coincide in how they down-weight the data term). None is the historical, byte-identical graph. spline_prior selects the monotone-spline weight prior (see _sample_activation_shape()).

y=None draws from the prior predictive (prior checks). pair_signs maps each pair to a prior family; pairs absent from the map default to "weak". Sign-constrained pairs (neg/pos) record a signed deterministic gamma_<ci>_<cj> so every pair has a uniform posterior site. activation selects the per-channel saturation family ("hill" default, "logistic" for exponential saturation).

likelihood selects the panel observation family: "normal" (default — the original graph, byte-identical), "studentt" for heavy-tailed robust continuous KPIs (adds a tail df nu ~ Gamma(2, 0.1) next to sigma; y ~ StudentT(nu, mu, sigma)), or "negbinomial" for count KPIs (phi ~ LogNormal(log 30, 2) replaces the sigma site and y ~ NegativeBinomial2(softplus(mu), phi); the identity-link mu composition a_geo[geo_idx] + surface (+ tau) is unchanged, with a softplus positivity guard). The summary block below stays Gaussian for both likelihoods: a lift ± se is a normal-approximation aggregate, and the geo intercept AND any national per-period constant cancel in the difference.

time_effect="national" (with period_idx (N,) / n_period) adds a zero-centered hierarchical national shock tau_t ~ Normal(0, sigma_tau) to every panel row of period t. Zero-centered partial pooling, NOT fixed effects — a free tau level is exactly collinear with the intercept hyper A (both are national constants); the hierarchical prior resolves it softly and the pre-period still pins a_geo. The tau sites are only sampled on the opt-in path, so the default graph is untouched.

y_loc/y_scale (defaults 0/1 = the original O(1) priors, byte-identical graph) put the intercept/noise/effect priors on the KPI’s natural scale: A ~ N(y_loc, 5·y_scale), sigma_a ~ HalfNormal(2·y_scale), sigma ~ HalfNormal(y_scale), and the beta/gamma scales multiply by y_scale too — the incremental response must reach y-magnitude, so the effect ceilings live on the y scale. fit() derives them from the data.

The panel likelihood is gated on spend.shape[0] > 0 — a summaries-only fit skips the geo-intercept plate and observation-noise sites entirely. The optional summary block (summary_test (M, K), summary_base (M, K), summary_scale (M,), summary_lift (M,), summary_se (M,)) adds one Gaussian observation per historical lift readout:

lift_m ~ Normal(scale_m * (R(test_m) - R(base_m)), se_m)

The geo intercept cancels in the difference, so summaries need no pre-period. Works with any activation.

class mmm_framework.continuous_learning.model.Posterior(samples, channels, pairs, pair_signs=<factory>, activation='hill', likelihood='normal', time_effect='none', spend_ref=None, diagnostics=<factory>)

Bases: object

Posterior draws plus the metadata the planner needs.

samples holds plain numpy arrays keyed by site name (beta, kappa, alpha, A, sigma_a, a_geo, sigma, and one gamma_<ci>_<cj> per pair; the intercept/noise sites are absent for a summaries-only fit). A likelihood="negbinomial" fit carries phi (the NB concentration) instead of sigma; a likelihood="studentt" fit carries nu (the tail df) alongside sigma; a time_effect="national" fit adds sigma_tau and tau ((draws, n_period)). spend_ref is the per-channel reference constant used to scale spend — convert with to_dollars() / to_scaled().

samples: dict[str, ndarray]

channels: list[str]

pairs: list[tuple[int, int]]

pair_signs: dict[tuple[int, int], str]

activation: str = 'hill'

likelihood: str = 'normal'

time_effect: str = 'none'

spend_ref: ndarray | None = None

diagnostics: dict[str, Any]

property n_channels: int

property n_draws: int

gamma_matrix(d)

Assemble the (K, K) upper-triangular gamma matrix for draw d.

Return type:

ndarray

property shape_names: tuple[str, ...]

The activation’s shape-parameter site names, in order.

draw_params(d)

Per-draw params for the surface functions, activation-agnostic.

Returns {beta, gamma, shape, act_fn, activation} where shape is a tuple of (K,) arrays in the activation’s parameter order and act_fn is its JAX activation. (For the Hill default the shape is (kappa, alpha).)

Return type:

dict[str, Any]

gamma_summary()

Per-pair posterior mean and 5/95 percentiles of the synergy.

Return type:

dict[str, dict[str, float]]

__init__(samples, channels, pairs, pair_signs=<factory>, activation='hill', likelihood='normal', time_effect='none', spend_ref=None, diagnostics=<factory>)

mmm_framework.continuous_learning.model.fit(data, *, channels, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', beta_scale=1.0, gamma_scale=0.8, num_warmup=500, num_samples=500, num_chains=2, seed=0, progress_bar=False, spend_ref=None, prior_scaling='auto', discount_half_life=None, spline_prior='iid')

Fit the response surface to a geo-week panel and/or summary readouts.

Parameters:

• data (dict[str, Any]) – the data contract — {"spend": (N, K), "geo_idx": (N,), "n_geo": int, "y": (N,)} (scaled spend, natural-unit y), optionally with "summaries": list[dict] (see the module header), an optional "period_idx": (N,) (0-based int; required when time_effect="national"), or summaries-only ({"summaries": [...], "n_geo": 0} / an empty (0, K) panel). At least one of panel rows / summaries required.

• channels (list[str]) – channel names, length K.

• pairs (list[tuple[int, int]] | None) – interaction pairs to model (default: all upper-triangular pairs).

• pair_signs (dict[tuple[int, int], str] | None) – sign-informed prior family per pair (default: all "weak"). Keys are normalized to (min, max) orientation; conflicting duplicate entries raise.

• likelihood (str) – panel observation family — "normal" (default, byte-identical to the original graph), "studentt" for heavy-tailed robust continuous KPIs (a learned tail df nu joins sigma; outlier geo-weeks stop dragging the surface), or "negbinomial" for count KPIs (y must be non-negative integer counts; the summary block stays Gaussian either way). Link caveat for the planner: with "negbinomial" the marginal readouts (marginal_roas, expected_regret, free-mode allocation) differentiate the LATENT surface R, while the observable count mean is softplus(mu) whose derivative is sigmoid(mu) — so the count-scale marginal response is sigmoid(mu) * dR/ds. Since sigmoid(mu) ≈ 1 whenever mu >> 1 (typical counts), the readouts are asymptotically exact; near mu ≈ 0 (per-row means of a few counts) they OVERSTATE the count-scale marginal response by up to 1/sigmoid(mu). fit() warns loudly when the observed count mean is below ~20. Fixed-budget Thompson argmax is unaffected (softplus is monotone).

• time_effect (str) – "none" (default) or "national" — a zero-centered hierarchical per-period national shock tau_t added to the panel mean (a non-empty panel requires data["period_idx"]; see model()). A summaries-only fit needs NO period identity: there are no panel rows for tau_t to index and a national per-period constant cancels in the lift difference, so the summary path is unchanged (no tau sites are sampled).

• gamma_scale (float) – interaction prior scale — the most consequential knob; audit it with prior_sensitivity() (guide §8.2).

• spend_ref (ndarray | None) – per-channel reference constant used to scale spend, carried on the returned Posterior for dollar mapping (see to_dollars()).

• prior_scaling (str) – "auto" (default) derives the intercept/noise/effect prior scales from the evidence so natural-unit KPIs at any magnitude fit sanely; "unit" reproduces the original O(1) priors exactly (see _resolve_prior_scaling()).

• discount_half_life (float | None) – information-decay half-life h in WEEKS (None — the default — is the historical static fit, every row weighing the same forever). When set, each panel row’s likelihood is discounted by its age: the observation scale is inflated to sigma * exp(0.5 * lambda * age) with lambda = ln2 / h — the math page’s Eq. 22 decay clock moved INTO the likelihood (West–Harrison discounting; the count family power-scales the log-likelihood by exp(-lambda * age) instead, since it has no noise scale). Row ages come from data["row_age"] ((N,) weeks, 0 = newest — LearningState maintains this automatically) or, failing that, are derived from data["period_idx"] as max - period; a discounted panel fit with neither raises. Summaries may carry an "age_weeks" key — their SE is inflated by the same clock. The effect: effective sample size SATURATES instead of growing without bound, so the posterior keeps an honest variance floor when media behaviour drifts (a static fit’s bands shrink like 1/sqrt(rows) while the truth moves — narrow and wrong). Use estimate_half_life() to measure h from the programme’s own wave-to-wave drift; a measured lambda ~ 0 recovers full pooling.

• spline_prior (str) – monotone-spline weight prior — "iid" (default, the historical independent LogNormal(0, 1) weights, byte-identical) or "pspline": a first-order random walk on the log-weights with a learned per-channel smoothness w_tau (P-spline / adaptive-smoothing construction). tau -> 0 collapses to the neutral near-linear ramp, so the EFFECTIVE basis dimension grows only as the designed cells earn it — the adaptive answer to the fixed-basis approximation-bias floor (math page Eq. 24) without hand-tuning the knot count. Only valid with activation="monotone_spline"; the public w1..wJ sites are unchanged (the acquisition layer and planner need no changes).

Return type:

Posterior

Returns:

A Posterior with merged-chain numpy samples, R-hat/ESS diagnostics on the key parameters, and diagnostics["evidence"] = {"n_rows", "n_summaries"}.

mmm_framework.continuous_learning.model.refit_fn_from_data(base_data, *, channels, pairs=None, pair_signs=None, activation='hill', likelihood='normal', beta_scale=1.0, gamma_scale=0.8, num_warmup=200, num_samples=200, num_chains=1, seed=7, prior_scaling='auto', time_effect='none')

Build a refit_fn(extra_spend, extra_geo_idx, extra_y) -> Posterior.

The knowledge-gradient acquisition (planner.knowledge_gradient()) fantasizes wave outcomes and refits; this closure appends the fantasy rows to base_data and runs a short NUTS chain. This is the expensive path — see guide §9.1 for the Laplace-update replacement in production.

Requires a panel base dataset: fantasy observations are geo-week rows, so a summaries-only data dict has nothing to append them to.

likelihood must match the posterior the knowledge gradient is scoring: the fantasy y values are drawn from that same observation family (Gaussian, Student-t, or NB counts), so refitting them under a different family would score the design against a model nobody will run.

time_effect other than "none" raises NotImplementedError: the knowledge-gradient fantasy rows carry no period identity, so a national tau_t cannot be refit on the augmented panel. A period_idx present in base_data is DROPPED from the refit base for the same reason (the refit models no time effect, so the index is inert either way).

Return type:

Callable[[ndarray, ndarray, ndarray], Posterior]

Design

Layer 2 — experimental design (central-composite geo cells).

A wave is a designed batch of geo cells run for a fixed window. The central-composite design (CCD) is the minimal set of scaled allocations that identifies the local response surface around an operating point:

• 1 center cell — the current operating allocation (the trust-region center).

• 2K axial cells — each channel scaled to (1 +/- delta) while the others hold at center; these give the gradient / main effects.

• 2 off-axis cells per probed pair — two channels moved jointly (1 +/- delta); the only way to recover d^2 R / ds_c ds_c' (the synergy gamma). Probe the decision-pivotal pairs only.

• K shutoff cells — one channel set to 0. These isolate the remaining terms and break the beta/gamma collinearity; without them beta attenuates.

delta is a multiplicative trust-region step (a spend-variation fraction): delta = 0.6 moves a channel +/-60% off its center. A shutoff is the -100% corner. Working multiplicatively keeps the design scale-free, so it behaves the same whether spend is scaled to O(1) or left in dollars.

mmm_framework.continuous_learning.design.central_composite(center, delta, probe_pairs)

Build the CCD cells (rows = cells, columns = channels).

Parameters:

• center (ndarray) – the operating allocation, shape (K,) (scaled spend).

• delta (float) – multiplicative trust-region step in (0, 1] (e.g. 0.6).

• probe_pairs (list[tuple[int, int]]) – channel pairs to identify the cross-partial for.

Return type:

ndarray

Returns:

A (n_cells, K) array of non-negative scaled allocations, with n_cells = 1 + 2K + 2 * len(probe_pairs) + K.

mmm_framework.continuous_learning.design.assign_geos(design, n_geo, rng, *, n_holdout=0, center=None, baseline=None)

Assign CCD cells to geos — shuffled round-robin, or stratified/blocked.

Without a baseline (the default, byte-identical to the historical behavior) cells are tiled round-robin and shuffled: balanced cell counts, zero covariate awareness, and the first n_holdout geos (positional) become holdouts.

With a baseline (a per-geo covariate such as the pre-period KPI level, matched-market style) the assignment is a classic blocked randomization: geos are sorted by baseline, walked in blocks of n_cells, and each block receives a random permutation of the cell indices (the ragged tail gets a random subset without replacement). Every cell’s geos are then spread evenly across the baseline distribution, so between-cell baseline means are nearly equal. Holdouts are carved FIRST and are stratum-aware too: exactly n_holdout evenly spaced positions in baseline-sorted order, so the status-quo counterfactual spans the baseline range AND honors the requested count (a strided pick would silently under-deliver whenever n_holdout does not divide n_geo evenly).

Parameters:

• design (ndarray) – CCD cells, shape (n_cells, K).

• n_geo (int) – number of geos.

• rng (Generator) – a numpy random generator (caller owns the seed).

• n_holdout (int) – hold n_holdout geos at center for the test window (a status-quo counterfactual). Requires center.

• center (ndarray | None) – the status-quo allocation for holdout geos, shape (K,).

• baseline (Union[ndarray, Sequence[float], None]) – optional per-geo covariate, length n_geo, positionally aligned with the geo indices (the caller resolves geo ids). None keeps the legacy shuffled round-robin path.

Return type:

tuple[ndarray, ndarray]

Returns:

(geo_alloc, cell_idx) where geo_alloc is (n_geo, K) (each geo’s test allocation) and cell_idx is (n_geo,) (the design-row index, or -1 for a holdout geo).

Planner

Layer 4 — planner / acquisition.

The decision is max_s  value * R(s) - 1^T s under one of two regimes:

• fixed budget (mode="fixed"): 1^T s = B — a pure simplex.

• free budget (mode="free"): 0 <= s_c <= cap — the total is itself a decision; each channel self-funds until its marginal ROAS hits 1.

Every routine here differentiates the same incremental surface the likelihood fits (via mmm_framework.continuous_learning.surface.grad_incremental), so the allocator can never disagree with the fitted model. The surface is non-concave (negative gamma), so the allocator multi-starts and keeps the best.

Acquisition / readouts:

• plan_from_posterior() — ONE Thompson pass producing every per-wave readout coherently (recommendation, funding line at it, warm-started regret from it) as a PlanResult; prefer it over calling the pieces below separately, which each re-sample.

• thompson_wave() — a posterior over the optimal split (the spread is the exploration signal; the mean is the recommendation).

• marginal_roas() — the funding line: a channel is funded where P(value * dR/ds_c > 1) > 0.5.

• expected_regret() — E[regret] from posterior uncertainty; with enbs() it drives the stopping rule.

• knowledge_gradient() — decision-aware EVSI for scoring candidate test designs (the expensive path; see guide §7.5/9.1).

Import-light apart from JAX (shared with the model) and SciPy SLSQP.

mmm_framework.continuous_learning.planner.response_grid(post, spend_matrix, draws)

Incremental response for every (spend row, draw) pair -> (G, D).

Activation-agnostic (reads post.activation): the visualizations use this to build the posterior mean/uncertainty surfaces for any activation family, so the same animation code renders a Hill world or a logistic one.

Return type:

ndarray

mmm_framework.continuous_learning.planner.allocate_under_sample(params, B, value, *, x0=None, n_starts=3, seed=0, mode='fixed', cap=None, group_budgets=None)

Optimal allocation + profit under one posterior draw’s params.

group_budgets fixes sub-budgets over channel groups (arms): [(indices, B_g), ...] adds sum(s[indices]) == B_g per group.

Return type:

tuple[ndarray, float]

mmm_framework.continuous_learning.planner.thompson_wave(post, B, value, *, q=300, seed=0, mode='fixed', cap=None, n_starts=3, group_budgets=None)

A posterior over the optimal split — solve the allocation for q draws.

Returns (allocs, profits, draws): allocs is (q, K) (the mean is the recommended allocation, the spread is the exploration signal), profits is (q,) (each draw’s optimum under itself), and draws are the sampled draw indices.

Return type:

tuple[ndarray, ndarray, ndarray]

mmm_framework.continuous_learning.planner.recommend_allocation(post, B, value, **kwargs)

The recommended allocation — the Thompson posterior mean split.

Return type:

ndarray

mmm_framework.continuous_learning.planner.marginal_roas(post, alloc, value, *, q=300, seed=1)

Posterior of value * dR/ds_c at alloc — the funding line.

Returns (mean_mroas, prob_above_line, draws_mroas). A channel is funded where prob_above_line > 0.5.

Return type:

tuple[ndarray, ndarray, ndarray]

mmm_framework.continuous_learning.planner.expected_regret(post, B, value, *, q=300, seed=2, mode='fixed', cap=None, n_starts=3, group_budgets=None)

Expected regret of acting on the consensus allocation (guide §7.4).

regret_d = profit(best-for-draw-d under d) - profit(consensus under d) >= 0 — the profit still on the table from posterior uncertainty. Each draw’s “best” is max(cold multistart, a solve warm-started FROM the consensus, profit at the consensus) (review fix F3): the warm-started second pass catches per-draw optima the cold multistart missed near the consensus, so the regret is less biased low and the ENBS stop fires less prematurely. Returns (E[regret], consensus_alloc, alloc_sd, optimal_profit_sd).

Return type:

tuple[float, ndarray, ndarray, float]

mmm_framework.continuous_learning.planner.posterior_optimal_allocation(post, B, value, *, q=200, seed=3, mode='fixed', cap=None, n_starts=4, group_budgets=None)

argmax_a E_post[profit(a)] and its value — optimize the mean surface.

Return type:

tuple[ndarray, float]

class mmm_framework.continuous_learning.planner.PlanResult(channels, B, value, mode, recommendation, consensus, allocs, profits, draws, mroas_mean, prob_above_line, mroas_draws, e_regret, alloc_sd, profit_sd)

Bases: object

Every per-wave decision readout, derived from ONE Thompson sample.

Historically recommend/funding/regret each re-sampled with a different seed (planner seeds 0/1/2), so the funded set and the reported regret referred to different consensus vectors and paid 2-3x the SLSQP cost. Here all readouts share the same q draws: recommendation is the Thompson mean (rescaled onto the budget simplex in fixed mode), consensus is that same vector, the funding line is evaluated at it, and the regret pass is warm-started from it.

channels: list[str]

B: float

value: float

mode: str

recommendation: ndarray

consensus: ndarray

allocs: ndarray

profits: ndarray

draws: ndarray

mroas_mean: ndarray

prob_above_line: ndarray

mroas_draws: ndarray

e_regret: float

alloc_sd: ndarray

profit_sd: float

to_dict()

A small JSON-safe snapshot (summary stats, not the draw matrices).

Return type:

dict[str, Any]

__init__(channels, B, value, mode, recommendation, consensus, allocs, profits, draws, mroas_mean, prob_above_line, mroas_draws, e_regret, alloc_sd, profit_sd)

mmm_framework.continuous_learning.planner.plan_from_posterior(post, B, value, *, q=300, seed=0, mode='fixed', cap=None, n_starts=3, group_budgets=None)

One Thompson pass producing every decision readout coherently (fix F4).

Runs thompson_wave() once; the recommendation is the draw-mean (rescaled to the simplex in fixed mode); the funding line and the warm-started regret pass reuse the SAME draw indices, so the funded set, the consensus, and E[regret] all describe one allocation.

Return type:

PlanResult

mmm_framework.continuous_learning.planner.enbs(e_regret, *, margin, population, wave_cost)

Expected net benefit of sampling = E[regret] * margin * population - cost.

E[regret] is profit-per-affected-unit on the response-curve scale; margin and population convert it to a total dollar value of resolving the uncertainty, which a wave must beat to be worth running.

Return type:

float

mmm_framework.continuous_learning.planner.should_stop(e_regret, *, margin, population, wave_cost)

(stop, enbs) — stop when no wave’s expected value clears its cost.

Return type:

tuple[bool, float]

mmm_framework.continuous_learning.planner.knowledge_gradient(post, candidate_design, refit_fn, B, value, *, n_fantasy=10, t_test=10, n_geo=None, noise=None, mode='fixed', cap=None, q=120, seed=4)

One-step-lookahead EVSI for a candidate test design (guide §7.5).

Return type:

float

``KG(d) = E_y[ max_a profit(a | posterior updated with fantasised y) ]

• max_a profit(a | current posterior)``

For each fantasy: draw params, simulate the candidate design’s outcomes, refit via refit_fn(extra_spend, extra_geo_idx, extra_y) -> Posterior, and re-optimize. Expensive — refit_fn runs a (short) NUTS chain per fantasy; in production swap it for a Laplace/conjugate update (guide §9.1).

noise=None (default) fantasizes with the posterior mean of sigma — the fitted observation noise — rather than a hard-coded constant. Fantasy outcomes are drawn from the posterior’s own observation family (_fantasy_outcomes()): Gaussian, Student-t (posterior-mean nu), or NB counts (posterior-mean phi; noise ignored). refit_fn must refit under that same family — build it with model.refit_fn_from_data()’s matching likelihood.

Requires a panel-fitted posterior: fantasies simulate geo-week outcomes off the geo intercepts (a_geo or the A/sigma_a hypers), which a summaries-only fit never samples.

Loop

The outer control loop: carry the posterior across waves and know when to stop.

LearningState accumulates every wave’s data and refits the response surface on all of it each pass — that is how the posterior is “carried” across waves (each wave borrows strength from all prior data). It exposes the per-wave decision readouts (recommend / funding line / regret / stop) over the current posterior.

run_closed_loop() drives a TrueWorld through the full loop — fit -> score -> run designed wave -> update -> ENBS stop — and is both the runnable demo and the closure/stopping test backbone. In a real deployment you would replace the synthetic simulate_wave collector with the actual geo holdout results; nothing else changes.

Re-test scheduling reuses the framework’s existing information-decay model (mmm_framework.planning.eig.reexperiment_due()) so a continuous-learning program and a model-anchored program agree on when evidence has gone stale.

mmm_framework.continuous_learning.loop.world_optimal_allocation(world, B, value, *, mode='fixed', cap=None, n_starts=8, seed=0)

The truth-optimal allocation + profit (the recovery/closure target).

Return type:

tuple[ndarray, float]

class mmm_framework.continuous_learning.loop.WaveRecord(wave, n_rows, e_regret, enbs, stop, recommendation, funded, mroas_mean, prob_above_line, profit_gap, profit_gap_rel, max_rhat, n_summaries=0, kg_used=False, chosen_delta=None)

Bases: object

A JSON-safe snapshot of one wave’s decision state.

kg_used / chosen_delta describe the designed wave whose data this record’s fit ingested: whether the Laplace knowledge-gradient selected the design (see select_next_design()) and which trust-region delta it chose. Both default to the fixed-design values so historical records (and the byte-stable default loop) are unchanged.

wave: int

n_rows: int

e_regret: float

enbs: float

stop: bool

recommendation: list[float]

funded: list[bool]

mroas_mean: list[float]

prob_above_line: list[float]

profit_gap: float

profit_gap_rel: float

max_rhat: float | None

n_summaries: int = 0

kg_used: bool = False

chosen_delta: float | None = None

__init__(wave, n_rows, e_regret, enbs, stop, recommendation, funded, mroas_mean, prob_above_line, profit_gap, profit_gap_rel, max_rhat, n_summaries=0, kg_used=False, chosen_delta=None)

class mmm_framework.continuous_learning.loop.LearningState(channels, center, B, value, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', mode='fixed', cap=None, beta_scale=1.0, gamma_scale=0.8, spend_ref=None, discount_half_life=None, spline_prior='iid', data=None, posterior=None, history=<factory>, summaries=<factory>, geo_ids=None, stationarity_reports=<factory>)

Bases: object

Accumulated experiment data + the current posterior + decision readouts.

spend_ref (dollars per scaled unit, per channel) is carried for the dollar boundary — convert with to_dollars() / to_scaled(); everything inside the state (center, panel spend, summaries) is scaled units. geo_ids (optional, from the data dict’s "geo_ids" key) pins the geo identity across waves — the misspecification study showed the loop diverges if geo baselines are silently re-drawn under the same geo_idx. likelihood / time_effect mirror model.fit()’s knobs (defaults reproduce the old behavior exactly) and are threaded on every refit; a time_effect != "none" program requires and accumulates a global period_idx across waves (see ingest()).

channels: list[str]

center: ndarray

B: float

value: float

pairs: list[tuple[int, int]] | None = None

pair_signs: dict[tuple[int, int], str] | None = None

activation: str = 'hill'

likelihood: str = 'normal'

time_effect: str = 'none'

mode: str = 'fixed'

cap: float | None = None

beta_scale: float = 1.0

gamma_scale: float = 0.8

spend_ref: ndarray | None = None

discount_half_life: float | None = None

spline_prior: str = 'iid'

data: dict[str, Any] | None = None

posterior: Posterior | None = None

history: list[WaveRecord]

summaries: list[dict[str, Any]]

geo_ids: list[str] | None = None

stationarity_reports: list[Any]

ingest(wave_data, *, check_stationarity=False)

Append a wave’s rows to the accumulated panel (same geos throughout).

wave_data may carry an optional "geo_ids" key (list[str] of length n_geo): the first wave pins the program’s geo identities and later waves must match them exactly (order included).

check_stationarity=True (opt-in; requires the wave to carry period_idx and cell_idx) runs the within-wave stationarity guard (wave_stationarity_check()) before ingesting: the report is appended to self.stationarity_reports and a flagged mid-wave regime change raises a UserWarning — the wave is still ingested, so the caller decides whether to censor_periods() and re-ingest, or extend/repeat the wave.

When the program models a national time effect (time_effect != "none"), every wave MUST carry a "period_idx" (wave-local, 0-based); it is accumulated with a per-wave offset (shifted by the accumulated maximum + 1) so two waves’ shocks never alias onto one global tau_t. Programs with time_effect="none" ignore any period_idx a wave carries (the accumulated data dict is unchanged).

Return type:

None

ingest_summaries(items)

Append summary observations (historical lift readouts, no panel).

Each item follows the summary schema in mmm_framework.continuous_learning.model (spend_test/ spend_base in SCALED units, total lift ± se in natural KPI units, scale = geo-periods aggregated). Validated eagerly.

Return type:

None

fit(**fit_kwargs)

Refit the surface on ALL accumulated evidence (carries the posterior).

The state’s discount_half_life / spline_prior are threaded to every refit (both overridable per call via fit_kwargs). When discounting is active, per-row ages are derived from the accumulated row_week column that ingest() maintains (age = newest week − row week, in weeks).

Return type:

Posterior

recommend(**kwargs)

Return type:

ndarray

funding(alloc=None, **kwargs)

regret(**kwargs)

plan(**kwargs)

All per-wave readouts from ONE Thompson pass (see fix F4).

Delegates to plan_from_posterior() with this program’s budget/value/mode/cap; pass q/seed/ group_budgets etc. through as keyword arguments.

Return type:

PlanResult

next_design(delta, *, center=None, probe_pairs=None)

Return type:

ndarray

recenter(alloc, *, min_frac=0.05)

Move the trust-region center to alloc, floored away from zero.

The CCD is multiplicative (review fix F5): recentering a channel onto ~0 zeroes every axial/off-axis/shutoff cell for it — no variation, so beta is unidentified and the spend floor’s zero gradient means the channel can never be resurrected by the loop. Each channel is therefore clamped to >= min_frac * (B / K) in the center’s scaled units (min_frac=0 restores the old unfloored behavior). Warns whenever a clamp fires.

Return type:

None

__init__(channels, center, B, value, pairs=None, pair_signs=None, activation='hill', likelihood='normal', time_effect='none', mode='fixed', cap=None, beta_scale=1.0, gamma_scale=0.8, spend_ref=None, discount_half_life=None, spline_prior='iid', data=None, posterior=None, history=<factory>, summaries=<factory>, geo_ids=None, stationarity_reports=<factory>)

mmm_framework.continuous_learning.loop.select_next_design(post, center, pairs, B, value, *, mode='fixed', cap=None, candidate_deltas=(0.3, 0.6, 0.9), candidate_probe_sets=None, n_geo=80, t_test=10, n_outcomes=32, seed=0, fallback_delta=0.6, cost_fn=None)

Score candidate CCDs with the Laplace knowledge-gradient; pick the best.

Candidates are central_composite(center, delta, probe_set) for every delta in candidate_deltas (clipped to (0, 1]) × every probe set in candidate_probe_sets (default: just pairs). Each candidate is scored with laplace_knowledge_gradient() — decision-aware EVSI, no MCMC — using the SAME seed for every candidate (common random numbers, so the Monte-Carlo argmax does not flap). The Fisher weights come from the fitted observation family (observation_unit_info()): Gaussian and Student-t posteriors need the sigma site (the noise scale), NegBinomial posteriors need phi plus the A/a_geo baseline. There is deliberately NO numeric fallback for missing sites: a summaries-only fit (whose beta/gamma live on the KPI’s natural scale under prior_scaling="auto") has no meaningful observation scale to score with, so a hard-coded guess would make every candidate’s Fisher information explode identically and the reported per-candidate scores carry no information while claiming kg_used=True.

cost_fn (optional) makes the selection cost-aware: each candidate is ranked by KG / cost_fn(cells) — decision value per dollar, the knowledge-gradient analogue of expected-improvement-per-unit-cost (EIpu; Lee, Perrone, Archambeau & Seeger 2020). Cell counts grow ~3 per extra arm plus 2 per probed pair, so a cost-blind KG systematically over-favors large designs; dividing by the design’s own cost is the greedy myopic version of the ENBS cost-benefit trade (planner.enbs), which keeps the acquisition and the stopping rule coherent. cost_fn receives the candidate’s cell matrix and must return a positive finite cost (e.g. lambda cells: fixed + per_cell * len(cells)). The raw score and the ranking score_per_cost are both recorded; cost_fn=None is byte-identical to the historical ranking.

Degrades gracefully: an activation without a transform spec, an unknown likelihood family, a posterior missing its observation-noise sites (summaries-only fit), a non-finite candidate score, a non-positive or non-finite candidate cost, a near-singular moment-matched covariance, or any scoring ValueError falls back to the fixed design central_composite(center, fallback_delta, pairs).

Returns:

on success {"kg_used": True, "chosen_delta", "chosen_probe_pairs", "kg_scores": [{"delta", "probe_pairs", "score"[, "cost", "score_per_cost"]}, ...], "sigma", "cost_aware", "surrogate"} (sigma is None for a count posterior, which has no Gaussian noise scale); on fallback {"kg_used": False, "reason", "surrogate"}. surrogate is the compact surrogate_validity() report (ok/khat/params_flagged…) — when ok is False the Gaussian surrogate is suspect for this posterior and the wave’s scores should be spot-checked against (or replaced by) the NUTS-refit knowledge_gradient().

Return type:

(cells, meta) — the chosen design array and a JSON-safe meta dict

mmm_framework.continuous_learning.loop.run_closed_loop(world, *, center, B, value, n_geo=80, t_pre=6, t_test=10, delta=0.6, noise=0.6, n_holdout=0, mode='fixed', cap=None, gamma_scale=0.8, pair_signs=None, margin=1.0, population=1.0, wave_cost=0.05, max_waves=4, recenter=True, planner_q=120, fit_kwargs=None, use_laplace_kg=False, candidate_deltas=(0.3, 0.6, 0.9), kg_n_outcomes=32, stratify_geos=False, stop_patience=1, discount_half_life=None, spline_prior='iid', seed=0)

Run the full loop against a known world and return the decision trace.

Each pass: refit on all data, compute the recommendation + funding line + expected regret, evaluate the ENBS stopping rule, then (unless stopping) recenter on the recommendation, run a fresh designed wave over the SAME geos, and ingest it. Returns the per-wave WaveRecord history plus the final recommendation and the truth-optimal target for comparison.

use_laplace_kg=True (opt-in; the default path is byte-identical to the historical loop) scores candidate_deltas with select_next_design() each wave and runs the EVSI-best design instead of the fixed-delta CCD; the choice is recorded on the NEXT wave’s WaveRecord (kg_used / chosen_delta). stratify_geos=True (opt-in) blocks the geo→cell randomization on the true per-geo baselines a_geo (see assign_geos()).

stop_patience (default 1 — byte-identical to the historical loop) requires that many consecutive ENBS <= 0 evaluations before the loop stops. Evaluating ENBS every wave is NOT frequentist optional stopping — a Bayesian expected-value rule is stopping-rule-invariant under the likelihood principle (Edwards, Lindman & Savage 1963; Berger & Wolpert 1988) — but that invariance assumes a correctly specified model, and the misspecification study shows intervals can be narrow and wrong. stop_patience=2 is the cheap belt-and-suspenders guard: one optimistically-low regret read (from a temporarily overconfident, possibly misspecified posterior) cannot end the programme on its own. Each WaveRecord’s stop stays the raw per-wave verdict.

Return type:

dict[str, Any]

mmm_framework.continuous_learning.loop.due_for_retest(sigma_post, weeks_elapsed, channel, sigma_exp, *, half_life_overrides=None)

Whether a channel’s evidence has decayed enough to warrant a re-test.

Thin reuse of mmm_framework.planning.eig.reexperiment_due() so the continuous-learning loop’s re-test trigger uses the framework’s existing information-decay model (sigma_eff^2(t) = sigma_post^2 exp(lambda t), lambda = ln2 / half_life) rather than a private one. Returns (due, current_eig_nats).

Return type:

tuple[bool, float]

mmm_framework.continuous_learning.loop.adstock_half_life(adstock_alpha)

Half-life (periods) of a geometric adstock with retention alpha.

alpha**h = 1/2  =>  h = ln(1/2) / ln(alpha). The advertising literature finds carryover half-lives are strongly category-specific (industry FMCG figures cluster near 2–2.5 weeks against academic estimates of 7–12 — Fry, Broadbent & Dixon 2000), which is the empirical argument against one global decay constant: use each channel’s measured carryover as the anchor for its own evidence half-life, via due_for_retest(..., half_life_overrides={channel: h}). The carryover half-life is a lower anchor — a finding decays at least as fast as the mechanism that produced it drifts, and audience/creative drift usually puts the evidence half-life at a multiple of it; when the programme has accumulated waves, prefer the measured estimate_half_life().

alpha <= 0 (no carryover) returns 0.0; alpha >= 1 (no decay) returns inf.

Return type:

float

mmm_framework.continuous_learning.loop.estimate_half_life(gaps_weeks, shifts, sigmas)

Estimate the information-decay half-life from realized posterior drift.

The decay clock (math page Eq. 22) inflates a finding’s variance as sigma_eff^2(t) = sigma_post^2 * exp(lambda * t) with lambda = ln2 / h — but h need not be guessed: the loop observes how fast its own beliefs drift. Between two fits of the same quantity separated by a gap t_j (weeks), the decay model says the realized mean shift is a draw with variance sigma_j^2 * (exp(lambda t_j) - 1) on top of the posterior spread, so z_j = (shift_j / sigma_j)^2 has expectation exp(lambda t_j) - 1. The method-of-moments estimate solves mean_j[z_j] = mean_j[exp(lambda t_j) - 1] for lambda (monotone — a bisection). This is the adaptive-discounting idea of the dynamic-linear-model literature: West & Harrison (1997, §6.3) model each component as losing a constant information fraction delta per period (delta = exp(-lambda) here, returned as "discount"), with separate discounts per component — i.e. estimate one h per channel, not one global constant.

Parameters:

• gaps_weeks (ndarray) – per-observation elapsed weeks between the two fits.

• shifts (ndarray) – realized posterior-mean shifts of the tracked quantity (e.g. a channel’s marginal ROAS) over each gap.

• sigmas (ndarray) – the posterior sd of the quantity at the earlier fit.

Return type:

dict[str, float]

Returns:

{"half_life": h, "lambda": lam, "discount": exp(-lam), "n": J} — half_life = inf when no excess drift is observed (the findings are not decaying measurably). Requires at least one observation with positive gap and sigma.

Acquisition

Fast acquisition — Laplace knowledge-gradient and pure-EIG (guide §9.1/9.2).

The reference mmm_framework.continuous_learning.planner.knowledge_gradient() refits the model with full NUTS once per fantasy — too slow to score many candidate designs. This module replaces that with a Gaussian-linear treatment that needs no MCMC:

1. Moment-match the current posterior over the surface parameters to a Gaussian N(mu, Sigma) — in an unconstrained reparameterization eta (ThetaMap): positive parameters (beta, kappa, lam, …) are matched in log space, bounded ones (alpha, the mixture weight w) through a scaled logit, and sign-constrained synergies through a sign-aware log. Skewed posteriors of positive parameters are far closer to Gaussian on the log scale, and every fantasy sample maps back to a VALID parameter vector by construction (no clipping).

2. Linearize the surface in eta around mu (the Laplace step). A candidate design’s Fisher information is then Lambda = sum_cells w_c * u_c * g_c g_c^T with g_c = d incremental / d eta and u_c the observation family’s unit Fisher information at cell c (1/sigma^2 for the Gaussian; (nu+1)/((nu+3) sigma^2) for the Student-t; sigmoid(eta)^2 / (m + m^2/phi) with m = softplus(eta) for the NegBinomial count family — the GLM weight through the softplus link). The posterior covariance after running the design is Sigma_post = (Sigma^-1 + Lambda)^-1.

From this one linear-algebra object two acquisitions fall out cheaply:

• Laplace knowledge-gradient (decision value, EVSI): the pre-posterior spread of the updated mean is V = Sigma - Sigma_post; sample fantasy means eta_m ~ N(mu, V), map them back to valid surface parameters, re-optimize the allocation for each, and average the uplift over max_a profit(a | mu). Milliseconds instead of minutes.

• Pure EIG (information, D-/D_s-optimality): the entropy reduction of a parameter block S is 0.5 * (logdet Sigma_SS - logdet Sigma_post_SS). Use the full block for D-optimality, or the gamma sub-block (target= "gamma") for D_s-optimality — synergy-targeted waves the exploit-heavy Thompson waves under-probe. (EIG differences are invariant to the fixed reparameterization — the transform’s log-Jacobian cancels between the prior and posterior terms — so the D/D_s orderings match the constrained-space treatment.)

Any activation registered in surface.ACTIVATIONS with a transform spec in SHAPE_TRANSFORMS is supported (Hill, logistic, hill_mixture, monotone_spline), as is any observation family in model.VALID_LIKELIHOODS. The geo intercept is profiled out by centering the cell gradients across cells (the per-geo baselines are randomized, so the identifying variation is between-cell), mirroring planning.identification.

Prefer the Laplace KG when the goal is the decision; use EIG only to shore up decision-pivotal interactions (guide §9.2).

Two robustness layers guard the surrogate itself:

• Numerical PSD safeguards. V = Sigma - Sigma_post is PSD in exact arithmetic (Lambda >= 0) but finite-precision inversion/subtraction can leave tiny negative eigenvalues in V or a log-det sub-block. V is symmetrized and eigenvalue-clipped before fantasy sampling, and every log-determinant goes through the same symmetrize-and-clip guard — the nearest-PSD projection of Higham (1988; 2002).

• Surrogate validity (surrogate_validity()). The Gaussian moment-match is an approximation to the carried NUTS posterior; it degrades when the posterior is skewed, heavy-tailed or ridge-shaped (negative-gamma cannibalisation ridges are the expected failure mode — Kuss & Rasmussen 2005 for the intuition). The diagnostic scores the fit of N(mu, Sigma) to the draws in the unconstrained space — per-parameter skew/kurtosis flags plus a generalized-Pareto tail-shape khat of the Mahalanobis exceedances (the same Zhang & Stephens 2009 estimator and the same 0.7 alarm bar PSIS uses; Vehtari, Simpson, Gelman, Yao & Gabry 2024). When it fires, prefer the NUTS-refit planner.knowledge_gradient() for that wave; loop.select_next_design() records the report in its meta.

class mmm_framework.continuous_learning.acquisition.Transform(label, unconstrain, constrain)

Bases: NamedTuple

A support <-> R bijection: numpy unconstrain, JAX constrain.

label: str

Alias for field number 0

unconstrain: Callable[[ndarray], ndarray]

Alias for field number 1

constrain: Callable[[Any], Any]

Alias for field number 2

mmm_framework.continuous_learning.acquisition.SHAPE_TRANSFORMS: dict[str, dict[str, Transform]] = {'hill': {'alpha': ('logit[0.5,5]', <function _interval.<locals>.unconstrain>, <function _interval.<locals>.constrain>), 'kappa': ('log', <function <lambda>>, jax.numpy.exp)}, 'hill_mixture': {'alpha1': ('logit[0.5,6]', <function _interval.<locals>.unconstrain>, <function _interval.<locals>.constrain>), 'alpha2': ('logit[0.5,5]', <function _interval.<locals>.unconstrain>, <function _interval.<locals>.constrain>), 'kappa1': ('log', <function <lambda>>, jax.numpy.exp), 'kappa2': ('log', <function <lambda>>, jax.numpy.exp), 'w': ('logit[0,1]', <function _interval.<locals>.unconstrain>, <function _interval.<locals>.constrain>)}, 'logistic': {'lam': ('log', <function <lambda>>, jax.numpy.exp)}, 'monotone_spline': {'w1': ('log', <function <lambda>>, jax.numpy.exp), 'w2': ('log', <function <lambda>>, jax.numpy.exp), 'w3': ('log', <function <lambda>>, jax.numpy.exp), 'w4': ('log', <function <lambda>>, jax.numpy.exp), 'w5': ('log', <function <lambda>>, jax.numpy.exp), 'w6': ('log', <function <lambda>>, jax.numpy.exp), 'w7': ('log', <function <lambda>>, jax.numpy.exp), 'w8': ('log', <function <lambda>>, jax.numpy.exp), 'w9': ('log', <function <lambda>>, jax.numpy.exp)}}

Per-activation transform spec for the shape parameters, matching the prior supports in model._sample_activation_shape(). Register a new activation family here (plus in surface.ACTIVATIONS and the model’s shape sampler) and the whole acquisition layer picks it up.

class mmm_framework.continuous_learning.acquisition.ThetaMap(channels, pairs, activation, site_transforms, gamma_transforms)

Bases: object

Layout + bijection between posterior sites and the packed eta vector.

The flat layout is [beta (K), *shape sites (K each), gamma_pairs (P)] — the Hill case reproduces the historical [beta, kappa, alpha, gamma] packing order. All moment matching, linearization and fantasy sampling happen in the unconstrained eta space; constrain_params() maps any eta back to a VALID planner params dict.

channels: tuple[str, ...]

pairs: tuple[tuple[int, int], ...]

activation: str

site_transforms: tuple[tuple[str, Transform], ...]

gamma_transforms: tuple[Transform, ...]

property k: int

property dim: int

property gamma_idx: list[int]

Positions of the synergy parameters within the packed vector.

unconstrain_draws(samples)

Posterior draws -> the unconstrained (n_draws, dim) matrix.

Return type:

ndarray

constrain_jax(eta)

eta -> (beta, shape_tuple, gamma_matrix) in JAX ops.

constrain_params(eta)

eta -> a planner params dict (allocate_under_sample format).

The transforms enforce every support constraint, so the result is always a valid surface parameterization — no clipping needed.

Return type:

dict[str, Any]

surface_from_eta(eta, spend)

Incremental response at spend under packed eta (JAX).

__init__(channels, pairs, activation, site_transforms, gamma_transforms)

mmm_framework.continuous_learning.acquisition.theta_map(post)

Build the ThetaMap for a posterior’s activation + pair signs.

Return type:

ThetaMap

mmm_framework.continuous_learning.acquisition.eta_grad(eta_bar, spend_rows, tmap)

d incremental / d eta at eta_bar for each spend row -> (n, dim).

The chain rule through the constraining bijection is automatic — JAX differentiates surface(constrain(eta)) directly.

Return type:

ndarray

mmm_framework.continuous_learning.acquisition.theta_moments(post, *, ridge=1e-06, tmap=None)

Gaussian (mu, Sigma) matched to the posterior — in eta space.

The match happens in the unconstrained reparameterization of theta_map(): positive parameters in log space, bounded ones through a scaled logit, sign-constrained synergies through a sign-aware log. Any activation with a SHAPE_TRANSFORMS entry is supported; an unknown family raises NotImplementedError (the planner’s decision readouts stay activation-agnostic).

Return type:

tuple[ndarray, ndarray]

mmm_framework.continuous_learning.acquisition.surrogate_validity(post, *, tmap=None, tail_prob=0.25, khat_threshold=0.7, skew_threshold=0.5, kurt_threshold=1.0)

Diagnose whether N(mu, Sigma) is a valid stand-in for the posterior.

The Laplace acquisitions replace the carried NUTS posterior with its Gaussian moment-match in the unconstrained eta space (Eq. 12 of the math page). This runs entirely off the existing draws — no densities, no refit — and checks the two ways that stand-in fails: :rtype: dict[str, Any]

• Shape: per-parameter skewness and excess kurtosis of the unconstrained draws (a Gaussian has both ~0; the transforms already absorb the generic positivity skew, so what remains is real).

• Tails / ridges: the generalized-Pareto shape khat of the Mahalanobis-distance exceedances above their 1 - tail_prob quantile. Under a genuinely Gaussian posterior the squared distances are chi-square (exponential-class tail, khat ~ 0); a skewed, heavy-tailed or ridge-shaped posterior (the negative-gamma cannibalisation ridge is the expected offender) inflates khat. The alarm bar is the same 0.7 PSIS uses (Vehtari et al. 2024), estimated with the same Zhang–Stephens fit.

Returns a JSON-safe dict: ok (False when khat clears the threshold, any parameter is flagged, or there are too few draws to tell), khat, params_flagged, max_abs_skew, max_excess_kurtosis, and the per_param detail. When ok is False, spot-check or replace the wave’s Laplace scores with the NUTS-refit knowledge_gradient() (the correction path in the literature is a skewness-aware expansion — Rue, Martino & Chopin 2009).

mmm_framework.continuous_learning.acquisition.observation_unit_info(post, design_cells, tmap, mu_eta, *, sigma=None)

Per-cell Fisher information of ONE observation, for the fitted family. :rtype: ndarray

• "normal": 1 / sigma^2 (homoskedastic — constant across cells).

• "studentt": (nu + 1) / ((nu + 3) sigma^2) with the posterior-mean tail df — the classic heavy-tail efficiency discount (< the Gaussian info at any finite nu, recovering it as nu -> inf).

• "negbinomial": the GLM weight through the model’s softplus link, sigmoid(eta_c)^2 / (m_c + m_c^2 / phi) with eta_c = baseline + R(s_c; mu), m_c = softplus(eta_c) and the posterior-mean phi. The baseline is the posterior-mean A (or the grand mean of a_geo when the hyper is absent). sigma is ignored — a count family has no Gaussian noise scale.

sigma=None reads the posterior-mean sigma site for the Gaussian families; a missing required site raises ValueError (summaries-only fits never sample the observation sites, and no fixed guess is meaningful across KPI scales).

mmm_framework.continuous_learning.acquisition.design_information(design_cells, theta_bar, *, sigma=None, k=None, pairs=None, tmap=None, unit_info=None, cell_weights=None, n_geo=80, t_test=10, residualize=True)

Fisher information Lambda (dim×dim) a design carries about eta.

Lambda = sum_c w_c u_c (g_c - g_bar)(g_c - g_bar)^T where g_c is the unconstrained-space parameter gradient at cell c, u_c the observation family’s per-cell unit information (unit_info; defaults to the homoskedastic Gaussian 1/sigma^2), and g_bar the weighted mean across cells (profiling the geo intercept). residualize=False skips the centering (treats the baseline as known).

theta_bar is the unconstrained eta vector from theta_moments(). Pass the matching tmap; when omitted, a default Hill map with identity gamma transforms is built from (k, pairs) — only correct for a Hill posterior with default (“weak”/”zero”) pair signs.

Return type:

ndarray

mmm_framework.continuous_learning.acquisition.gaussian_eig(sigma0, lam, *, idx=None)

EIG (nats) of a parameter block under a Gaussian-linear update.

0.5 * (logdet Sigma0_SS - logdet Sigma_post_SS) for the block idx (all parameters if None) — the entropy the design removes. Both log-dets run through the _logdet_psd() symmetrize-and-clip guard: Sigma_post is PSD in exact arithmetic (Lambda >= 0) but the finite-precision double inversion can leave a marginal sub-block — equivalently the Schur complement Lambda_{S|rest} — with tiny negative eigenvalues.

Return type:

float

mmm_framework.continuous_learning.acquisition.design_eig(post, design_cells, *, sigma=None, target='all', n_geo=80, t_test=10, cell_weights=None)

Pure information gain of a design (D-optimal target="all" or D_s-optimal target="gamma" over the synergy sub-block).

Works for any registered activation and any fitted observation family; sigma=None derives the noise scale (Gaussian families) from the posterior. sigma is ignored for a count posterior — the GLM weights come from phi and the baseline instead.

Return type:

float

mmm_framework.continuous_learning.acquisition.laplace_knowledge_gradient(post, design_cells, B, value, *, sigma=None, n_geo=80, t_test=10, n_outcomes=64, mode='fixed', cap=None, n_starts=3, cell_weights=None, seed=0)

Decision-aware EVSI of a candidate design, with no refit (guide §9.1).

Gaussian-linear surrogate for mmm_framework.continuous_learning.planner.knowledge_gradient(): the pre-posterior spread of the updated mean is V = Sigma - Sigma_post; fantasy means eta_m ~ N(mu, V) are mapped back to valid surface parameters, re-optimized and averaged against the current best. Orders designs the same way as the NUTS KG at a fraction of the cost.

Works for any registered activation and fitted observation family (the Fisher weights come from observation_unit_info()); sigma=None derives the Gaussian-family noise scale from the posterior.

Each fantasy’s re-optimization is warm-started at the base (mu) allocation, in addition to the usual uniform/random multi-starts. Fantasy means are draws local to mu (V is a reduction in spread from the prior), so the true optimum is almost always near the base one; without the anchor, a non-concave surface (gamma can be negative) evaluated with only a handful of random restarts can converge to a different local optimum depending on the platform’s BLAS/LAPACK build (SLSQP’s floating- point path through the same starting point is not bitwise-portable), occasionally flipping a close two-candidate comparison. Anchoring removes that basin-hopping risk at no extra cost (it replaces one of the existing random starts, per _starts()).

Return type:

float

Dgp

Synthetic geo-RSM world with causal ground truth (recovery + fantasies).

A TrueWorld is a known response surface (beta, kappa, alpha, gamma) plus geo-intercept hyperparameters. simulate_panel() runs a central-composite wave against it and returns the data contract of guide §3 plus the answer key — this is BOTH the recovery harness (fit and compare to truth) and the engine that fantasizes wave outcomes inside the knowledge-gradient acquisition.

The DGP evaluates the same mmm_framework.continuous_learning.surface() functions the model fits, so a successful recovery test certifies the exact map the planner later exploits — there is no second, drifting implementation of the surface.

Identification (guide §3.2) is built in: a pre-period where every geo shares the center allocation (pins each geo intercept), then a test-period of designed cross-sectional variation (the CCD cells). Without the pre-period the geo intercept and the incremental response are collinear.

class mmm_framework.continuous_learning.dgp.TrueWorld(beta, gamma_pairs, channels, kappa=None, alpha=None, activation='hill', shape=<factory>, pairs=<factory>, a_level=4.0, sigma_a=1.0, phi_true=10.0, nu_true=4.0)

Bases: object

A known response surface for recovery / closure testing.

gamma_pairs is aligned to pairs (one synergy value per pair); the matrix form is assembled on demand. A/sigma_a set the geo-intercept distribution. The activation is pluggable: pass Hill kappa/alpha (the default and back-compatible path), or set activation + shape for another family (e.g. activation="logistic", shape={"lam": …}). phi_true is the NegativeBinomial concentration used when a simulation asks for noise_family="negbinomial"; nu_true is the Student-t tail df used by noise_family="studentt" (each ignored otherwise).

beta: ndarray

gamma_pairs: ndarray

channels: list[str]

kappa: ndarray | None = None

alpha: ndarray | None = None

activation: str = 'hill'

shape: dict[str, ndarray]

pairs: list[tuple[int, int]]

a_level: float = 4.0

sigma_a: float = 1.0

phi_true: float = 10.0

nu_true: float = 4.0

property n_channels: int

act_fn()

The JAX activation for this world’s family.

shape_tuple()

Shape-parameter arrays in the activation’s parameter order.

Return type:

tuple

gamma_matrix()

Return type:

ndarray

response_mean(spend_matrix)

Mean incremental response for a (N, K) spend panel (no noise).

Return type:

ndarray

answer_key()

Return type:

dict[str, object]

__init__(beta, gamma_pairs, channels, kappa=None, alpha=None, activation='hill', shape=<factory>, pairs=<factory>, a_level=4.0, sigma_a=1.0, phi_true=10.0, nu_true=4.0)

mmm_framework.continuous_learning.dgp.make_world(seed=0, channels=None)

A reproducible, non-trivial 4-channel world for tests and the demo.

Channels ["Chatter", "Pulse", "Orbit", "Vibe"] (generic social-network surfaces) with a mix of strong and weak main effects, distinct half- saturations and Hill shapes, and a synergy structure that includes a genuine cannibalization (negative gamma) and a complementarity (positive gamma) so recovery has to get signs right.

Return type:

TrueWorld

mmm_framework.continuous_learning.dgp.make_world_logistic(seed=0, channels=None)

A known logistic (exponential-saturation) world, sibling of make_world().

Same channels and synergy structure, but each channel saturates as f(s) = 1 - exp(-lam * s) — a concave curve with no S-shape — rather than a Hill. Demonstrates that the whole loop (fit, plan, stop, animate) is activation-agnostic. The saturation rates lam are chosen so half- saturation ln(2)/lam lands in the same O(1) scaled range as the Hill kappa above, keeping the two worlds comparable in scale.

Return type:

TrueWorld

mmm_framework.continuous_learning.dgp.make_world_hill_mixture(seed=0, channels=None)

A known weighted-sum-of-two-Hills world (a misspecification stress test).

Each channel’s true response is w·Hill(κ1,α1) + (1-w)·Hill(κ2,α2) with a low-κ, high-α early component (a soft activation threshold) plus a high-κ later component — a two-phase / shoulder shape a single Hill can only average over and a logistic (concave, no inflection) cannot represent at all. Fit it with activation="hill" (mild misspecification) or "logistic" (severe) to see what a wrong response family does to the recovered curve, the funding line, and the recommended allocation.

Return type:

TrueWorld

mmm_framework.continuous_learning.dgp.draw_geo_intercepts(world, n_geo, rng)

Sample the per-geo baseline intercepts a_geo ~ Normal(A, sigma_a).

Return type:

ndarray

mmm_framework.continuous_learning.dgp.drift_world(world, *, rate=0.08, drift_beta=True, drift_shape=True, drift_gamma=False, seed=0)

A slightly different TrueWorld — media behaviour drifting.

Applies one step of multiplicative log-normal jitter, x -> x * exp(Normal(0, rate)), to the channel ceilings beta and the activation shape params (bounded params are clipped back to their family’s support; drift_gamma=True also jitters the synergy magnitudes, sign-preserving). Applying it once per wave makes the response surface a geometric random walk — the changing-media- behaviour regime the information discount (fit(discount_half_life=...)) is built for. Channels, pairs and the geo-intercept hypers are untouched, so pair it with simulate_wave() over a fixed a_geo: the baselines stay pinned while the RESPONSE moves.

Return type:

TrueWorld

mmm_framework.continuous_learning.dgp.simulate_panel(world, center, *, n_geo=80, t_pre=6, t_test=10, delta=0.6, probe_pairs=None, noise=0.6, noise_family='normal', tau_scale=0.0, n_holdout=0, a_geo=None, adstock_alpha=None, adstock_l_max=8, stratify=False, seed=0)

Simulate a CCD wave against a known world.

Returns the data contract (spend, geo_idx, n_geo, y, period_idx) plus design, geo_alloc, cell_idx, a_geo, tau_true and answer_key. Pass an existing a_geo (and the same seed offset) to simulate a later wave over the same geos with their baselines intact.

period_idx is always provided (np.repeat(arange(t_pre + t_test), n_geo) — the row order is week-major, pre block then test block); a fit only uses it when it opts into time_effect="national".

noise_family="negbinomial" draws count outcomes (gamma–Poisson with mean softplus(mu) and concentration world.phi_true); the default "normal" path is byte-identical to the old Gaussian draw. tau_scale > 0 adds true national per-period shocks tau_t ~ Normal(0, tau_scale) to mu (0.0 leaves y byte-identical and draws nothing).

adstock_alpha adds carryover: the response is driven by the geometric-adstocked spend series (within each geo over weeks), while the returned spend stays raw. Fitting on the raw panel is then biased; the preprocess.adstock_prepass recovers it (guide §9.4).

stratify=True (opt-in — the default keeps the historical rng stream byte-identical) resolves a_geo BEFORE the geo assignment and blocks the randomization on it (assign_geos(..., baseline=a_geo)), matching the production practice of stratifying on the pre-period KPI level.

Return type:

dict[str, object]

mmm_framework.continuous_learning.dgp.simulate_wave(world, design, a_geo, *, t_test=10, center=None, n_holdout=0, noise=0.6, noise_family='normal', tau_scale=0.0, stratify=False, seed=0)

Simulate a single test window for a recentered design over the SAME geos.

Used by the closed loop to generate wave t > 0 outcomes: the geos and their baselines (a_geo) persist, only the test allocations change. period_idx (wave-LOCAL, np.repeat(arange(t_test), n_geo)) is always provided; ingest() applies the cross-wave offset. noise_family / tau_scale behave as in simulate_panel() (defaults byte-identical to the old draw). stratify=True (opt-in) blocks the geo→cell randomization on the known per-geo baselines a_geo.

Return type:

dict[str, object]