Build a Marketing Mix Model

Fit a Bayesian MMM to decompose sales and optimize ad spend across channels.

Introduction

30 Second Summary

Every company spends money on advertising across multiple channels. Few can confidently answer which ads actually drove sales and which were wasted budget.

In this project, you will build a Marketing Mix Model that measures how advertising spend drives weekly sales and recommends how to allocate budget across channels for maximum return.

What You'll Build

You will build an end-to-end Bayesian Marketing Mix Model using PyMC-Marketing that decomposes sales data into channel-level contributions and optimizes budget allocation.

By the end of this project, you'll have:

  • A fitted Bayesian MMM that decomposes weekly sales into channel, seasonality, and baseline contributions.
  • Channel-level ROAS estimates with uncertainty intervals showing which ads deliver the most revenue per dollar spent.
  • An optimized budget allocation recommendation that maximizes expected sales across marketing channels.
  • Secret Mission: Add channel-specific budget constraints and quantify their cost in expected sales.

Are there any prerequisites?

Basic Python and pandas familiarity is helpful. No Bayesian statistics background is required.

Set Up Your Python Environment

Before you can build a Marketing Mix Model, you need a Python environment with the right libraries installed.

You will use Google Colab as your notebook environment. Colab runs Python in the cloud so you do not need to install anything locally.

In this step, get ready to:
  • Open Google Colab and create a new notebook.
  • Install PyMC-Marketing and verify the installation.
  • Import all required libraries and configure plot settings.
Open Google Colab and create a new notebook

Google Colab is a free, browser-based notebook environment from Google. It lets you write and run Python code without installing anything on your machine.

  • Open your browser and navigate to colab.research.google.com.
  • Sign in with your Google account if prompted.
  • Click New Notebook in the bottom-right of the pop-up dialog.
  • Click the notebook title at the top left (it defaults to Untitled0.ipynb) and rename it to marketing_mix_model.ipynb.

What is Google Colab?

Colab is a hosted Jupyter notebook service. It gives you free access to Python with many data science libraries pre-installed.

Your notebooks are stored in Google Drive. You can share them just like a Google Doc.

Install PyMC-Marketing

The core library for this project is PyMC-Marketing. It provides pre-built Marketing Mix Model classes with adstock, saturation, and budget optimization built in.

  • In the first code cell of your notebook, install PyMC-Marketing by running this command:
!pip install pymc-marketing==0.19.4

What does this command do?

The ! prefix tells Colab to run a terminal command instead of Python code. pip install downloads the package from PyPI and installs it into your Colab environment.

The ==0.19.4 pins the exact version so your results match this project. PyMC-Marketing requires Python 3.12 or higher, which Colab already provides.

  • Press Shift+Enter to run the cell.
  • Wait for the installation to complete. You should see output scrolling through the installed dependencies, ending with a success message.

Installation taking a long time or failing?

The install can take 2-3 minutes because PyMC-Marketing has many dependencies. If you see a red error about version conflicts, try restarting your runtime by clicking Runtime in the menu bar, then Restart session. Then re-run the install cell.

help me fix my PyMC-Marketing installation error in Google Colab

Import libraries and configure settings

With PyMC-Marketing installed, you need to import the libraries you will use throughout the project. This includes tools for data manipulation, visualization, and Bayesian modeling.

  • Click the + Code button at the top left to create a new code cell below the install cell.
  • Paste the following imports into the new cell:
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

What are these libraries?

  • arviz provides tools for visualizing and diagnosing Bayesian model results.
  • matplotlib and seaborn handle plotting and data visualization.
  • numpy and pandas handle numerical arrays and tabular data.
  • Prior defines the Bayesian priors that shape the model's initial beliefs.
  • GeometricAdstock and LogisticSaturation model how ad effects carry over time and how returns diminish with higher spend.
  • MMM is the main Marketing Mix Model class that ties everything together.
  • Below the imports (in the same cell), add the plot configuration and random seed:
az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

What does this configuration do?

az.style.use("arviz-darkgrid") sets a clean visual theme for all Bayesian diagnostic plots.

The rcParams lines set a consistent figure size and resolution for all charts. The seed and rng variables create a reproducible random number generator so your model results match this walkthrough.

  • Press Shift+Enter to run the cell.
  • You should see no output if everything imported successfully. A blank output means all libraries loaded without errors.

Seeing an ImportError?

If you see an error like ModuleNotFoundError, the install cell may not have completed. Scroll back up and confirm the first cell finished running. Then click Runtime and select Restart session. Re-run both cells in order (install first, imports second).

help me fix my import errors for PyMC-Marketing in Google Colab

✔️ Awesome, I've got everything!

Great. Double check that you have saved your notebook. Colab auto-saves to Google Drive, but you can also press Ctrl+S (Windows) or Cmd+S (macOS) to save manually.

ⓧ I'd like to double check the full code

Here is everything your notebook should contain at this point across two cells:

Cell 1 (install):

!pip install pymc-marketing==0.19.4

Cell 2 (imports and config):

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

Your Python environment is ready. Next up, you will load real marketing data and explore how advertising spend and sales relate over time.

Load and Explore the Marketing Data

Your environment is set up and ready to go. Now you need data to model.

A Marketing Mix Model needs historical records of advertising spend alongside a business outcome like sales. In this step, you'll load a sample dataset from the PyMC-Marketing repository and explore what the data looks like before feeding it to the model.

In this step, get ready to:
  • Load a sample marketing dataset with 179 weeks of sales and channel spend.
  • Examine the data structure including columns, date range, and variable types.
  • Visualize the time series of sales and channel spend to spot patterns.
Load the sample marketing dataset

