The Mathematics of Digital Media Optimization: From Response Curves to Constrained Budget Allocation

A deep technical walkthrough of the models, math, and code that move you from “how much did each channel drive?” to “where should the next dollar go?” — with runnable implementations, a fully worked example, the Bayesian and causal layers, and reference patterns in Papilon and production Snowflake.

Why this is an optimization problem, not a reporting problem

Most media measurement stops at attribution: a dashboard reporting that paid search “got credit” for X conversions last month. That is backward-looking, and worse, it is usually wrong about the only thing that matters for planning — what happens when you move the next dollar.

The real question is forward-looking and constrained:

Given a fixed budget B, how should I split it across channels so that total predicted revenue is maximized — and would a larger B even pay for itself?

Framed that way, this is a constrained optimization problem sitting on top of a stack of statistical models. The interesting work lives in five layers:

Transforms — encoding carryover (adstock) and diminishing returns (saturation) so spend maps to response realistically.
Estimation — fitting a response curve per channel, ideally jointly and with honest uncertainty.
Optimization — solving for the allocation that maximizes total response under budget and business constraints.
Inference — making sure the curves reflect causal effects, not demand artifacts.
Validation — proving the model predicts out-of-sample and agrees with experiments.

This article works through each layer with the math and runnable numpy/scipy code, includes a fully worked synthetic example with real computed outputs, then maps everything onto the open-source landscape (Google Meridian, PyMC-Marketing, Meta Robyn) and two reference implementations — Papilon, the end-to-end causal MMM framework I maintain, and a lightweight Snowflake channel optimizer for production allocation.

The data foundation

Before any modeling, you need a modeling-ready table. For channel-level MMM that means a single tidy frame at a consistent grain — weekly is the standard, daily if you have the signal — with one row per period and aligned columns:

Outcome: revenue (or conversions, or whatever KPI you optimize)
Media: spend per channel (paid_search, paid_social, display, video, ...) — spend is fine; impressions or clicks are better if you have them, because they decouple cost inflation from delivery
Controls: seasonality flags, holidays, price index, promotions, macro covariates, competitor activity
Baseline drivers: anything that moves revenue independent of media (organic demand proxies, email sends)

The single biggest predictor of how long a model takes is data quality, not model choice. A clean weekly table with consistent channel taxonomy and aligned outcome columns gets you to a credible model in a day or two; messy inputs turn into 200 hours of debugging data problems disguised as model problems.

Layer 1: The two transforms that matter

Raw spend does not map linearly to revenue. Two non-linearities dominate, and every serious media model encodes both.

Carryover (adstock)

A dollar of advertising today keeps working tomorrow — brand impressions decay but do not vanish on impact. The standard model is geometric adstock, where effective media accumulates a decayed share of prior periods:

adstock(t) = spend(t) + λ · adstock(t−1),    0 ≤ λ < 1

λ is the carryover rate — the fraction of last period's effect surviving into this one. High-funnel channels (display, video, CTV) carry over heavily (λ ≈ 0.6–0.8); high-intent channels (branded search) decay almost instantly (λ ≈ 0–0.2).

import numpy as np

def geometric_adstock(spend, decay):
    """Media effect persists and decays geometrically over time."""
    out = np.zeros_like(spend, dtype=float)
    out[0] = spend[0]
    for t in range(1, len(spend)):
        out[t] = spend[t] + decay * out[t - 1]
    return out

One subtlety worth internalizing: under constant weekly spend s, geometric adstock converges to a steady-state level of s / (1 − λ). That multiplier matters when you optimize a recurring weekly budget rather than a one-time pulse — you are really optimizing the steady-state response, not the single-week response. (The worked example below uses exactly this correction.)

Delayed adstock. Geometric adstock peaks immediately and only decays. Many channels — especially upper-funnel video — peak after a lag. The delayed (or Weibull) form lets the peak land at lag θ:

w(l) = λ^((l − θ)²),    l = 0, 1, …, L
adstock(t) = Σ_l  w(l) · spend(t − l)  /  Σ_l w(l)

