MMM（マーケティングミックスモデル）を個人で実装する完全手順（Python）
メディア出稿・売上・外的変数からMMMをゼロから実装する手順を、PyMCの動くコードで通します。アドストックと飽和の自作、ベイズ推定による寄与度分解、予算配分シミュレーションまで。OSSライブラリの比較ではなく「自分で組む」ことに絞った実装記事です。
by Shin
MMMを個人で実装するには何が必要なのか。答えは「アドストック・飽和の2つの変換と、ベイズ回帰、それだけ」です。専用ライブラリを使わずとも、出稿データと売上があれば自分で組めます。この記事では合成データで全コードを動かしながら、前処理→推定→寄与度分解→予算配分シミュレーションまでを一気通貫で通します。
対象は /lab レイヤー、つまり手を動かす実装者です。MMMを「なぜ経営に使うのか」「何が測れて何が測れないのか」という意思決定側の話は なぜMMMが経営に効くのか に分けてあります。ここでは判断の話はせず、コードだけで進めます。
検証環境は Python 3 系・PyMC 5 系以降（手元では 6 系でも動作確認済み）・scipy。本記事のコードはすべて手元で実行し、合成データでのパラメータ回復を確認しています。
データ準備（メディア出稿・売上・外的変数）
MMMに必要な最小構成は3種類です。第一に目的変数（週次の売上やCV）、第二に各チャネルの出稿量（金額でもインプレッションでもよいが単位は固定）、第三に外的変数（トレンド・季節性・プロモ・祝日・天候など、売上を動かすが広告ではない要因）。粒度は週次が定番です。日次はノイズが多く、月次は観測点が足りません。
外的変数を入れない MMM は、広告に「広告以外の効果」まで背負わせてしまい、寄与度を過大評価します。ここを削るとモデルが嘘をつくので、最初から枠を作っておきます。
説明のため、真のパラメータを既知にした合成データを作ります。手元の実装が「正しく真値を取り戻せるか（パラメータ回復）」を確認できるので、自作MMMの検証では最初にやるべき手順です。
import numpy as np
import pandas as pd
rng = np.random.default_rng(42)
N = 130 # 週数（約2.5年）
dates = pd.date_range("2024-01-07", periods=N, freq="W-SUN")
t = np.arange(N)
# メディア出稿（3チャネル）。social/displayは断続的なパターンを混ぜる
search = np.abs(rng.normal(100, 30, N)) * (0.5 + 0.5 * np.sin(t / 6))
social = np.abs(rng.normal(60, 25, N))
display = np.abs(rng.normal(40, 20, N)) * (rng.random(N) > 0.3)
X_raw = np.vstack([search, social, display]).T # (N, 3)
channels = ["search", "social", "display"]
# 外的変数：トレンド・季節性・プロモ
trend = 0.4 * t
seasonal = 8 * np.sin(2 * np.pi * t / 52)
promo = (rng.random(N) > 0.85).astype(float)季節性はここでは単純な正弦波で代用していますが、実データでは年次のフーリエ項（sin/cos を数本）かダミー変数で表現します。週番号 0..51 を周期 52 のフーリエ基底に展開するのが扱いやすいです。
アドストック／飽和（収穫逓減）のモデリング
MMMの核は2つの非線形変換です。どちらも「広告は素直に効かない」という現実を式に落としたものです。
アドストック（キャリーオーバー）。今週の出稿は今週だけでなく数週間効き続けます。幾何アドストックは過去の出稿を減衰させながら足し合わせます。
θ は減衰率で、値域は 0 < θ < 1。重みを総和1に正規化しておくと、定常的な出稿のもとでアドストック後の水準が出稿量と同じスケールに保たれ、後段の係数解釈が楽になります。L は打ち切り窓（週）で、8 週もあれば実務上は十分です。
飽和（収穫逓減）。出稿を倍にしても売上は倍にならない。最初の1万円と100万円目の1万円では効きが違う、という現象です。Hill 関数（ロジスティック飽和の一般化）が標準です。
K は半飽和点（最大効果の50%に達する出稿量）、s は曲線の立ち上がりの鋭さです。
順序は「アドストック → 飽和」。この2つの変換は交換しません（順序で結果が変わります）。出稿が各期で小さく累積で効くケースは先にアドストック、出稿が一部の期に集中する on/off 型は先に飽和、というのが Google の MMM 論文以来の経験則ですが、実務の大半は前者なので「アドストック先」を既定にします。
まず純粋な numpy で両変換を書いておきます。合成データの生成にも、後段の検算にも使います。
def adstock_np(x, theta, L=8):
w = theta ** np.arange(L)
w = w / w.sum()
out = np.zeros_like(x)
for i in range(len(x)):
lo = max(0, i - L + 1)
seg = x[lo:i + 1][::-1] # 直近を先頭に
out[i] = (seg * w[:len(seg)]).sum()
return out
def hill_np(x, K, s):
return x ** s / (x ** s + K ** s)真のパラメータを置いて売上を合成します。ベースライン（無施策でも立つ売上）を必ず分けて持つのがポイントです。
true_theta = np.array([0.30, 0.50, 0.15]) # アドストック減衰
true_K = np.array([120., 70., 45.]) # 半飽和点
true_slope = np.array([1.30, 1.00, 0.80]) # Hillの鋭さ
true_beta = np.array([18.0, 12.0, 7.0]) # チャネル係数
true_base = 60.0
L = 8
contrib_true = np.zeros((N, 3))
for j in range(3):
a = adstock_np(X_raw[:, j], true_theta[j], L)
contrib_true[:, j] = true_beta[j] * hill_np(a, true_K[j], true_slope[j])
mu = true_base + trend + seasonal + 9 * promo + contrib_true.sum(1)
sales = mu + rng.normal(0, 4, N) # 観測ノイズ実装（前処理 → 推定 → 寄与度分解）
ここからが本体です。PyMC で、上の2変換を pytensor 上に書き直し、ベイズ回帰として推定します。ベイズにする理由は「点推定ではなく区間が出る」こと。予算判断は「効いたか/効かなかったか」ではなく「どれくらいの確からしさで効いたか」で下すので、不確実性が出ないMMMは意思決定に使いにくいのです。
まずアドストックと Hill を pytensor で。numpy 版と同じ式を、自動微分が通る形で書き直します。
import pytensor.tensor as pt
import pymc as pm
import arviz as az
def geometric_adstock(x, theta, L):
w = pt.power(theta, pt.arange(L, dtype="float64"))
w = w / pt.sum(w)
padded = pt.concatenate([pt.zeros(L - 1), x])
# 各 t について直近 L 期を並べたラグ行列を作る
idx = pt.arange(x.shape[0])[:, None] + pt.arange(L)[None, ::-1]
lagged = padded[idx] # (N, L)
return pt.dot(lagged, w)
def hill(x, K, s):
xs = pt.power(x, s)
return xs / (xs + pt.power(K, s))モデル本体。事前分布に「広告効果は非負」「θ は0〜1」といった常識を埋め込みます。これがベイズの効きどころで、後述の通り K・s は弱識別になりやすいため、事前情報で締めておかないと区間が無意味に広がります。
coords = {"channel": channels, "obs": np.arange(N)}
with pm.Model(coords=coords) as model:
x_data = pm.Data("x_data", X_raw) # (N, 3)
trend_d = pm.Data("trend_d", t.astype("float64"))
season_d = pm.Data("season_d", seasonal)
promo_d = pm.Data("promo_d", promo)
base = pm.Normal("base", 60, 20)
b_trend = pm.Normal("b_trend", 0, 1)
b_seas = pm.Normal("b_season", 0, 1)
b_promo = pm.Normal("b_promo", 0, 10)
theta = pm.Beta("theta", 2, 2, dims="channel") # 0〜1
K = pm.HalfNormal("K", 100, dims="channel") # 非負
slope = pm.Gamma("slope", 3, 2, dims="channel") # 正
beta = pm.HalfNormal("beta", 15, dims="channel") # 非負（広告は売上を減らさない）
contribs = []
for j in range(len(channels)):
a = geometric_adstock(x_data[:, j], theta[j], L)
contribs.append(beta[j] * hill(a, K[j], slope[j]))
media = pt.stack(contribs, axis=1) # (N, 3)
pm.Deterministic("media_contrib", media, dims=("obs", "channel"))
mu = (base + b_trend * trend_d + b_seas * season_d
+ b_promo * promo_d + pt.sum(media, axis=1))
sigma = pm.HalfNormal("sigma", 10)
pm.Normal("y", mu=mu, sigma=sigma, observed=sales, dims="obs")
idata = pm.sample(1000, tune=1000, chains=4, target_accept=0.9,
random_seed=1)推定が終わったら、まず収束を確認します。r_hat が 1.01 以下、発散（divergences）がほぼゼロ、各パラメータの有効サンプル数が十分か。ここを飛ばして寄与度を読むのは禁物です。
print(az.summary(idata, var_names=["theta", "K", "slope", "beta", "base"]))寄与度分解は media_contrib の事後分布を平均して取り出します。「期間全体で各チャネルが売上を何ポイント押し上げたか」を、点ではなく区間つきで出せるのがベイズMMMの実利です。
# 事後平均の寄与（N×3）→ 期間合計のチャネル別寄与
contrib_post = idata.posterior["media_contrib"].mean(("chain", "draw")).values
contrib_sum = contrib_post.sum(axis=0) # チャネル別の累積寄与
base_sum = idata.posterior["base"].mean().item() * N
decomp = pd.DataFrame({
"channel": channels,
"contribution": contrib_sum,
"share_%": 100 * contrib_sum / sales.sum(),
})
print(decomp.round(2))
print(f"baseline share: {100 * base_sum / sales.sum():.1f}%")予算配分シミュレーション
寄与度がわかったら「次の予算をどう振るか」です。ここでは推定された応答曲線を使い、総予算固定のもとで売上期待を最大化する配分を、scipy の制約付き最適化（SLSQP）で求めます。
定常的な出稿のもとではアドストックの正規化重みの総和が1なので、一定出稿 s に対するアドストック後の水準は s に等しくなります。したがって週あたりの増分応答は、チャネルごとに beta * Hill(s) で近似できます。
from scipy.optimize import minimize
# 事後平均のパラメータ（idata から取り出した想定値）
theta_hat = np.array([0.22, 0.41, 0.24])
K_hat = np.array([105., 74., 97.])
slope_hat = np.array([1.52, 1.35, 1.69])
beta_hat = np.array([18.8, 12.8, 12.0])
def response(spend):
spend = np.asarray(spend, float)
xs = spend ** slope_hat
return beta_hat * (xs / (xs + K_hat ** slope_hat)) # チャネル別 週次増分
current = np.array([100., 60., 40.]) # 現状の週次配分
total = current.sum()
# 推定が確かなレンジ（±50%）に探索を制限し、外挿を避ける
bounds = list(zip(current * 0.5, current * 1.5))
cons = ({"type": "eq", "fun": lambda s: s.sum() - total},)
res = minimize(lambda s: -response(s).sum(), current,
method="SLSQP", bounds=bounds, constraints=cons)
opt = res.x
lift = response(opt).sum() / response(current).sum() - 1
print("最適配分:", opt.round(1), " 合計:", opt.sum().round(1))
print(f"現状比リフト: {lift * 100:.1f}%")ここで重要なのは探索範囲を現状の±50%程度に縛ることです。MMMは観測した出稿レンジの外では応答曲線を保証しません。「使ったことのない予算水準」へ最適化が飛ぶと、モデルが知らない領域へ外挿してしまい、提案が空想になります。
もう一つ実務で効くのが**限界効果（marginal ROAS）**の確認です。現在の出稿点で各チャネルの応答曲線の傾きを見れば、「次の1単位をどこに足すと一番効くか」が一目でわかります。
eps = 1e-3
mroas = (response(current + eps) - response(current)) / eps
for c, m in zip(channels, mroas):
print(f"{c}: 限界効果 {m:.4f}（出稿1単位あたり増分売上）")なお、この記事の最適化は「固定の応答曲線に対する一回の解」です。事後分布の不確実性まで取り込んで「探索と活用」を回す配分最適化（ベイズ最適化）は範囲が別なので、ベイズ最適化で広告配分を最適化する実装 に分けています。MMMの推定結果を「次の一手」に変える運用は、そちらが本題です。
限界効果（曲線の傾き）が大きいチャネルへ寄せると合計応答が伸びる。やがて 収穫逓減で傾きが揃い、再配分の旨味は消える。MMMが見せるのはこの「均衡点」で、 劇的なリフトではなく勘で振っていた配分への根拠づけが本体。
まとめ
MMMの自作に必要なのは、アドストックと飽和という2つの変換、それをベイズ回帰に載せる骨格、収束と弱識別への目配り、そして外挿を避けた配分最適化——この4点でした。専用ライブラリは生産性を上げますが、一度自分で組んでおくと、出てきた寄与度のどこを信じ、どこを疑うべきかが体でわかります。
「そもそもMMMで何が決められて何が決められないのか」という上流の判断は なぜMMMが経営に効くのか に、既製OSSとの使い分けは Robyn vs LightweightMMM vs 自作 にあります。