PyMC-Marketing provides a sample CSV file hosted on GitHub. It contains 179 weeks of simulated data with two advertising channels, control variables, and a target sales column.

  • In your Colab notebook, create a new code cell below your setup code and paste the following:
# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

What does this code do?

This cell loads the dataset directly from the PyMC-Marketing GitHub repository into a pandas DataFrame. The parse_dates argument converts the date_week column from text into proper datetime objects.

The print statements show the dataset's dimensions, column names, and the time period it covers. The data.head() call displays the first five rows so you can see real values.

  • Run the cell by pressing Shift+Enter.
  • You should see output showing Data shape: (179, 8), a list of column names, the date range, and a table preview of the first five rows.

What do the columns mean?

  • date_week is the weekly timestamp for each observation.
  • x1 and x2 are the weekly advertising spend for two marketing channels.
  • event_1 and event_2 are binary control variables representing special events (like promotions or holidays).
  • t is a trend variable capturing long-term growth or decline.
  • y is the target variable: weekly sales that the model will learn to decompose.

Seeing a connection error?

If you get a URL or connection error, check that your Colab notebook has internet access. Try clicking Runtime then Reconnect in the top menu and run the cell again.

If you see a NameError for pd, make sure you ran the Step 1 setup cell first.

help me debug a data loading error in my Colab notebook

Visualize sales and channel spend

Numbers alone don't tell you much about patterns in the data. A time series plot reveals seasonality, trends, and how spend varies week to week.

  • Create a new code cell and paste the following visualization code:
# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

What does this code do?

This creates a 3-panel figure with a shared x-axis (time). The top panel shows your target variable (weekly sales). The middle and bottom panels show spend for each advertising channel.

The sharex=True parameter aligns all three charts on the same time axis so you can visually compare spend spikes to sales peaks.

  • Run the cell by pressing Shift+Enter.
  • You should see a 3-panel chart. The top panel shows weekly sales fluctuating over time. The middle panel shows channel x1 with higher spend values. The bottom panel shows channel x2 with lower spend values.

Plot not showing?

If you see an empty output or an error, make sure you ran the data loading cell above first. The data variable must exist before you can plot it.

If the chart appears but looks squished, check that you included figsize=(12, 9) in the plt.subplots() call.

help me fix my matplotlib plot in Google Colab

✔️ Awesome, I've got everything!

Great. Double check that both cells ran without errors and that you can see the 3-panel chart.

ⓧ I'd like to double check the full code

Here is the complete notebook code up to this point. Compare your cells against this:

# Marketing Mix Model with PyMC-Marketing
# Run in Google Colab

# --- Step 1: Setup ---
# !pip install pymc-marketing==0.19.4

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

You now have a clear picture of what the data looks like. Next up, you'll configure the MMM to learn how each channel's spend translates into sales.

Build and Configure the MMM

You have data showing weekly sales and marketing spend across two channels. The next question is: how much of that sales activity can actually be attributed to each channel's advertising?

In this step, you will build a Bayesian Marketing Mix Model that accounts for two key advertising dynamics: carry-over effects (ads continue influencing sales for weeks after they run) and diminishing returns (doubling spend does not double results).

In this step, get ready to:
  • Set informative priors based on channel spend share.
  • Define the model with GeometricAdstock and LogisticSaturation.
  • Run a prior predictive check to verify the configuration is reasonable.
Set informative priors based on channel spend share

Before building the model, you need to separate your data into features (what the model sees) and a target variable (what the model predicts). You also need to tell the model how confident to be about each channel's contribution before it sees any data.

  • Create a new code cell and separate your features from the target variable by pasting this code:
# Prepare features and target
X = data.drop("y", axis=1)
y = data["y"]

print(f"Features: {X.columns.tolist()}")

What does this code do?

This separates your DataFrame into two pieces that the model needs.

  • X contains everything the model can use as input: channel spend, events, trend, and date.
  • y contains the weekly sales values the model will learn to predict.
  • Run the cell. You should see a printed list of feature column names including date_week, x1, x2, and the control/event columns.

Seeing a NameError?

If you see NameError: name 'data' is not defined, make sure you ran the cells from Step 2 that load the CSV file. Colab cells must be run in order.

help me fix a NameError for the data variable in my Colab notebook

Next, you will calculate how much each channel spends relative to the total. This spend share becomes the basis for your model's prior beliefs about channel importance.

  • Add a new code cell and calculate the spend share for prior specification by pasting this code:
# Calculate spend share for prior specification
total_spend_per_channel = data[["x1", "x2"]].sum(axis=0)
spend_share = total_spend_per_channel / total_spend_per_channel.sum()
print(f"\nSpend share per channel:\n{spend_share}")

n_channels = 2
prior_sigma = n_channels * spend_share.to_numpy()
print(f"\nPrior sigma for channels: {prior_sigma}")

Why use spend share for priors?

In Bayesian inference, priors encode your beliefs before seeing data. A channel that historically received more budget probably has a larger effect on sales.

By scaling the prior's width (prior_sigma) by spend share, you tell the model that higher-spend channels are allowed a wider range of possible effect sizes. This is an informative prior that guides the model toward sensible results.

  • Run the cell. You should see the spend share for each channel (they add up to 1.0) and the resulting prior_sigma array with two values.

Spend shares don't add up to 1.0?

Double-check that you are summing columns x1 and x2 with axis=0 (column-wise). Using axis=1 would sum across rows, giving you 179 values instead of 2.

help me debug my spend share calculation

Now you will assemble the full model configuration dictionary. This tells PyMC-Marketing what prior distributions to use for each model component.

  • Add a new code cell and define the model configuration by pasting this code:
# Define model configuration
my_model_config = {
    "intercept": Prior("Normal", mu=0.5, sigma=0.2),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=0.05, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=0.2, dims="fourier_mode"),
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=6), dims="date"),
}

my_sampler_config = {"progressbar": True}

What do these priors mean?

  • intercept: baseline sales when all marketing spend is zero. Centered at 0.5 (scaled units).
  • saturation_beta: how strongly each channel drives sales. Uses the spend-share-based sigma you just calculated.
  • gamma_control: effect of control variables (events, trend). Tight sigma of 0.05 keeps these small.
  • gamma_fourier: seasonality coefficients using a Laplace prior that encourages sparse (small) seasonal effects.
  • likelihood: the noise model for observed sales. A half-normal on sigma means the model expects moderate observation noise.
  • Run the cell. You won't see any output yet. The next chunk uses this configuration to build the model.
Define the model with GeometricAdstock and LogisticSaturation

The model configuration defines your priors. Now you will define the model architecture itself, choosing how advertising effects decay over time and how they saturate at high spend levels.

  • Add a new code cell and initialize the MMM instance by pasting this code:
# Initialize MMM
mmm = MMM(
    model_config=my_model_config,
    sampler_config=my_sampler_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=["x1", "x2"],
    control_columns=["event_1", "event_2", "t"],
    yearly_seasonality=2,
)

What are adstock and saturation?

GeometricAdstock models carry-over effects. When you run an ad this week, its impact doesn't vanish instantly. It decays geometrically over l_max=8 weeks. Think of it like a half-life for advertising memory.

LogisticSaturation models diminishing returns. The first thousand dollars of ad spend generates much more impact than the hundredth thousand. The logistic curve flattens as spend increases, capturing this real-world dynamic.

yearly_seasonality=2 adds 2 Fourier modes to capture repeating seasonal patterns in sales (like holiday spikes).

  • Run the cell. You won't see output yet. The next step builds the model graph and generates predictions from your priors.

Seeing an ImportError for MMM?

Make sure you imported from the correct path: from pymc_marketing.mmm.multidimensional import MMM. Older tutorials may reference a different import location that no longer works in version 0.19.4.

help me fix my MMM import error

Run a prior predictive check

A prior predictive check generates fake data from your priors before the model sees real data. If the simulated sales look wildly different from observed sales, your priors are misconfigured.

  • Add a new code cell and build the model, then generate prior predictive samples by pasting this code:
# Build model and run prior predictive check
mmm.build_model(X, y)
mmm.sample_prior_predictive(X, y, samples=1_000, random_seed=rng)

# Plot prior predictive
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.prior.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.prior.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Prior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

What does this plot show?

  • The black line is your actual observed sales data.
  • The blue line is the mean prediction from your priors (before the model has learned anything).
  • The shaded bands show 94% and 50% Highest Density Intervals (HDI). These represent the range of sales values your priors consider plausible.
  • If the observed data falls within the shaded bands, your priors are reasonable. The model can start from this belief and refine it with data.
  • Run the cell. This may take 10-20 seconds to sample. You should see a plot with a black line (observed sales) overlaid by a blue mean prediction line and shaded blue bands.
  • Confirm that the black observed line falls within the wider 94% HDI band. This means your priors are not ruling out the real data.

Observed data falls outside the bands?

If the black line sits far outside the shaded regions, the priors are too tight. Check that prior_sigma was computed correctly (it should be n_channels * spend_share.to_numpy(), not just spend_share alone).

Also verify you passed random_seed=rng to sample_prior_predictive. Without it your results will differ from the expected output, though they should still be reasonable.

help me debug my prior predictive check

✔️ Awesome, I've got everything!

Great. Double check you have all the cells from this step saved in your notebook.

ⓧ I'd like to double check the full code

Here is the complete notebook code from Step 1 through Step 3. Compare your cells against this:

# Marketing Mix Model with PyMC-Marketing
# Run in Google Colab

# --- Step 1: Setup ---
# !pip install pymc-marketing==0.19.4

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

# --- Step 3: Build and Configure the MMM ---
# Prepare features and target
X = data.drop("y", axis=1)
y = data["y"]

print(f"Features: {X.columns.tolist()}")

# Calculate spend share for prior specification
total_spend_per_channel = data[["x1", "x2"]].sum(axis=0)
spend_share = total_spend_per_channel / total_spend_per_channel.sum()
print(f"\nSpend share per channel:\n{spend_share}")

n_channels = 2
prior_sigma = n_channels * spend_share.to_numpy()
print(f"\nPrior sigma for channels: {prior_sigma}")

# Define model configuration
my_model_config = {
    "intercept": Prior("Normal", mu=0.5, sigma=0.2),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=0.05, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=0.2, dims="fourier_mode"),
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=6), dims="date"),
}

my_sampler_config = {"progressbar": True}

# Initialize MMM
mmm = MMM(
    model_config=my_model_config,
    sampler_config=my_sampler_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=["x1", "x2"],
    control_columns=["event_1", "event_2", "t"],
    yearly_seasonality=2,
)

# Build model and run prior predictive check
mmm.build_model(X, y)
mmm.sample_prior_predictive(X, y, samples=1_000, random_seed=rng)

# Plot prior predictive
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.prior.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.prior.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Prior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

Your model is configured and your priors pass the sanity check. Next up, you will fit this model to the real data using MCMC sampling and check whether it converges properly.

Fit the Model and Run Diagnostics

Your model is configured and the prior predictive check confirmed the setup is reasonable. Now it's time to fit the model to the actual data so it can learn the relationship between advertising spend and sales.