def delayed_adstock(spend, decay, theta, L=12):
    """Carryover with a delayed peak at lag theta (Weibull-style weights)."""
    lags = np.arange(L + 1)
    weights = decay ** ((lags - theta) ** 2)
    weights /= weights.sum()
    return np.convolve(spend, weights, mode="full")[: len(spend)]

Robyn supports both geometric and Weibull adstock; Meridian and PyMC-Marketing default to flexible geometric/delayed forms. Start geometric, reach for delayed only when a channel’s fit clearly demands a lagged peak.

Diminishing returns (saturation)

The second non-linearity is the one that drives allocation. Doubling paid search spend does not double conversions once you are already bidding on every relevant query — each marginal dollar buys less than the last. This is saturation, and the workhorse is the Hill curve (from biochemistry — receptor binding — a reminder that good math travels):

response(x) = Vmax · xⁿ / (Kⁿ + xⁿ)

Vmax — asymptotic ceiling (maximum achievable response from the channel)
K — half-saturation point (spend yielding half of Vmax)
n — shape/steepness; n = 1 gives the Michaelis–Menten curve, n > 1 gives an S-shape with a slow start

def hill(spend, vmax, k, n):
    """Hill saturation: response with diminishing returns. n=1 -> Michaelis-Menten."""
    spend = np.maximum(spend, 1e-9)
    return vmax * spend**n / (k**n + spend**n)

The common alternatives:

def logistic_saturation(x, lam):                 # Robyn / PyMC default flavor
    return (1 - np.exp(-lam * x)) / (1 + np.exp(-lam * x))

def log_saturation(x, beta):                     # simple, concave everywhere
    return beta * np.log1p(x)

Hill is the most flexible because n lets the data decide whether the curve is concave-from-zero (each dollar already past peak efficiency) or S-shaped (a threshold of spend before the channel “turns on”). The shape is the entire game for reallocation:

A channel below its inflection point has a steep curve — high marginal return, feed it. A channel past saturation has a flat curve — low marginal return, trim it.

Layer 2: Estimating the response curves

For intuition, fit each Hill curve directly from historical (adstocked spend, revenue contribution) with non-linear least squares:

from scipy.optimize import curve_fit

def fit_response_curve(adstocked_spend, contribution):
    """Fit a Hill saturation curve to one channel's adstocked history."""
    p0 = [contribution.max() * 1.4, np.median(adstocked_spend), 1.0]
    bounds = ([0, 0, 0.3], [np.inf, np.inf, 4.0])    # keep params identifiable
    params, _ = curve_fit(
        hill, adstocked_spend, contribution,
        p0=p0, bounds=bounds, maxfev=20000,
    )
    return params  # -> (vmax, k, n)

In production you never fit channels in isolation — you fit them jointly inside one regression so the baseline, seasonality, and controls are shared and the channels compete for explanatory credit. The full structural model:

revenue(t) = baseline(t)
           + Σ_k  β_k · saturationₖ( adstockₖ( spendₖ(t) ) )
           + Σ_j  γ_j · controlⱼ(t)
           + trend(t) + seasonality(t)
           + ε(t)

This is exactly what Meridian, PyMC-Marketing, and Robyn implement. The per-channel fit above is the conceptual building block — the joint model is the production object, because only the joint fit can correctly partition revenue across correlated channels.

Layer 3: The optimization

With a fitted response curve Rₖ(·) per channel, allocation becomes a clean constrained program:

maximize    Σ_k  R_k(s_k)
   s
subject to  Σ_k s_k = B                 (spend the whole budget)
            s_k^min ≤ s_k ≤ s_k^max     (per-channel floors / caps)

Because Hill curves are smooth and the constraints are linear, Sequential Least Squares Programming (SLSQP) is the natural solver — it handles equality and inequality constraints with bounds and converges fast at this scale.

from scipy.optimize import minimize