Fitting a Bayesian model means running MCMC sampling to explore the space of possible parameter values. After fitting, you need to verify the sampler converged properly and that the model actually captures the patterns in your data.

In this step, get ready to:
  • Fit the MMM using MCMC sampling with 4 chains.
  • Check for divergences and examine trace plots to confirm convergence.
  • Run a posterior predictive check to validate model fit against observed sales.
Fit the model using MCMC sampling

MCMC (Markov Chain Monte Carlo) is the algorithm that explores different parameter values and figures out which combinations best explain the observed sales data. Running multiple chains in parallel helps verify the sampler is finding consistent answers.

  • Create a new code cell in your notebook and fit the model by running this code:
# --- Step 4: Fit the Model and Run Diagnostics ---
_ = mmm.fit(
    X=X,
    y=y,
    chains=4,
    target_accept=0.85,
    random_seed=rng,
)

What does this code do?

This runs the MCMC sampler on your model. Here's what each parameter controls:

  • chains=4 runs 4 independent sampling chains. If they all converge to similar values, you can trust the results.
  • target_accept=0.85 controls how carefully the sampler explores. Higher values reduce divergences but slow things down.
  • random_seed=rng ensures reproducible results across runs.
  • Run the cell. This takes 2-5 minutes depending on your Colab runtime.

You should see a progress bar for each of the 4 chains. Once complete, the cell finishes without error.

Sampling taking too long?

If your Colab runtime is slow, make sure you are using a GPU runtime. Click Runtime in the top menu, then Change runtime type, and select a GPU option. You will need to rerun all cells from the top.

Help me debug slow MCMC sampling in my PyMC-Marketing model on Google Colab

Check for divergences and examine trace plots

Divergences are a warning sign that the sampler struggled to explore part of the parameter space. Zero divergences means the sampler ran smoothly and you can trust the parameter estimates.

  • Check for divergences by running this code in a new cell:
# Check divergences
n_divergences = mmm.idata["sample_stats"]["diverging"].sum().item()
print(f"Number of divergences: {n_divergences}")
if n_divergences == 0:
    print("No divergences - sampling was successful!")
else:
    print("Warning: Model had divergences. Consider increasing target_accept.")

What does this code do?

  • mmm.idata["sample_stats"]["diverging"] accesses the sampler diagnostics stored in the inference data object.
  • .sum().item() counts the total number of divergent transitions across all chains and converts it to a plain Python integer.
  • The conditional print tells you whether the sampling succeeded or if you need to adjust settings.
  • Run the cell. You should see Number of divergences: 0 followed by the success message.

Next, trace plots let you visually confirm that all 4 chains explored the same region of parameter space. Overlapping traces mean the chains agree on the answer.

  • Generate trace plots for the key model parameters by running this code in a new cell:
# Trace plots
_ = az.plot_trace(
    data=mmm.fit_result,
    var_names=["saturation_beta", "saturation_lam", "adstock_alpha"],
    compact=True,
    backend_kwargs={"figsize": (10, 6), "layout": "constrained"},
)
plt.gcf().suptitle("Trace Plots for Key Parameters", fontsize=16)
plt.show()

What does this code do?

  • az.plot_trace() from ArviZ creates two panels per parameter: the posterior distribution (left) and the sampling trace over time (right).
  • var_names selects the three most important parameters: channel effect sizes (saturation_beta), saturation rate (saturation_lam), and adstock decay (adstock_alpha).
  • compact=True overlays all chains on the same axes so you can see if they agree.
  • Run the cell. You should see trace plots with smooth, overlapping distributions on the left and well-mixed (caterpillar-like) traces on the right.

Seeing divergences or poorly mixed traces?

If you see divergences, try increasing target_accept to 0.90 or 0.95 in the mmm.fit() call and rerunning. If traces show chains stuck in different places (not overlapping), rerun the fit cell. Sometimes the sampler needs a second attempt.

Help me interpret my MCMC trace plots and fix convergence issues

Run a posterior predictive check

A posterior predictive check asks: "If the model's learned parameters are correct, can it reproduce the sales data it was trained on?" This is the ultimate test of whether your model captured the real patterns in the data.

  • Sample from the posterior predictive distribution and plot the results by running this code in a new cell:
# Posterior predictive check
mmm.sample_posterior_predictive(X, extend_idata=True, combined=True)

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.posterior_predictive.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.posterior_predictive.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Posterior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

What does this code do?

  • mmm.sample_posterior_predictive() generates predicted sales using the learned parameter values. extend_idata=True stores the predictions alongside the fitted parameters.
  • The black line shows actual observed sales. The blue line shows the model's mean prediction.
  • The shaded bands are HDI (Highest Density Intervals). The 94% band captures nearly all uncertainty, while the 50% band shows the most likely range.
  • Run the cell. You should see a plot where the blue prediction closely tracks the black observed line, with the shaded uncertainty bands wrapping tightly around the data.

If the black observed line sits within the shaded bands for most of the time series, your model is fitting well. The tighter the bands compared to the prior predictive check, the more the model learned from the data.

Posterior predictions don't match observed sales?

If the mean prediction line diverges significantly from the observed data, or the black line frequently falls outside the 94% HDI bands, the model may not have converged. Go back and check your trace plots for mixing issues. Also verify that the prior predictive check from Step 3 looked reasonable.

Help me understand why my posterior predictive check looks wrong

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is all the code from Step 1 through Step 4 combined. Compare this to your notebook cells:

# Marketing Mix Model with PyMC-Marketing
# Run in Google Colab

# --- Step 1: Setup ---
# !pip install pymc-marketing==0.19.4

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