def optimize_budget(curves, decays, total_budget, bounds_per_channel):
    """
    curves: {channel: (vmax, k, n)}; decays: {channel: lambda}
    Optimizes STEADY-STATE response for a recurring weekly budget.
    """
    channels = list(curves.keys())

    def neg_total_revenue(s):
        total = 0.0
        for i, ch in enumerate(channels):
            ss_spend = s[i] / (1 - decays[ch])          # steady-state adstock
            total += hill(ss_spend, *curves[ch])
        return -total

    constraints = [{"type": "eq", "fun": lambda s: np.sum(s) - total_budget}]
    bounds = [bounds_per_channel[ch] for ch in channels]
    x0 = np.full(len(channels), total_budget / len(channels))

    result = minimize(neg_total_revenue, x0, method="SLSQP",
                      bounds=bounds, constraints=constraints)
    return dict(zip(channels, result.x)), -result.fun

The economic intuition: equal marginal returns

Understand why the optimizer lands where it does, because it gives you a sanity check that does not require trusting the solver. Form the Lagrangian for the equality-constrained problem:

L = Σ_k R_k(s_k) − μ ( Σ_k s_k − B )

Setting ∂L/∂s_k = 0 gives the KKT optimality condition for every interior channel:

dR_k/ds_k = μ     for all k not pinned at a bound

In plain language: at the optimum, the marginal return of the last dollar is equal across every channel. If branded search returns $0.40 on its next dollar while display returns $0.90, you are not optimal — money should flow from search to display until their marginal returns meet. The multiplier μ is the portfolio marginal ROAS at budget B, which also answers "should I raise the total budget?" (raise it while μ > 1).

def marginal_roas(spend, vmax, k, n, decay, eps=1.0):
    f = lambda x: hill(x / (1 - decay), vmax, k, n)
    return (f(spend + eps) - f(spend)) / eps

A correct solution has near-equal marginal_roas across all interior channels — the single most useful diagnostic you can run.

Beyond a single objective

Two extensions matter in practice.

Risk-adjusted allocation. With a Bayesian model you have a posterior over response, not a point estimate. Optimize a risk-penalized objective using posterior draws, trading expected revenue against its variance (a mean–variance / Markowitz framing borrowed from portfolio theory):

maximize   E[ R(s) ] − γ · SD[ R(s) ]

where γ is your risk aversion. This naturally down-weights channels whose effects are large in expectation but highly uncertain.

Flighting (allocation across time). Expand the decision from s_k to s_{k,t} — spend per channel and week — subject to a total over the horizon. Because adstock couples periods, the optimizer will discover pulsing strategies (concentrate spend) versus even spreading depending on each channel's decay and saturation. This is a larger but structurally identical SLSQP problem.

A fully worked example (real computed outputs)

To make it concrete, here is a synthetic but realistic run: 104 weeks, three channels with known adstock and Hill parameters, fit and then optimized at a fixed weekly budget. All numbers below are actual program output.

The recovered curves (fit on adstocked spend):

paid_search   decay=0.20   Vmax=90,000    K=18,000   n=1.30
paid_social   decay=0.50   Vmax=120,000   K=40,000   n=1.60
display       decay=0.70   Vmax=70,000    K=55,000   n=1.20

The current allocation (recent 13-week average), total weekly budget $72,314:

paid_search   $20,292
paid_social   $24,977
display       $27,045

Running SLSQP with ±guardrails (no channel moves more than its bounds allow) on the same total budget:

OPTIMIZED (same $72,314 total):
  paid_search   $22,050   (delta  +1,758)
  paid_social   $33,218   (delta  +8,241)
  display       $17,046   (delta −10,000)

Predicted weekly revenue:  current $170,508  ->  optimized $175,936   (+3.2%)

The optimizer pulls $10K out of display — high carryover but low ceiling and already past its useful spend — and pushes it into paid social, which has the highest ceiling and is sitting below saturation. Same budget, +3.2% revenue, purely from reallocation.

And the proof it found a real optimum — marginal ROAS at the solution:

paid_search   mROAS = 1.230
paid_social   mROAS = 1.231
display       mROAS = 1.231

The marginal returns equalize at ≈1.23 across all three interior channels — exactly the KKT condition dR_k/ds_k = μ. That equalization is the optimum; if these numbers diverged, money would still be misallocated.

Layer 4: Why you go Bayesian

The least-squares fit above is a fine teaching tool and a workable v1. It has two real weaknesses the modern toolset is built to fix.