# --- Step 3: Build and Configure the MMM ---
# Prepare features and target
X = data.drop("y", axis=1)
y = data["y"]

print(f"Features: {X.columns.tolist()}")

# Calculate spend share for prior specification
total_spend_per_channel = data[["x1", "x2"]].sum(axis=0)
spend_share = total_spend_per_channel / total_spend_per_channel.sum()
print(f"\nSpend share per channel:\n{spend_share}")

n_channels = 2
prior_sigma = n_channels * spend_share.to_numpy()
print(f"\nPrior sigma for channels: {prior_sigma}")

# Define model configuration
my_model_config = {
    "intercept": Prior("Normal", mu=0.5, sigma=0.2),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=0.05, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=0.2, dims="fourier_mode"),
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=6), dims="date"),
}

my_sampler_config = {"progressbar": True}

# Initialize MMM
mmm = MMM(
    model_config=my_model_config,
    sampler_config=my_sampler_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=["x1", "x2"],
    control_columns=["event_1", "event_2", "t"],
    yearly_seasonality=2,
)

# Build model and run prior predictive check
mmm.build_model(X, y)
mmm.sample_prior_predictive(X, y, samples=1_000, random_seed=rng)

# Plot prior predictive
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.prior.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.prior.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Prior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

# --- Step 4: Fit the Model and Run Diagnostics ---
_ = mmm.fit(
    X=X,
    y=y,
    chains=4,
    target_accept=0.85,
    random_seed=rng,
)

# Check divergences
n_divergences = mmm.idata["sample_stats"]["diverging"].sum().item()
print(f"Number of divergences: {n_divergences}")
if n_divergences == 0:
    print("No divergences - sampling was successful!")
else:
    print("Warning: Model had divergences. Consider increasing target_accept.")

# Trace plots
_ = az.plot_trace(
    data=mmm.fit_result,
    var_names=["saturation_beta", "saturation_lam", "adstock_alpha"],
    compact=True,
    backend_kwargs={"figsize": (10, 6), "layout": "constrained"},
)
plt.gcf().suptitle("Trace Plots for Key Parameters", fontsize=16)
plt.show()

# Posterior predictive check
mmm.sample_posterior_predictive(X, extend_idata=True, combined=True)

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.posterior_predictive.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.posterior_predictive.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Posterior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

Your model is fitted and validated. The diagnostics confirm it converged properly and can reproduce the observed sales patterns. Next up, you'll use this fitted model to decompose sales into channel contributions and estimate how much return each advertising channel delivers per dollar spent.

Analyze Channel Contributions and ROAS

Your model is fitted, validated, and ready to answer the question every marketing team asks: where is our money actually working?

In this step, you will decompose total sales into the contribution from each channel, seasonality, and baseline demand. You will also visualize how each channel's effectiveness diminishes as spend increases.

In this step, get ready to:
  • Generate a waterfall chart decomposing sales into components.
  • Plot each channel's contribution share with uncertainty.
  • Visualize saturation curves to see diminishing returns per channel.
Generate the waterfall decomposition chart

A fitted MMM knows how much each component contributed to total sales. The waterfall chart breaks the total into additive pieces: baseline demand, each channel's contribution, and seasonal effects.

  • In a new code cell, generate the waterfall decomposition by running:
# Waterfall decomposition
mmm.plot.waterfall_components_decomposition(original_scale=True)
plt.show()

What does this code do?

The waterfall_components_decomposition method sums up each component's contribution across all 179 weeks and displays them as a waterfall chart.

Setting original_scale=True means the y-axis shows actual sales units rather than the standardized scale the model uses internally.

You should see a waterfall chart with bars for the intercept (baseline), each channel (x1, x2), seasonality, and control variables, adding up to total predicted sales.

Seeing an empty or misshapen plot?

Make sure you ran the model fitting cell in Step 4 completely before this cell. The waterfall depends on posterior samples being present in mmm.idata.

If the plot renders but looks cut off, try adding plt.tight_layout() before plt.show().

help me debug my waterfall decomposition chart

Plot channel contribution share

The waterfall shows totals, but stakeholders want percentages. The contribution share plot shows what fraction of media-driven sales each channel is responsible for, complete with uncertainty intervals.

  • In a new code cell, plot the channel contribution share by running:
# Channel contribution share
fig = mmm.plot.channel_contribution_share_hdi(figsize=(8, 3))
plt.show()

What does this code do?

The channel_contribution_share_hdi method calculates each channel's percentage share of total channel contribution. The HDI (Highest Density Interval) bars show the uncertainty range around each estimate.

Because the model is Bayesian, you get a full distribution over contribution share rather than a single point estimate. The HDI bars communicate how confident the model is.

You should see a horizontal bar chart with one bar per channel (x1 and x2), each showing its percentage share of total channel-driven sales along with uncertainty bands.

Plot looks squished or labels overlap?

The figsize=(8, 3) parameter controls the plot dimensions. If labels are cut off, try increasing the second number (height) to 4 or 5.

help me fix my channel contribution share plot

Visualize saturation curves for each channel

Saturation curves answer a critical question: at what point does spending more on a channel stop producing proportional returns? A channel that has already hit its saturation point is wasting budget that could be reallocated elsewhere.

  • In a new code cell, visualize the saturation curves by running:
# Saturation curves (direct contribution)
fig = mmm.plot.saturation_scatterplot(original_scale=True)
plt.suptitle("Saturation Curves: Spend vs. Contribution", fontsize=16, y=1.02)
plt.show()

What does this code do?

The saturation_scatterplot method plots the learned saturation curve for each channel. Each point represents one week's spend (x-axis) vs. the contribution the model attributes to that spend (y-axis).