Honest uncertainty. A point estimate of Vmax or channel ROI hides how much you actually know. With correlated channels and ~100 weekly observations, confidence in any single contribution can be wide — and a recommendation that ignores that is overconfident. Bayesian MMM produces a full posterior distribution per parameter, so every contribution and response curve carries credible intervals, and you can encode platform benchmarks or prior lift-test results as priors instead of pretending the model knows nothing.

A representative Bayesian specification:

revenue(t) ~ Normal( μ(t), σ )
μ(t)       = α + trend(t)
              + Σ_k β_k · saturationₖ( adstockₖ( xₖ(t) ; λ_k ) )
              + Σ_j γ_j · controlⱼ(t)

β_k        ~ HalfNormal(·)            # channel effects constrained positive
λ_k        ~ Beta(·, ·)               # adstock decay in [0,1)
K_k, n_k   ~ Gamma/HalfNormal(·)      # saturation shape priors
σ          ~ HalfNormal(·)

The posterior is estimated via MCMC (Hamiltonian Monte Carlo / NUTS) or variational inference. In PyMC-Marketing that is a few lines:

from pymc_marketing.mmm import MMM, GeometricAdstock, LogisticSaturation

mmm = MMM(
    date_column="date",
    channel_columns=["paid_search", "paid_social", "display"],
    control_columns=["holiday", "price_index"],
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
)
mmm.fit(X, y)                          # NUTS posterior (NumPyro / Nutpie backend)
mmm.sample_posterior_predictive(X)
allocation = mmm.optimize_budget(budget=1_000_000, num_periods=13)

Google Meridian follows the same Bayesian-causal philosophy with a geo-hierarchical model and TensorFlow-Probability sampling; the workflow is build InputData → ModelSpec (with priors) → sample_posterior() → BudgetOptimizer. Its distinctive feature is shipping interactive scenario-planning dashboards (Looker Studio) that non-technical stakeholders use for planning without touching code. (Check the current Meridian docs for exact API — it requires Python 3.11–3.13 and a GPU is recommended.)

The current open-source landscape:

Capability Google Meridian PyMC-Marketing Meta Robyn Papilon Estimation Bayesian MCMC (TF-Probability) Bayesian (NumPyro / BlackJAX / Nutpie) Ridge regression + Nevergrad evolutionary search Bayesian hierarchical Language Python Python R Python Adstock + saturation Yes Yes Yes Yes Time-varying coefficients Yes Yes Limited Yes Experiment calibration (lift) Yes (geo) Yes Limited Yes Causal discovery / inference — — — Yes Agent-based / complex-systems sim — — — Yes Stakeholder dashboards Yes (Looker) Build-your-own Static PNGs Build-your-own Best for Google-heavy spend, geo data, exec dashboards Maximum flexibility, custom priors R teams, strong docs/community Causal + simulation under one roof

Practical notes from current releases: Google deprecated LightweightMMM in early 2025 and Meridian (v1.7 as of 2026) is its supported successor. PyMC-Marketing is the most-downloaded Python option and the most flexible, with swappable sampling backends and benchmark-leading sampling speed. Robyn remains the popular R choice but uses ridge regression with evolutionary hyperparameter search rather than a full posterior — point estimates, not credible intervals.

Layer 5: The causal layer (the part most models skip)

Regression conflates the conditional expectation E[revenue | spend] with the interventional quantity you actually want, E[revenue | do(spend)). The gap is not academic. When spend responds to demand — you bid harder during high season, you retarget users already intending to buy — a naive model attributes the demand's lift to the spend, inflating ROI on exactly the channels you most want to judge honestly.

Concretely, suppose demand Z drives both spend X and revenue Y (X ← Z → Y). Then Z is a confounder, and the backdoor path biases the naive estimate upward. The fix is to identify and adjust for the confounder (or block the backdoor path), formalized by Pearl's backdoor criterion:

P(Y | do(X)) = Σ_z  P(Y | X, Z = z) · P(Z = z)

There are three practical routes:

Encode structure. Specify a DAG, identify confounders, and adjust for them — the explicit causal-inference path.
Calibrate against experiments. Use geo lift tests or holdout experiments as ground-truth incrementality, and fit the MMM so its channel contributions match the experimental estimates. Meridian and PyMC-Marketing both support this; it is the single most powerful thing you can do to make an MMM trustworthy.
Triangulate. Reconcile MMM (top-down, causal-ish), multi-touch attribution (bottom-up, correlational), and incrementality experiments (gold standard but narrow). Where they disagree, dig in.