A curve that flattens early means the channel saturates quickly. Additional spend beyond the flattening point produces very little incremental contribution. This is the diminishing returns effect that the LogisticSaturation transformation captures.

You should see one subplot per channel, each showing a curved line (the saturation function) with scatter points (actual weekly observations) overlaid. A steeper, earlier-flattening curve indicates faster saturation.

Title overlapping with the subplot?

The y=1.02 parameter in plt.suptitle positions the title slightly above the figure. If it still overlaps, increase the value to 1.05.

help me fix my saturation curve plot layout

✔️ Awesome, I've got everything!

Great. Double check you have all three visualizations displayed: waterfall, contribution share, and saturation curves.

ⓧ I'd like to double check the full code

Here is the complete notebook code up to and including this step:

# Marketing Mix Model with PyMC-Marketing
# Run in Google Colab

# --- Step 1: Setup ---
# !pip install pymc-marketing==0.19.4

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

# --- Step 3: Build and Configure the MMM ---
# Prepare features and target
X = data.drop("y", axis=1)
y = data["y"]

print(f"Features: {X.columns.tolist()}")

# Calculate spend share for prior specification
total_spend_per_channel = data[["x1", "x2"]].sum(axis=0)
spend_share = total_spend_per_channel / total_spend_per_channel.sum()
print(f"\nSpend share per channel:\n{spend_share}")

n_channels = 2
prior_sigma = n_channels * spend_share.to_numpy()
print(f"\nPrior sigma for channels: {prior_sigma}")

# Define model configuration
my_model_config = {
    "intercept": Prior("Normal", mu=0.5, sigma=0.2),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=0.05, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=0.2, dims="fourier_mode"),
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=6), dims="date"),
}

my_sampler_config = {"progressbar": True}

# Initialize MMM
mmm = MMM(
    model_config=my_model_config,
    sampler_config=my_sampler_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=["x1", "x2"],
    control_columns=["event_1", "event_2", "t"],
    yearly_seasonality=2,
)

# Build model and run prior predictive check
mmm.build_model(X, y)
mmm.sample_prior_predictive(X, y, samples=1_000, random_seed=rng)

# Plot prior predictive
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.prior.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.prior.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Prior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

# --- Step 4: Fit the Model and Run Diagnostics ---
_ = mmm.fit(
    X=X,
    y=y,
    chains=4,
    target_accept=0.85,
    random_seed=rng,
)

# Check divergences
n_divergences = mmm.idata["sample_stats"]["diverging"].sum().item()
print(f"Number of divergences: {n_divergences}")
if n_divergences == 0:
    print("No divergences - sampling was successful!")
else:
    print("Warning: Model had divergences. Consider increasing target_accept.")

# Trace plots
_ = az.plot_trace(
    data=mmm.fit_result,
    var_names=["saturation_beta", "saturation_lam", "adstock_alpha"],
    compact=True,
    backend_kwargs={"figsize": (10, 6), "layout": "constrained"},
)
plt.gcf().suptitle("Trace Plots for Key Parameters", fontsize=16)
plt.show()

# Posterior predictive check
mmm.sample_posterior_predictive(X, extend_idata=True, combined=True)

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.posterior_predictive.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.posterior_predictive.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Posterior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

# --- Step 5: Analyze Channel Contributions and ROAS ---
# Waterfall decomposition
mmm.plot.waterfall_components_decomposition(original_scale=True)
plt.show()

# Channel contribution share
fig = mmm.plot.channel_contribution_share_hdi(figsize=(8, 3))
plt.show()

# Saturation curves (direct contribution)
fig = mmm.plot.saturation_scatterplot(original_scale=True)
plt.suptitle("Saturation Curves: Spend vs. Contribution", fontsize=16, y=1.02)
plt.show()

You now have a clear picture of how each channel drives sales and where diminishing returns kick in. Next up, you will use this information to optimize how budget is allocated across channels.

Optimize Budget Allocation

Your model now tells you how much each channel contributes to sales. The natural next question is: given a fixed total budget, how should you split spending across channels to maximize total sales?

In this step, you will use PyMC-Marketing's built-in budget optimizer to find the allocation that maximizes expected response. The optimizer uses SLSQP (Sequential Least Squares Quadratic Programming) under the hood to solve this constrained optimization problem.

In this step, get ready to:
  • Run the budget optimizer to find the optimal spend allocation across channels.
  • Compare the optimized allocation to a naive equal-split strategy.
  • Visualize the expected response distribution under the optimized budget.
Run the budget optimizer

The optimizer needs two inputs: the total budget to allocate and the number of time periods to optimize over. You will use the total historical spend as the budget so the comparison is apples-to-apples.

  • In your Colab notebook, create a new code cell and paste the following code to define the budget and run the optimizer:
# --- Step 6: Optimize Budget Allocation ---
# Define total budget (sum of historical spend)
total_budget = data[["x1", "x2"]].sum().sum()
num_periods = len(data)

print(f"Total historical budget: {total_budget:.2f}")
print(f"Number of periods: {num_periods}")

# Optimize budget allocation
allocation_strategy, optimization_result = mmm.optimize_budget(
    budget=total_budget,
    num_periods=num_periods,
)

print(f"\nOptimal allocation:\n{allocation_strategy}")
print(f"\nOptimization success: {optimization_result.success}")

What does this code do?

total_budget sums all historical spend across both channels. This gives you the same total dollars to redistribute.