This causal layer is what Papilon was built around — putting Bayesian MMM and explicit causal inference under one API, because in practice you want both the response curves and a defensible answer to “is this effect real or just confounded by demand?”

import papilon as pp

# Recover a plausible causal graph from the data
dag = pp.discover_causal_structure(df)

# Estimate an interventional effect, not just a correlation
effect = pp.estimate_effect(
    df,
    treatment="paid_search_spend",
    outcome="revenue",
    # adjustment set inferred from the discovered DAG
)

This matters most for the channels entangled with demand — branded search and retargeting are the usual offenders, where a model wants to over-credit spend that was really just harvesting existing intent.

Validation: proving the model before you trust the budget

An MMM that fits the training data perfectly and predicts nothing out-of-sample is worse than no model. Minimum bar before any recommendation ships:

Time-series cross-validation (rolling origin). Never random K-fold — that leaks the future into the past. Train on weeks 1–t, predict t+1…t+h, roll forward, average the error. Report out-of-sample NRMSE / MAPE.
Holdout. Reserve the most recent 8–13 weeks, fit on the rest, score the holdout.
Posterior predictive checks (Bayesian). Simulate from the posterior and confirm the replicated series looks like the observed one. Watch R-hat (≈1.0) and divergences for sampling health.
Calibration against experiments. The decisive test — do the model’s channel contributions agree with geo lift results? If the MMM says paid social ROI is 4× but a lift test says 1.5×, trust the experiment and re-fit with it as a prior.
Face validity. Decay rates, saturation shapes, and ROIs should be defensible to a media planner. A branded-search adstock of 0.9 or a display ROI above paid search should trigger scrutiny, not celebration.

def rolling_origin_cv(df, fit_fn, predict_fn, min_train=52, horizon=4):
    errors = []
    for t in range(min_train, len(df) - horizon):
        model = fit_fn(df.iloc[:t])
        pred  = predict_fn(model, df.iloc[t : t + horizon])
        actual = df["revenue"].iloc[t : t + horizon].values
        nrmse = np.sqrt(np.mean((pred - actual) ** 2)) / actual.mean()
        errors.append(nrmse)
    return np.mean(errors), np.std(errors)

The portfolio question: is B even the right budget?

Channel split assumes the total is fixed. It usually should not be. Sweep the budget B, re-optimize at each level, and plot predicted revenue versus spend — the portfolio response curve. Its slope at any point is the portfolio marginal ROAS (the shadow price μ from earlier):

def efficient_frontier(curves, decays, bounds_fn, lo, hi, steps=25):
    frontier = []
    for B in np.linspace(lo, hi, steps):
        alloc, rev = optimize_budget(curves, decays, B, bounds_fn(B))
        frontier.append((B, rev))
    B, R = np.array(frontier).T
    marginal = np.gradient(R, B)        # dRevenue/dBudget at each total
    return B, R, marginal

Keep raising the budget while the portfolio marginal ROAS exceeds your target (≥ 1.0 for revenue-neutral, higher for a margin hurdle). The point where the marginal curve crosses your threshold is the defensible total budget — a far stronger answer than “spend what we spent last year.”

Reference implementation: Papilon

Papilon — the framework I maintain under Vector1 Research — puts the Bayesian and causal layers under one API. The MMM path mirrors the math above but fits jointly with Bayesian inference and runs constrained optimization on the posterior curves:

import papilon as pp
from papilon.mmm import MMM

model = MMM(
    data=df,
    target="revenue",
    media=["paid_search", "paid_social", "display", "video"],
    date_col="date",
)
model.fit()                                   # Bayesian posteriors per channel
rois = model.get_roi()                        # ROI with credible intervals
allocation = model.optimize_budget(total=1_000_000)

Where it diverges from a pure MMM library is the surrounding stack — causal discovery and inference (above), Monte Carlo / bootstrap simulation for scenario uncertainty, agent-based and economic input-output modeling for second-order effects, and an orchestration layer for scheduled retraining.

Productionizing: the Snowflake allocation pattern

A trained model is only useful if the allocation refreshes against live performance data. The lightweight pattern I run inside a Snowflake Python notebook keeps the whole loop next to the ROAS and conversion data where it already lives:

Pull channel history (spend, revenue, conversions, roas, ctr) at weekly grain from a clean modeling table.
Transform each channel’s spend through adstock.
Fit a response curve per channel (or load posteriors from the trained MMM).
Optimize next period’s allocation with SLSQP under the live budget and per-channel guardrails.
Write a recommendations table back to the warehouse — recommended spend, predicted revenue, marginal ROAS, delta vs. current — so it lands where planners already work (Snowflake → Tableau / Looker).

# Inside a Snowflake Python notebook
allocation, predicted_revenue = optimize_budget(
    curves=fitted_curves,                 # {channel: (vmax, k, n)}
    decays=fitted_decays,                 # {channel: lambda}
    total_budget=current_quarter_budget,
    bounds_per_channel=channel_guardrails,
)

recommendations = (
    pd.DataFrame({"channel": list(allocation),
                  "recommended_spend": list(allocation.values())})
    .assign(
        current_spend = lambda d: d.channel.map(current_spend),
        delta         = lambda d: d.recommended_spend - d.current_spend,
        marginal_roas = lambda d: d.apply(
            lambda r: marginal_roas(r.recommended_spend,
                                    *fitted_curves[r.channel],
                                    fitted_decays[r.channel]), axis=1),
    )
)
# session.write_pandas(recommendations, "MEDIA_ALLOCATION_RECS")

The guardrails are where the model meets reality: you cannot 10× paid social in a week even if the math says so, and some channels carry contractual minimums. Encoding those as bounds keeps recommendations actionable instead of theoretically optimal. Retrain on a schedule (weekly is typical) — a response curve fit on six-month-old data quietly drifts away from current auction economics and creative fatigue.

Putting it together

The throughline: media optimization is curve-fitting plus constrained optimization, and the answer is only as good as the curves. Get adstock and saturation right, fit them jointly with honest uncertainty, confirm the effects are causal rather than confounded, validate out-of-sample and against experiments — and the allocation step is almost mechanical, SLSQP equalizing marginal returns across channels under your budget. In the worked example that discipline turned a pure reallocation into a +3.2% revenue gain at zero incremental spend.

Where frameworks differ is how much rigor they bring to the curves and the inference. A curve_fit + minimize script gets you a defensible v1 in an afternoon. Meridian, PyMC-Marketing, and Papilon earn their complexity when the budgets are large enough that overconfident point estimates cost real money — and when "is this effect causal?" stops being academic and starts deciding where millions of dollars go.

If there is one place to invest beyond the basics, it is the causal and experimental-calibration layer. The next several years of media measurement are not about better attribution — they are about better counterfactuals.

All code is illustrative and runnable as written with numpy and scipy; the worked-example figures are actual program output. Framework snippets reflect each library's public interface — consult current documentation for versions and installation (Meridian requires Python 3.11–3.13 and a GPU is recommended). Replace the Papilon GitHub placeholder link with the live repository URL before publishing.

Brian Curry is a Kansas City–based cognitive systems architect and applied AI researcher working at the intersection of AI architecture, knowledge engineering, and intelligent-system design. He is the founder of Vector1 Research, an applied AI lab publishing on agent systems, AI evaluation, cognitive architecture, and the macroeconomics of AI-driven transformation. He has led applied AI and data initiatives across major enterprises including Koch Industries, Tractor Supply, Vail Resorts, Hallmark, Garmin, AT&T, and McClatchy, and founded KC AI Lab, a 1,500-member AI/ML community.

His open-source work includes Papilon (causal inference and counterfactual simulation), PyCausalSim (causal discovery via agent-based simulation), MeaningFlow (semantic modeling and ontology), and Memory-Node Encapsulation (episodic memory for agentic systems). His current research focuses on agent evaluation, knowledge engineering, and AI architecture in regulated, operationally-critical domains.

Contact: brian at vector1.ai · LinkedIn