mmm.optimize_budget() uses SLSQP to find the channel allocation that maximizes total expected sales, given the fitted model's learned saturation curves and adstock parameters. It returns an allocation_strategy (an xarray with optimal spend per channel) and an optimization_result (a scipy OptimizeResult confirming whether the solver converged).

  • Run the cell.
  • You should see the total historical budget printed, followed by the optimal allocation per channel and Optimization success: True.

Optimization success is False?

If the optimizer did not converge, make sure your model was fitted successfully in Step 4 (zero divergences). A poorly fitted model can produce an unstable optimization surface.

Help me debug why my budget optimization didn't converge in PyMC-Marketing.

Compare against a naive equal split

A naive strategy would split the total budget equally: 50% to channel x1 and 50% to channel x2. Look at your allocation_strategy output and compare.

How do I read the optimal allocation?

The allocation_strategy xarray shows the total spend the optimizer recommends for each channel across all periods. If the optimizer shifts budget toward one channel, it means that channel has more room on its saturation curve to generate incremental sales.

A naive 50/50 split ignores these saturation dynamics. The optimizer exploits the fact that each channel saturates at different rates to squeeze more total sales from the same budget.

Visualize the optimized response distribution

The optimizer found the best split, but how confident can you be in the expected outcome? Sampling from the model's posterior lets you see the full distribution of possible responses under the optimized budget.

  • Create a new code cell and paste the following to sample the response distribution and visualize the result:
# Sample response under optimized allocation
optimized_response = mmm.sample_response_distribution(
    allocation_strategy=allocation_strategy,
    time_granularity="weekly",
    num_periods=num_periods,
    noise_level=0.01,
)

# Visualize budget allocation
fig, ax = mmm.plot.budget_allocation(
    samples=optimized_response, figsize=(12, 8)
)
plt.show()

print("\n--- Summary ---")
print(f"Model fitted with {n_divergences} divergences")
print(f"Total budget: {total_budget:.2f}")
print(f"Optimized allocation: {allocation_strategy.to_numpy()}")
print("Budget optimization complete!")

What does this code do?

mmm.sample_response_distribution() draws samples from the model's posterior to simulate what sales would look like under the optimized allocation. The noise_level=0.01 adds a small amount of noise to make the simulation realistic.

mmm.plot.budget_allocation() produces a visualization showing the optimized spend per channel and the expected response distribution. The summary block at the end prints a clean recap of the full model run.

  • Run the cell.
  • You should see a budget allocation chart showing spend per channel and the expected response distribution, followed by a summary printout confirming the total budget and optimized allocation array.

Don't see the plot?

Make sure you ran the optimization cell first. The allocation_strategy and optimized_response variables must exist before calling the plot function.

Help me debug why my budget allocation plot isn't showing up.

✔️ Awesome, I've got everything!

Great. Double check you have both cells saved in your notebook.

ⓧ I'd like to double check the full code

Here is the complete notebook code up to and including this step:

# Marketing Mix Model with PyMC-Marketing
# Run in Google Colab

# --- Step 1: Setup ---
# !pip install pymc-marketing==0.19.4

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pymc_extras.prior import Prior

from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100

seed = sum(map(ord, "mmm"))
rng = np.random.default_rng(seed=seed)

# --- Step 2: Load and Explore Data ---
url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(url, parse_dates=["date_week"])

print(f"Data shape: {data.shape}")
print(f"Columns: {data.columns.tolist()}")
print(f"Date range: {data['date_week'].min()} to {data['date_week'].max()}")
data.head()

# Visualize sales and channel spend
fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)

axes[0].plot(data["date_week"], data["y"], color="black", linewidth=2)
axes[0].set(ylabel="Sales", title="Target Variable: Weekly Sales")

axes[1].plot(data["date_week"], data["x1"], color="C0", linewidth=2)
axes[1].set(ylabel="Spend", title="Channel x1 (Higher Spend)")

axes[2].plot(data["date_week"], data["x2"], color="C1", linewidth=2)
axes[2].set(xlabel="Date", ylabel="Spend", title="Channel x2 (Lower Spend)")

plt.tight_layout()
plt.show()

# --- Step 3: Build and Configure the MMM ---
# Prepare features and target
X = data.drop("y", axis=1)
y = data["y"]

print(f"Features: {X.columns.tolist()}")

# Calculate spend share for prior specification
total_spend_per_channel = data[["x1", "x2"]].sum(axis=0)
spend_share = total_spend_per_channel / total_spend_per_channel.sum()
print(f"\nSpend share per channel:\n{spend_share}")

n_channels = 2
prior_sigma = n_channels * spend_share.to_numpy()
print(f"\nPrior sigma for channels: {prior_sigma}")

# Define model configuration
my_model_config = {
    "intercept": Prior("Normal", mu=0.5, sigma=0.2),
    "saturation_beta": Prior("HalfNormal", sigma=prior_sigma, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=0.05, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=0.2, dims="fourier_mode"),
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=6), dims="date"),
}

my_sampler_config = {"progressbar": True}

# Initialize MMM
mmm = MMM(
    model_config=my_model_config,
    sampler_config=my_sampler_config,
    date_column="date_week",
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    channel_columns=["x1", "x2"],
    control_columns=["event_1", "event_2", "t"],
    yearly_seasonality=2,
)

# Build model and run prior predictive check
mmm.build_model(X, y)
mmm.sample_prior_predictive(X, y, samples=1_000, random_seed=rng)

# Plot prior predictive
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.prior.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.prior.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Prior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

# --- Step 4: Fit the Model and Run Diagnostics ---
_ = mmm.fit(
    X=X,
    y=y,
    chains=4,
    target_accept=0.85,
    random_seed=rng,
)

# Check divergences
n_divergences = mmm.idata["sample_stats"]["diverging"].sum().item()
print(f"Number of divergences: {n_divergences}")
if n_divergences == 0:
    print("No divergences - sampling was successful!")
else:
    print("Warning: Model had divergences. Consider increasing target_accept.")

# Trace plots
_ = az.plot_trace(
    data=mmm.fit_result,
    var_names=["saturation_beta", "saturation_lam", "adstock_alpha"],
    compact=True,
    backend_kwargs={"figsize": (10, 6), "layout": "constrained"},
)
plt.gcf().suptitle("Trace Plots for Key Parameters", fontsize=16)
plt.show()

# Posterior predictive check
mmm.sample_posterior_predictive(X, extend_idata=True, combined=True)

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(mmm.prior.date, mmm.y, color="black", label="Observed")
ax.plot(
    mmm.prior.date,
    mmm.posterior_predictive.y_original_scale.mean(("chain", "draw")),
    color="C0",
    label="Mean Prediction",
)
for hdi in (0.94, 0.5):
    az.plot_hdi(
        mmm.prior.date,
        mmm.posterior_predictive.y_original_scale,
        hdi_prob=hdi,
        smooth=False,
        ax=ax,
        color="C0",
        fill_kwargs=dict(alpha=0.3 if hdi == 0.94 else 0.5, label=f"{hdi:.0%} HDI"),
    )
ax.legend()
ax.set(title="Posterior Predictive Check", xlabel="Date", ylabel="Sales")
plt.show()

# --- Step 5: Analyze Channel Contributions and ROAS ---
# Waterfall decomposition
mmm.plot.waterfall_components_decomposition(original_scale=True)
plt.show()

# Channel contribution share
fig = mmm.plot.channel_contribution_share_hdi(figsize=(8, 3))
plt.show()

# Saturation curves (direct contribution)
fig = mmm.plot.saturation_scatterplot(original_scale=True)
plt.suptitle("Saturation Curves: Spend vs. Contribution", fontsize=16, y=1.02)
plt.show()

# --- Step 6: Optimize Budget Allocation ---
# Define total budget (sum of historical spend)
total_budget = data[["x1", "x2"]].sum().sum()
num_periods = len(data)

print(f"Total historical budget: {total_budget:.2f}")
print(f"Number of periods: {num_periods}")

# Optimize budget allocation
allocation_strategy, optimization_result = mmm.optimize_budget(
    budget=total_budget,
    num_periods=num_periods,
)

print(f"\nOptimal allocation:\n{allocation_strategy}")
print(f"\nOptimization success: {optimization_result.success}")

# Sample response under optimized allocation
optimized_response = mmm.sample_response_distribution(
    allocation_strategy=allocation_strategy,
    time_granularity="weekly",
    num_periods=num_periods,
    noise_level=0.01,
)

# Visualize budget allocation
fig, ax = mmm.plot.budget_allocation(
    samples=optimized_response, figsize=(12, 8)
)
plt.show()

print("\n--- Summary ---")
print(f"Model fitted with {n_divergences} divergences")
print(f"Total budget: {total_budget:.2f}")
print(f"Optimized allocation: {allocation_strategy.to_numpy()}")
print("Budget optimization complete!")

You now have a complete Marketing Mix Model that measures advertising effectiveness and recommends how to allocate budget. Next up in the Secret Mission, you will add real-world business constraints to this optimization and quantify their cost.

Secret mission

Add Business Constraints to Budget Optimization

In the real world, budget optimization is never unconstrained. Contracts lock in minimum spends, channel capacity limits cap how much you can push through a single platform, and leadership mandates minimum investment in brand channels regardless of short-term ROAS. Can you add per-channel constraints and quantify how much those constraints cost in expected sales?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely in Google Colab, so there are no ongoing costs. Decide whether to keep your notebook, disconnect the runtime to free up compute, or delete everything entirely.

Resources you used:

  • Google Colab notebook (stored in your Google Drive).
  • Google Colab runtime (ephemeral compute instance for running code).

Keep everything running

No action needed. Choose this if you plan to revisit the notebook, tweak priors, or experiment with different budget scenarios.

  • Your notebook stays in Google Drive at no cost.
  • The Colab runtime will eventually disconnect on its own after a period of inactivity.

Pause - I'll come back to this later

Disconnect the Colab runtime to free up compute resources. Your notebook remains saved in Google Drive for next time.

  • In your Colab notebook, click Runtime in the top menu bar.
  • Click Disconnect and delete runtime.
  • Your notebook code and outputs remain saved. The installed packages (like PyMC-Marketing) will need to be reinstalled when you reconnect.

Delete - I don't want to use this again

Remove all project resources completely. This deletes your notebook from Google Drive.

  • In your Colab notebook, click Runtime in the top menu bar, then click Disconnect and delete runtime.
  • Open Google Drive at drive.google.com.
  • Navigate to the Colab Notebooks folder.
  • Right-click your Marketing Mix Model notebook and select Move to trash.
  • Empty your trash to permanently delete the file.

That's a Wrap!

That's a Wrap!

Nice work! You just built the same type of model that marketing data scientists use at Fortune 500 retailers to drive millions in media optimization decisions.

You've learned how to:

  • Build a Bayesian Marketing Mix Model using PyMC-Marketing that decomposes sales into channel contributions, with full uncertainty quantification on ROAS estimates.
  • Interpret adstock (carry-over effects) and saturation (diminishing returns) curves to understand how advertising effectiveness changes with spend levels and time.
  • Run budget optimization using SLSQP to find the allocation that maximizes expected sales across marketing channels.
  • Add channel-specific budget constraints to the optimizer and quantify the cost of those constraints in expected sales.

Ready to quiz yourself?