import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
from scipy.special import gammaln, betaln
import matplotlib.pyplot as plt
from IPython.display import display_markdown
%config InlineBackend.figure_formats = ['svg']
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
pd.set_option("display.precision", 2)
PANEL_SIZE = 5021 # static panel
YEAR = 1 # weeks 1-52
NSTAR = 170 # truncation for the sums over n (justified below)NBD-Dirichlet Applied to Consumer Panel Data
Applies the Dirichlet model to the edible grocery panel used in Analysing Buyer Behaviour Using Consumer Panel Data — 5,021 static panellists, two years, five brands (Alpha, Bravo, Charlie, Delta, Other).
The derivations, the literature and the marketing argument are in NBD-Dirichlet Model of Consumer Buying Behavior, which reproduces the 1984 toothpaste example.
1 Formula reference card
Four assumption levels: category purchases Poisson(\(\mu\)) × gamma across consumers → NBD; brand choice multinomial(\(\mathbf p\)) × Dirichlet(\(\boldsymbol\alpha\)) across consumers → Dirichlet-multinomial; the two independent.
Three parameters: \(M\) (per-capita category rate), \(K\) (diversity in how often), \(S = \sum_j\alpha_j\) (diversity in which brand). Brand \(j\) has \(\alpha_j = S\mu_j\) where \(\mu_j\) is market share.
| Quantity | Formula |
|---|---|
| Category incidence | \(P_n = \left(1+\tfrac{MT}{K}\right)^{-K}\tfrac{\Gamma(K+n)}{n!\,\Gamma(K)}\left(\tfrac{MT}{MT+K}\right)^{n}\) |
| Fit \(M\) | \(M = B \times W\) (category penetration × category buy rate) |
| Fit \(K\) (“mean and zeros”) | solve \(1-B = (1+M/K)^{-K}\) |
| Brand choice given \(n\) | \(p(r\mid n)=\binom{n}{r}\dfrac{B(\alpha_j+r,\;S-\alpha_j+n-r)}{B(\alpha_j,\;S-\alpha_j)}\) |
| — the \(r=0\) case | \(p(0\mid n)=\prod_{a=0}^{n-1}\dfrac{S-\alpha_j+a}{S+a}\), \(\;p(0\mid 0)=1\) |
| Fit \(S\) | per brand, solve \(1-\sum_n P_n\,p(0\mid n) = b_j^{\text{obs}}\); pool by share-weighted mean, dropping outliers |
| Penetration | \(b_j = 1-\sum_{n\ge0}P_n\,p(0\mid n)\) |
| Purchases per brand buyer | \(w_j = \tfrac{1}{b_j}\sum_{n\ge1}P_n\sum_{r=1}^{n}r\,p(r\mid n)\) |
| Category purchases per brand buyer | \(w_{P,j} = \tfrac{1}{b_j}\sum_{n\ge1}n\,P_n\,[1-p(0\mid n)]\) |
| Brand purchase frequency dist. | \(f_j(r)=\sum_{n\ge r}P_n\,p(r\mid n)\) |
| Restricted to category buyers in \(R\) | \(b_{j\mid R}=1-\tfrac{\sum_{n\in R}P_n p(0\mid n)}{\sum_{n\in R}P_n}\) |
| Duplication | \(b_{(j+k)}=1-\sum_n P_n\,p_{j+k}(0\mid n)\); \(\;b_{jk}=b_j+b_k-b_{(j+k)}\); \(\;b_{j\mid k}=b_{jk}/b_k\) |
| Time extrapolation | replace \(M \to MT\). \(K\) and \(S\) are time-invariant |
Additivity is what makes all of this computable: merging brands adds their \(\alpha\), so “brand \(j\) vs everything else” is a two-category Dirichlet-multinomial, i.e. a beta-binomial.
2 Imports
3 Deriving the model inputs from the panel
The Dirichlet needs four things. The panel gives transactions, so we aggregate.
A purchase occasion is a unique (panel_id, trans_id). A trip on which two brands were bought counts as one category occasion but as one occasion for each brand — the convention used in the source notebook.
df = pd.read_csv("data/panel-datasets/edible_grocery.csv", encoding="utf-8-sig")
yr = df[(df.week > (YEAR - 1) * 52) & (df.week <= YEAR * 52)]
brand_occ = yr[["panel_id", "trans_id", "brand"]].drop_duplicates()
cat_occ = yr[["panel_id", "trans_id"]].drop_duplicates()
brand_name = np.array(sorted(brand_occ.brand.unique()))
nbrand = len(brand_name)
grp = brand_occ.groupby("brand")
occasions = grp.size().loc[brand_name]
buyers = grp.panel_id.nunique().loc[brand_name]
brand_share = (occasions / occasions.sum()).to_numpy()
brand_pen_obs = (buyers / PANEL_SIZE).to_numpy()
cat_pen = cat_occ.panel_id.nunique() / PANEL_SIZE
cat_buyrate = len(cat_occ) / cat_occ.panel_id.nunique()
pd.DataFrame({"occasions": occasions, "buyers": buyers,
"share": brand_share, "pen (obs)": brand_pen_obs},
index=brand_name)| occasions | buyers | share | pen (obs) | |
|---|---|---|---|---|
| Alpha | 9060 | 2624 | 0.44 | 0.52 |
| Bravo | 8255 | 2562 | 0.40 | 0.51 |
| Charlie | 1882 | 813 | 0.09 | 0.16 |
| Delta | 859 | 380 | 0.04 | 0.08 |
| Other | 422 | 176 | 0.02 | 0.04 |
Code
display_markdown(
f"Category penetration **{cat_pen:.4f}**, buy rate **{cat_buyrate:.4f}** "
f"→ $M$ = **{cat_pen * cat_buyrate:.4f}**\n\n"
f"Brand occasions {occasions.sum():,} vs category occasions {len(cat_occ):,} "
f"— **{100 * (occasions.sum() - len(cat_occ)) / len(cat_occ):.2f}%** of trips "
f"involved more than one brand.",
raw=True,
)Category penetration 0.9110, buy rate 4.3791 → \(M\) = 3.9892
Brand occasions 20,478 vs category occasions 20,030 — 2.24% of trips involved more than one brand.
A modelling wrinkle to note up front. The Dirichlet assumes every category purchase goes to exactly one brand, so \(\sum_j r_j = n\). Here 2.2% of trips break that — a household bought Alpha and Delta on one trip. \(M\) is therefore built from category occasions while shares are built from brand occasions, and the two bases differ by 2.2%. Small enough to proceed, but it is a real (if minor) violation of the model’s own arithmetic, not just a rounding matter.
Compare with the toothpaste example: category penetration is 0.91 here against 0.56, and two brands hold 85% of the category. This is a much more concentrated, much more widely bought category — a genuinely different test.
4 Fitting
4.1 \(M\) and \(K\)
def estimate_K(M, cat_pen, max_K=30.0):
"""Mean-and-zeros: solve (1 + M/K)^-K = 1 - B for K."""
cp = np.log(1.0 - cat_pen)
return minimize_scalar(lambda K: (K * np.log(1.0 + M / K) + cp) ** 2,
bounds=(1e-4, max_K), method="bounded").x
M0 = cat_pen * cat_buyrate
K = estimate_K(M0, cat_pen)
print(f"M = {M0:.4f} K = {K:.4f}")
print(f"existence condition M > -log(P0): {M0:.3f} > {-np.log(1 - cat_pen):.3f} -> "
f"{M0 > -np.log(1 - cat_pen)}")
print(f"implied P(0) = {(1 + M0 / K) ** -K:.4f} observed = {1 - cat_pen:.4f}")M = 3.9892 K = 2.6031
existence condition M > -log(P0): 3.989 > 2.419 -> True
implied P(0) = 0.0890 observed = 0.0890\(K = 2.6\) against 0.78 for toothpaste. Higher \(K\) means less dispersion in category buying rates — consistent with a category that 91% of households buy.
4.2 \(P_n\) and the truncation check
def Pn(M, K, n):
"""NBD probability of n category purchases, computed in logs."""
n = np.atleast_1d(np.asarray(n, dtype=int))
a = np.arange(int(n.max()))
logratio = np.concatenate(([0.0], np.cumsum(np.log(K + a) - np.log(1.0 + a))))
return np.exp(-K * np.log(1.0 + M / K)
+ logratio[n]
+ np.where(n > 0, n * np.log(M / (M + K)), 0.0))
def check_truncation(M, K, nstar=NSTAR):
ns = np.arange(nstar + 1)
P = Pn(M, K, ns)
return {"sum P(n)": P.sum(), "implied mean": (ns * P).sum(), "target M": M}
pd.DataFrame([check_truncation(M0 * t, K) | {"t": t} for t in (1, 2, 4)]).set_index("t")| sum P(n) | implied mean | target M | |
|---|---|---|---|
| t | |||
| 1 | 1.0 | 3.99 | 3.99 |
| 2 | 1.0 | 7.98 | 7.98 |
| 4 | 1.0 | 15.96 | 15.96 |
The default nstar = 50 in the R package is not safe here. With \(M \approx 4\) per year, the two-year and four-year sums need far more mass; NSTAR = 170 keeps \(\sum P_n\) at 1.000 and recovers the mean at every horizon used below. This diagnostic is the way to choose it — a too-small nstar silently biases every measure downward.
4.3 The brand-choice kernel
The whole beta-binomial is precomputed once per brand as a matrix \(Q[n, r] = p(r \mid n)\), upper-triangular since \(r \le n\).
def pr_matrix(alpha, S, nstar=NSTAR):
"""(nstar+1, nstar+1) matrix of p(r | n); zero where r > n.
log p(r|n) = log C(n,r) + log B(a+r, S-a+n-r) - log B(a, S-a)
"""
n = np.arange(nstar + 1)[:, None]
r = np.arange(nstar + 1)[None, :]
ok = r <= n
nr = np.where(ok, n - r, 0) # mask before gammaln to avoid poles
log_q = (gammaln(n + 1) - gammaln(r + 1) - gammaln(nr + 1)
+ betaln(alpha + r, S - alpha + nr) - betaln(alpha, S - alpha))
return np.where(ok, np.exp(log_q), 0.0)
def p_zero_vec(alpha, S, nstar=NSTAR):
"""p(0 | n) for n = 0..nstar via the telescoped product, as a cumulative sum."""
a = np.arange(nstar)
return np.concatenate(([1.0], np.exp(np.cumsum(np.log(S - alpha + a)
- np.log(S + a)))))4.4 \(S\)
def S_for_brand(j, M, K, share, pen_obs, nstar=NSTAR, max_S=30.0):
P = Pn(M, K, np.arange(nstar + 1))
def obj(S):
return (1.0 - np.sum(P * p_zero_vec(S * share[j], S, nstar)) - pen_obs[j]) ** 2
return minimize_scalar(obj, bounds=(1e-9, max_S), method="bounded").x
def fivenum(x):
"""R's stats::fivenum — Tukey hinges, not interpolated quartiles."""
x = np.sort(np.asarray(x, dtype=float))
n = len(x)
n4 = np.floor((n + 3) / 2) / 2
d = np.array([1.0, n4, (n + 1) / 2, n + 1 - n4, float(n)])
return 0.5 * (x[np.floor(d).astype(int) - 1] + x[np.ceil(d).astype(int) - 1])
def pool_S(S_all, share):
"""Share-weighted mean, dropping (a) 1.5*IQR outliers and (b) anything above
the upper boxplot notch. Both tests come from R's boxplot.stats(), which the
NBDdirichlet package calls but does not implement itself."""
_, hinge_lo, median, hinge_hi, _ = fivenum(S_all)
iqr = hinge_hi - hinge_lo
fence = (S_all < hinge_lo - 1.5 * iqr) | (S_all > hinge_hi + 1.5 * iqr)
notch_hi = median + 1.58 * iqr / np.sqrt(len(S_all))
drop = fence | (S_all > notch_hi)
return np.average(S_all[~drop], weights=share[~drop]), drop, notch_hi
S_all = np.array([S_for_brand(j, M0, K, brand_share, brand_pen_obs)
for j in range(nbrand)])
S, dropped, notch_hi = pool_S(S_all, brand_share)
display_markdown(
f"Upper notch **{notch_hi:.4f}** → dropped **"
+ (", ".join(brand_name[dropped]) if dropped.any() else "none")
+ f"** → pooled **S = {S:.4f}**",
raw=True,
)
pd.DataFrame({"share": brand_share, "pen (obs)": brand_pen_obs,
"S_j": S_all, "dropped": dropped}, index=brand_name)Upper notch 1.1689 → dropped Charlie → pooled S = 0.6171
| share | pen (obs) | S_j | dropped | |
|---|---|---|---|---|
| Alpha | 0.44 | 0.52 | 0.50 | False |
| Bravo | 0.40 | 0.51 | 0.68 | False |
| Charlie | 0.09 | 0.16 | 1.17 | True |
| Delta | 0.04 | 0.08 | 1.11 | False |
| Other | 0.02 | 0.04 | 0.87 | False |
\(S = 0.62\), roughly half the toothpaste value of 1.30. Low \(S\) means a polarised category: consumers have strongly differentiated brand propensities rather than picking near-randomly in proportion to share. That matches the source notebook’s finding that over two-thirds of category buyers bought only one brand all year.
Code
display_markdown(
f"| | $M$ | $K$ | $S$ |\n|---|---|---|---|\n"
f"| edible grocery (year, 5 brands) | {M0:.3f} | {K:.3f} | {S:.3f} |\n"
f"| toothpaste (quarter, 8 brands) | 1.456 | 0.778 | 1.295 |",
raw=True,
)| \(M\) | \(K\) | \(S\) | |
|---|---|---|---|
| edible grocery (year, 5 brands) | 3.989 | 2.603 | 0.617 |
| toothpaste (quarter, 8 brands) | 1.456 | 0.778 | 1.295 |
5 Brand performance measures
def brand_measures(M, K, S, alpha, nstar=NSTAR, limit=None):
"""(penetration, unnormalised E[brand purchases], unnormalised E[cat purchases])."""
ns = np.arange(nstar + 1) if limit is None else np.asarray(limit, dtype=int)
P = Pn(M, K, ns)
Q = pr_matrix(alpha, S, nstar)[ns]
r = np.arange(nstar + 1)
pen = 1.0 - np.sum(P * Q[:, 0])
return pen, np.sum(P * (Q * r).sum(axis=1)), np.sum(ns * P * (1.0 - Q[:, 0]))
def buy_table(M, K, S, share, names, nstar=NSTAR):
rows = []
for j in range(len(share)):
b, e_brand, e_cat = brand_measures(M, K, S, S * share[j], nstar)
rows.append([b, e_brand / b, e_cat / b])
return pd.DataFrame(rows, index=names,
columns=["pen.brand", "pur.brand", "pur.cat"])
buy = buy_table(M0, K, S, brand_share, brand_name)
buy.round(2)| pen.brand | pur.brand | pur.cat | |
|---|---|---|---|
| Alpha | 0.54 | 3.27 | 4.75 |
| Bravo | 0.50 | 3.20 | 4.78 |
| Charlie | 0.14 | 2.70 | 5.01 |
| Delta | 0.06 | 2.63 | 5.05 |
| Other | 0.03 | 2.60 | 5.07 |
5.1 Observed versus theoretical
The point of using real panel data rather than the paper’s published summary is that we can compute the observeds ourselves.
# panellist x brand matrix of purchase occasions
X = (brand_occ.groupby(["panel_id", "brand"]).size()
.unstack(fill_value=0).reindex(columns=brand_name, fill_value=0))
n_cat = cat_occ.groupby("panel_id").size()
X = X.reindex(n_cat.index, fill_value=0)
ever = (X > 0).astype(int)
obs = pd.DataFrame({
"b (O)": ever.sum() / PANEL_SIZE,
"w (O)": X.sum() / ever.sum(),
"wp (O)": [n_cat[ever[b] == 1].mean() for b in brand_name],
}, index=brand_name)
ot = pd.DataFrame({
"share": brand_share,
"b (O)": obs["b (O)"], "b (T)": buy["pen.brand"],
"w (O)": obs["w (O)"], "w (T)": buy["pur.brand"],
"wp (O)": obs["wp (O)"], "wp (T)": buy["pur.cat"],
})
ot.loc["Average"] = ot.mean()
ot.round(2)| share | b (O) | b (T) | w (O) | w (T) | wp (O) | wp (T) | |
|---|---|---|---|---|---|---|---|
| Alpha | 0.44 | 0.52 | 0.54 | 3.45 | 3.27 | 4.88 | 4.75 |
| Bravo | 0.40 | 0.51 | 0.50 | 3.22 | 3.20 | 4.59 | 4.78 |
| Charlie | 0.09 | 0.16 | 0.14 | 2.31 | 2.70 | 5.24 | 5.01 |
| Delta | 0.04 | 0.08 | 0.06 | 2.26 | 2.63 | 5.58 | 5.05 |
| Other | 0.02 | 0.04 | 0.03 | 2.40 | 2.60 | 5.96 | 5.07 |
| Average | 0.20 | 0.26 | 0.25 | 2.73 | 2.88 | 5.25 | 4.93 |
Read it the way Ehrenberg would — down the columns, at two significant figures.
- Penetration is close for the two big brands (0.52/0.54 and 0.51/0.50) but the model under-predicts the small brands: Charlie 0.16 observed against 0.14 theoretical, Delta 0.08 against 0.06.
- \(w\) goes the other way: the model over-predicts loyalty for small brands (Charlie 2.3 observed vs 2.7 theoretical).
- \(w_P\) shows the natural-monopoly trend rising as share falls in both columns, but the observed trend is steeper (4.9 → 6.0) than the theoretical (4.8 → 5.1).
So the fit is materially worse than the toothpaste example, where penetration came within a point across all eight brands. That is worth taking seriously rather than reporting an average and moving on.
5.2 Double Jeopardy
Code
dj = pd.DataFrame({
"share": brand_share,
"b (O)": obs["b (O)"], "w (O)": obs["w (O)"],
"w(1-b) (O)": obs["w (O)"] * (1 - obs["b (O)"]),
"w(1-b) (T)": buy["pur.brand"] * (1 - buy["pen.brand"]),
})
dj.round(2)| share | b (O) | w (O) | w(1-b) (O) | w(1-b) (T) | |
|---|---|---|---|---|---|
| Alpha | 0.44 | 0.52 | 3.45 | 1.65 | 1.51 |
| Bravo | 0.40 | 0.51 | 3.22 | 1.58 | 1.60 |
| Charlie | 0.09 | 0.16 | 2.31 | 1.94 | 2.34 |
| Delta | 0.04 | 0.08 | 2.26 | 2.09 | 2.46 |
| Other | 0.02 | 0.04 | 2.40 | 2.31 | 2.52 |
Both the observed and theoretical w(1-b) columns are far flatter than either \(b\) or \(w\) alone — the Double Jeopardy relation holds in the data. The observed column runs 1.6–2.2 while share varies twenty-fold; the theoretical runs 1.5–2.5. Neither is constant, but the residual variation is a fraction of the variation the relation absorbs.
6 The summary tables
6.1 Purchase frequency distribution
def freq_table(M, K, S, share, names, cutoff=5, nstar=NSTAR):
P = Pn(M, K, np.arange(nstar + 1))
rows = []
for j in range(len(share)):
f = (P[:, None] * pr_matrix(S * share[j], S, nstar)).sum(axis=0)
rows.append(list(f[:cutoff + 1]) + [f[cutoff + 1:].sum()])
cols = [str(i) for i in range(cutoff + 1)] + [f"{cutoff + 1}+"]
return pd.DataFrame(rows, index=names, columns=cols)
theo_freq = freq_table(M0, K, S, brand_share, brand_name)
theo_freq.round(2)| 0 | 1 | 2 | 3 | 4 | 5 | 6+ | |
|---|---|---|---|---|---|---|---|
| Alpha | 0.46 | 0.16 | 0.11 | 0.08 | 0.06 | 0.04 | 0.09 |
| Bravo | 0.50 | 0.16 | 0.10 | 0.07 | 0.05 | 0.04 | 0.08 |
| Charlie | 0.86 | 0.05 | 0.03 | 0.02 | 0.01 | 0.01 | 0.01 |
| Delta | 0.94 | 0.03 | 0.01 | 0.01 | 0.01 | 0.00 | 0.01 |
| Other | 0.97 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 |
Code
obs_freq = pd.DataFrame(
{b: X[b].clip(upper=6).value_counts().reindex(range(7), fill_value=0)
for b in brand_name}).T
obs_freq[0] += PANEL_SIZE - len(X) # households that bought no category at all
obs_freq = obs_freq / PANEL_SIZE
obs_freq.columns = [str(i) for i in range(6)] + ["6+"]
pd.concat({"Observed": obs_freq.round(2), "Theoretical": theo_freq.round(2)},
axis=1).swaplevel(axis=1).sort_index(axis=1)| 0 | 1 | 2 | 3 | 4 | 5 | 6+ | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Observed | Theoretical | Observed | Theoretical | Observed | Theoretical | Observed | Theoretical | Observed | Theoretical | Observed | Theoretical | Observed | Theoretical | |
| Alpha | 0.48 | 0.46 | 0.15 | 0.16 | 0.10 | 0.11 | 0.08 | 0.08 | 0.06 | 0.06 | 0.05 | 0.04 | 0.09 | 0.09 |
| Bravo | 0.49 | 0.50 | 0.14 | 0.16 | 0.10 | 0.10 | 0.08 | 0.07 | 0.06 | 0.05 | 0.04 | 0.04 | 0.08 | 0.08 |
| Charlie | 0.84 | 0.86 | 0.08 | 0.05 | 0.03 | 0.03 | 0.02 | 0.02 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 |
| Delta | 0.92 | 0.94 | 0.04 | 0.03 | 0.01 | 0.01 | 0.01 | 0.01 | 0.00 | 0.01 | 0.00 | 0.00 | 0.01 | 0.01 |
| Other | 0.96 | 0.97 | 0.02 | 0.01 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
The reverse-J shape is reproduced for every brand. The visible miss is at \(r = 1\) for the small brands — Charlie has 8% of households buying once against a predicted 5%. Too many light buyers relative to the model, the same signature Ehrenberg discusses as an “excess of occasional buyers”.
6.2 Heavy versus light category buyers
def heavy_table(M, K, S, share, names, limit, nstar=NSTAR):
limit = np.asarray(limit, dtype=int)
P_sum = Pn(M, K, limit).sum()
rows = []
for j in range(len(share)):
b_lim, e_brand, _ = brand_measures(M, K, S, S * share[j], nstar, limit=limit)
not_buying = 1.0 - b_lim
rows.append([1.0 - not_buying / P_sum, e_brand / (P_sum - not_buying)])
return pd.DataFrame(rows, index=names,
columns=["Penetration", "Avg Purchase Freq"])
pd.concat({"light (1-6)": heavy_table(M0, K, S, brand_share, brand_name, range(1, 7)),
"heavy (7+)": heavy_table(M0, K, S, brand_share, brand_name,
range(7, NSTAR + 1))}, axis=1).round(2)| light (1-6) | heavy (7+) | |||
|---|---|---|---|---|
| Penetration | Avg Purchase Freq | Penetration | Avg Purchase Freq | |
| Alpha | 0.57 | 2.45 | 0.69 | 5.90 |
| Bravo | 0.53 | 2.41 | 0.65 | 5.71 |
| Charlie | 0.14 | 2.10 | 0.19 | 4.39 |
| Delta | 0.06 | 2.05 | 0.09 | 4.21 |
| Other | 0.03 | 2.03 | 0.05 | 4.13 |
Heavy category buyers are more likely to buy every brand, and to buy each more often — a selection effect, not a segment.
6.3 Duplication of purchase — where the model breaks
def dup_table(M, K, S, share, names, focal, nstar=NSTAR):
P = Pn(M, K, np.arange(nstar + 1))
b_k = brand_measures(M, K, S, S * share[focal], nstar)[0]
out = np.zeros(len(share))
out[focal] = 1.0
for j in range(len(share)):
if j == focal:
continue
# additivity: the composite brand (j+k) has alpha = alpha_j + alpha_k
b_comp = 1.0 - np.sum(P * pr_matrix(S * (share[focal] + share[j]), S, nstar)[:, 0])
b_j = brand_measures(M, K, S, S * share[j], nstar)[0]
out[j] = (b_j + b_k - b_comp) / b_k
return pd.Series(out, index=names, name=names[focal])
theo_dup = pd.concat([dup_table(M0, K, S, brand_share, brand_name, f)
for f in range(nbrand)], axis=1).T
theo_dup.round(2)| Alpha | Bravo | Charlie | Delta | Other | |
|---|---|---|---|---|---|
| Alpha | 1.00 | 0.39 | 0.1 | 0.05 | 0.02 |
| Bravo | 0.41 | 1.00 | 0.1 | 0.05 | 0.02 |
| Charlie | 0.40 | 0.37 | 1.0 | 0.05 | 0.02 |
| Delta | 0.40 | 0.37 | 0.1 | 1.00 | 0.02 |
| Other | 0.40 | 0.37 | 0.1 | 0.05 | 1.00 |
Code
counts = ever.T.values @ ever.values # duplication count table
obs_dup = pd.DataFrame(counts / np.diag(counts)[:, None],
index=brand_name, columns=brand_name)
obs_dup.round(2)| Alpha | Bravo | Charlie | Delta | Other | |
|---|---|---|---|---|---|
| Alpha | 1.00 | 0.34 | 0.15 | 0.09 | 0.03 |
| Bravo | 0.35 | 1.00 | 0.15 | 0.05 | 0.04 |
| Charlie | 0.50 | 0.47 | 1.00 | 0.14 | 0.03 |
| Delta | 0.63 | 0.37 | 0.31 | 1.00 | 0.03 |
| Other | 0.39 | 0.62 | 0.15 | 0.06 | 1.00 |
Code
resid = (obs_dup - theo_dup).where(~np.eye(nbrand, dtype=bool))
resid.round(2)| Alpha | Bravo | Charlie | Delta | Other | |
|---|---|---|---|---|---|
| Alpha | NaN | -0.04 | 0.05 | 0.04 | 0.00 |
| Bravo | -0.06 | NaN | 0.05 | 0.01 | 0.02 |
| Charlie | 0.09 | 0.10 | NaN | 0.10 | 0.01 |
| Delta | 0.23 | -0.00 | 0.21 | NaN | 0.00 |
| Other | -0.01 | 0.25 | 0.05 | 0.01 | NaN |
The residuals are structured rather than noisy — but this table overstates part of the story, and it is worth being careful about why.
Two things are mixed together here. The Dirichlet’s own \(b_j\) is off for the small brands (§ Observed versus theoretical), and any duplication computed from a wrong \(b_j\) inherits that error. To separate a genuine duplication anomaly from penetration miss propagating, compare against the Duplication of Purchase Law fitted to the observed penetrations, with no Dirichlet involved:
\[b_{jk} = D\,b_j\,b_k\]
Ehrenberg (Repeat-Buying §10.5, footnote) recommends exactly this symmetric form over the conditional \(b_{j\mid k}\), because it is symmetric and homoscedastic — the residuals do not grow with the size of the brands being compared.
off = ~np.eye(nbrand, dtype=bool)
b_obs = obs["b (O)"].values
b_jk = (ever.T.values @ ever.values) / PANEL_SIZE # joint penetration
D = b_jk[off].sum() / np.outer(b_obs, b_obs)[off].sum()
sym_resid = pd.DataFrame(100 * (b_jk - D * np.outer(b_obs, b_obs)),
index=brand_name, columns=brand_name).where(off)
display_markdown(f"Single duplication coefficient **D = {D:.3f}** — "
"residuals below in **percentage points of the panel**", raw=True)
sym_resid.round(2)Single duplication coefficient D = 0.837 — residuals below in percentage points of the panel
| Alpha | Bravo | Charlie | Delta | Other | |
|---|---|---|---|---|---|
| Alpha | NaN | -4.48 | 0.94 | 1.45 | -0.16 |
| Bravo | -4.48 | NaN | 0.69 | -0.44 | 0.69 |
| Charlie | 0.94 | 0.69 | NaN | 1.28 | 0.04 |
| Delta | 1.45 | -0.44 | 1.28 | NaN | -0.02 |
| Other | -0.16 | 0.69 | 0.04 | -0.02 | NaN |
On this scale one residual dominates everything else:
- Alpha ↔︎ Bravo: −4.5 points. The two leaders share far fewer buyers than their sizes imply.
- Every other pair falls between −0.4 and +1.5, averaging about +0.5.
So the honest reading is narrower than “leaders under-duplicate, followers over-duplicate”. The follower-side excess is mostly the Dirichlet’s penetration miss showing through; against the observed penetrations it largely disappears. The one robust anomaly is that Alpha and Bravo behave as near-mutually-exclusive alternatives.
That corroborates something found independently in the source notebook: the weekly volume shares of Alpha and Bravo are strongly negatively correlated. Two brands trading share back and forth, while sharing few buyers, is one phenomenon seen two ways.
6.4 Sole-brand buyers — the sharpest miss
The proportion of the population buying only brand \(j\) all year is \(\sum_{n\ge1}P_n\,p(n\mid n)\), and since \(p(n\mid n)=B(\alpha_j+n,\,S-\alpha_j)/B(\alpha_j,\,S-\alpha_j)\) it is cheap to compute. The sum must start at \(n=1\) — \(p(0\mid 0)=1\), so including \(n=0\) would count every category non-buyer as a “sole buyer”.
from scipy.special import betaln as _betaln
pos = np.arange(1, NSTAR + 1)
P_pos = Pn(M0, K, np.arange(NSTAR + 1))[1:]
def sole_theoretical(alpha, S):
return np.sum(P_pos * np.exp(_betaln(alpha + pos, S - alpha)
- _betaln(alpha, S - alpha)))
n_brands_bought = ever.sum(axis=1)
sole_obs = np.array([((ever[b] == 1) & (n_brands_bought == 1)).sum() / PANEL_SIZE
for b in brand_name])
sole_theo = np.array([sole_theoretical(S * s, S) for s in brand_share])
pd.DataFrame({
"share": brand_share,
"sole (O)": sole_obs, "sole (T)": sole_theo,
"sole/pen (O)": sole_obs / obs["b (O)"].values,
"sole/pen (T)": sole_theo / buy["pen.brand"].values,
}, index=brand_name).map("{:.3f}".format)| share | sole (O) | sole (T) | sole/pen (O) | sole/pen (T) | |
|---|---|---|---|---|---|
| Alpha | 0.442 | 0.270 | 0.272 | 0.517 | 0.504 |
| Bravo | 0.403 | 0.270 | 0.241 | 0.529 | 0.479 |
| Charlie | 0.092 | 0.031 | 0.044 | 0.194 | 0.324 |
| Delta | 0.042 | 0.014 | 0.019 | 0.187 | 0.304 |
| Other | 0.021 | 0.005 | 0.009 | 0.131 | 0.296 |
Look at the last two columns — the share of a brand’s own buyers who bought nothing else all year:
- Leaders: essentially exact. Alpha 0.52 observed against 0.50 theoretical; Bravo 0.53 against 0.48.
- Followers: over-predicted by about 60%. Charlie 0.19 observed against 0.32 theoretical; Delta 0.19 against 0.30; Other 0.13 against 0.30.
The model thinks small brands should have a proportionally loyal core. They do not.
7 Diagnosis: which assumptions actually fail here
Four independent measurements now point the same way for the small brands. Each is a different quantity, and all four are consistent with one story:
| Measure | Observed | Theoretical | Direction |
|---|---|---|---|
| penetration \(b\) (Charlie) | 0.16 | 0.14 | more buyers than predicted |
| purchase rate \(w\) (Charlie) | 2.3 | 2.7 | who buy it less often |
| category rate \(w_P\) (Charlie) | 5.2 | 5.0 | and buy the category more |
| sole loyalty (Charlie) | 0.19 | 0.32 | and are far less exclusive |
That is a coherent description: the small brands are a repertoire add-on for heavy category buyers, not a proportionally scaled-down version of a leader. Meanwhile the two leaders behave like near-monogamous alternatives — about half of each one’s buyers bought nothing else all year, and they share 4.5 points fewer buyers than their sizes imply.
7.1 Ruling things out first
Before blaming a modelling assumption it is worth checking the cheap explanations.
def year_summary(y):
d = df[(df.week > (y - 1) * 52) & (df.week <= y * 52)]
o = d[["panel_id", "trans_id", "brand"]].drop_duplicates()
c = d[["panel_id", "trans_id"]].drop_duplicates()
g = o.groupby("brand")
return pd.DataFrame({"share": (g.size() / g.size().sum()).loc[brand_name],
"pen": (g.panel_id.nunique() / PANEL_SIZE).loc[brand_name]}), \
c.panel_id.nunique() / PANEL_SIZE
y1, cp1 = year_summary(1)
y2, cp2 = year_summary(2)
display_markdown(f"Category penetration: year 1 **{cp1:.3f}**, year 2 **{cp2:.3f}**. "
f"Largest share change **{(y2.share - y1.share).abs().max():.3f}**.",
raw=True)
pd.concat({"year 1": y1, "year 2": y2}, axis=1).round(3)Category penetration: year 1 0.911, year 2 0.912. Largest share change 0.018.
| year 1 | year 2 | |||
|---|---|---|---|---|
| share | pen | share | pen | |
| brand | ||||
| Alpha | 0.44 | 0.52 | 0.46 | 0.55 |
| Bravo | 0.40 | 0.51 | 0.41 | 0.50 |
| Charlie | 0.09 | 0.16 | 0.08 | 0.14 |
| Delta | 0.04 | 0.08 | 0.04 | 0.07 |
| Other | 0.02 | 0.04 | 0.02 | 0.04 |
Non-stationarity is ruled out. Category penetration is identical across the two years (0.911 → 0.912) and no brand’s share moves more than 1.8 points. Whatever is wrong, it is not that the market is in motion. (Charlie drifts down 9.2% → 7.7%, which is the largest movement and worth noting, but far too small to explain a 60% miss on sole loyalty.)
The multi-brand-trip wrinkle is too small. 2.2% of occasions allocate to more than one brand, violating \(\sum_j r_j = n\). Real, but an order of magnitude smaller than the effects above.
7.2 The assumption that fails: A2
Mapping onto the five assumptions:
| Assumption | Verdict here | |
|---|---|---|
| B1 | Poisson category incidence | Plausible. Frequently-bought grocery, 4 purchases/year, no seasonality or stockpiling visible in the annual aggregates. |
| B2 | Gamma-distributed rates | Plausible. \(K=2.6\) reproduces observed non-buyers exactly by construction; the category frequency distribution is well-behaved. |
| A1 | Multinomial brand choice | Mostly holds, apart from the 2.2% of multi-brand trips. |
| A2 | Dirichlet-distributed choice vectors | Fails. This is the one. |
| C | Rate ⟂ choice | Partly implicated — small brands are disproportionately bought by heavy category buyers, which is a mild dependence between \(\mu\) and \(\mathbf p\). |
The Dirichlet is the unique distribution expressing “independence except for the constraint \(\sum_j p_j = 1\)”. Fitting one says the category is unsegmented: your propensity for Alpha carries no information about how you split the remainder.
Here it plainly does. The population looks closer to a two-camp mixture — a mass of households with \(\mathbf p\) concentrated on Alpha, another concentrated on Bravo, with the small brands picked up occasionally by heavy buyers from either camp. A single Dirichlet cannot represent that: to reproduce the leaders’ high sole loyalty it must drive \(S\) down (we got \(S=0.62\)), but a low \(S\) then predicts polarised loyalty for every brand, including the followers — which is exactly the over-prediction we see.
This shows up in the fit before duplication is ever computed. The per-brand \(S_j\) run 0.50, 0.68, 1.17, 1.11, 0.87 — ordered so that the leaders want a low \(S\) and the followers a high one. In the toothpaste example the \(S_j\) scatter is just as wide but unordered with respect to share. Systematic ordering of \(S_j\) by brand size is the early warning; the duplication residuals are the confirmation.
Terminology worth getting right. This is not “partitioning” in the usual sense of decaf-vs-regular or diet-vs-sugar, where a product attribute splits the category and within-partition duplication runs high. Here the leaders duplicate low. What is segmented is the consumer base, not the product space — two loyal camps rather than two sub-markets. The distinction matters for what you would do about it.
7.3 Cross-checking B1/B2 — is the incidence layer really innocent?
The verdict above rests on the NBD layer being well specified. That is worth testing rather than asserting, and there is a cheap check: fit \(K\) two independent ways. If the gamma-mixed Poisson is right they should agree.
n_by_hh = cat_occ.groupby("panel_id").size()
full = np.concatenate([n_by_hh.values, np.zeros(PANEL_SIZE - len(n_by_hh))])
K_moments = full.mean() ** 2 / (full.var(ddof=1) - full.mean())
display_markdown(
f"| estimator | $K$ |\n|---|---|\n"
f"| mean-and-zeros (used above) | {K:.4f} |\n"
f"| method of moments | {K_moments:.4f} |\n\n"
f"disagreement **{abs(K - K_moments) / K * 100:.0f}%** — "
f"observed variance {full.var(ddof=1):.2f} vs NBD-implied "
f"{full.mean() * (1 + full.mean() / K):.2f}",
raw=True,
)| estimator | \(K\) |
|---|---|
| mean-and-zeros (used above) | 2.6031 |
| method of moments | 2.5858 |
disagreement 1% — observed variance 10.14 vs NBD-implied 10.10
They agree to 1%. Two estimators using completely different features of the data — one the proportion of non-buyers, the other the variance — land in the same place. The incidence half of the model is fine. Whatever is wrong is on the brand-choice side, which is what the rest of the diagnosis says.
7.4 Where the Dirichlet breaks in general
Everything above is a claim about this data, tested against it. The list below is the broader catalogue of known failure modes, mapped onto the assumptions. None of these are diagnosed here — they are the checklist to run against the next category.
B1 — non-Poisson incidence. Strongly seasonal or episodic categories (sun care, antifreeze, Christmas goods); long-interpurchase-cycle durables (cars, appliances, mattresses), where replacement timing is a hazard process rather than a constant rate and you cannot observe enough repeat events anyway; categories with stockpiling, purchase acceleration or genuine satiation; and any analysis window short relative to the interpurchase interval — Ehrenberg’s “very short periods” boundary.
A1 — constrained choice sets. The model assumes every consumer could have chosen any brand. In practice incomplete distribution is the big one: a brand on shelf in 40% of stores is not being rejected by the other 60%, it is not on offer. This shows up as under-predicted penetration and is the standard explanation for deviations from Double Jeopardy. Same mechanism for regional brands, retailer-exclusive lines, and category-level restrictions (pet food for non-pet-owners, alcohol for abstainers).
The usual remedy is to fit to the category-user base rather than the population. Worth knowing what that costs:
M_buyers = n_by_hh.mean()
K_buyers = M_buyers ** 2 / (n_by_hh.var(ddof=1) - M_buyers)
display_markdown(
f"| population | $B$ | $M$ | $K$ | mean-and-zeros usable? |\n|---|---|---|---|---|\n"
f"| whole panel ({PANEL_SIZE:,}) | {cat_pen:.3f} | {M0:.3f} | {K:.3f} | yes |\n"
f"| category buyers ({len(n_by_hh):,}) | 1.000 | {M_buyers:.3f} | {K_buyers:.3f} "
f"| **no** — $\\log(1-B)=\\log 0$ |",
raw=True,
)| population | \(B\) | \(M\) | \(K\) | mean-and-zeros usable? |
|---|---|---|---|---|
| whole panel (5,021) | 0.911 | 3.989 | 2.603 | yes |
| category buyers (4,574) | 1.000 | 4.379 | 3.798 | no — \(\log(1-B)=\log 0\) |
Two consequences, both easy to trip over. Restricting to category buyers sets \(B=1\), which makes “mean and zeros” undefined — you must fall back to the method of moments. And \(K\) moves from 2.60 to 3.80, a 46% change, purely from redefining who counts as a potential buyer. Here the choice barely matters for the brand-level conclusions because penetration is already 91%; in a category bought by 20% of households it would dominate everything.
Contractual and subscription settings. Banking, insurance, telecoms, SaaS. Switching is discrete and costly, so repertoire buying is replaced by monogamy plus churn, and loyalty is structurally rather than stochastically high. This is the contractual column of the Fader–Hardie taxonomy, and it is the wrong tool entirely — not a misfit but a category error.
Small-\(N\), heavily concentrated B2B. Gamma heterogeneity across a handful of large accounts is a poor description, and \(S\) estimation becomes unstable — recall that \(S\) here is already a share-weighted average over per-brand estimates with outliers dropped, which needs enough brands to be meaningful.
7.5 What kind of product does this imply?
Reading the parameters as Ehrenberg would, as descriptors rather than fitted noise:
- \(K = 2.6\) (vs 0.78 for toothpaste) — a near-universal staple. 91% of households buy it, and they differ relatively little in how often.
- \(S = 0.62\) (vs 1.30) — polarised choice. Two-thirds of category buyers bought one brand all year.
- Effectively a duopoly with a fringe: 85% share in two brands that do not share buyers, plus three small brands functioning as occasional alternates.
Categories that behave this way tend to be ones where households settle on a staple and repurchase it without much deliberation, and where the two leaders are close functional substitutes differentiated by something habitual — taste, format, or household convention — rather than by occasion. The Dirichlet’s home ground is the opposite: a repertoire category where buyers rotate across several acceptable brands.
8 Alternatives within the Ehrenberg–Bass tradition
The tradition’s response to a Dirichlet misfit is rarely “use a fancier model”. It is usually to keep the model as the norm and treat the gap as the finding, or to apply a more targeted tool. Options, roughly in order of how well they suit this data:
8.1 Partitioned duplication coefficients — the direct fix here
Ehrenberg’s own remedy for exactly this situation (Repeat-Buying §10.5, Tables 10.8–10.9): abandon a single \(D\) and fit one per pair-type. In his worked example two manufacturer groups gave \(D = 3.5\) and \(2.5\) within, \(1.5\) across.
import itertools
group = {"Alpha": "leader", "Bravo": "leader",
"Charlie": "follower", "Delta": "follower", "Other": "follower"}
pair_key = lambda i, j: tuple(sorted((group[brand_name[i]], group[brand_name[j]])))
cells = {}
for key in [("leader", "leader"), ("follower", "leader"), ("follower", "follower")]:
num = den = 0.0
for i, j in itertools.combinations(range(nbrand), 2):
if pair_key(i, j) == key:
num += b_jk[i, j]
den += b_obs[i] * b_obs[j]
cells[key] = num / den
D_matrix = np.array([[cells[pair_key(i, j)] if i != j else 1.0
for j in range(nbrand)] for i in range(nbrand)])
r_single = (b_jk - D * np.outer(b_obs, b_obs))[off]
r_part = (b_jk - D_matrix * np.outer(b_obs, b_obs))[off]
display_markdown(
"| coefficient | value |\n|---|---|\n"
f"| single $D$, all pairs | {D:.3f} |\n"
f"| $D$ within leaders | {cells[('leader', 'leader')]:.3f} |\n"
f"| $D$ leader × follower | {cells[('follower', 'leader')]:.3f} |\n"
f"| $D$ within followers | {cells[('follower', 'follower')]:.3f} |\n\n"
f"mean abs. residual: single $D$ **{100 * np.abs(r_single).mean():.2f} pp** → "
f"partitioned **{100 * np.abs(r_part).mean():.2f} pp** "
f"(**{100 * (1 - np.abs(r_part).mean() / np.abs(r_single).mean()):.0f}%** reduction)",
raw=True,
)| coefficient | value |
|---|---|
| single \(D\), all pairs | 0.837 |
| \(D\) within leaders | 0.669 |
| \(D\) leader × follower | 0.950 |
| \(D\) within followers | 1.471 |
mean abs. residual: single \(D\) 1.02 pp → partitioned 0.40 pp (61% reduction)
Three numbers replace one and cut the mean absolute residual by 61%. They are also directly interpretable: the leaders avoid each other (\(D=0.67\)), leader-and-follower are near-independent (\(D=0.95\)), and the small brands cluster (\(D=1.47\) — heavy category buyers collecting alternates).
8.2 Other options, and when they apply
| Situation | Tool | Source |
|---|---|---|
| Non-stationarity | Conditional Trend Analysis — split buyers by prior-period purchase level and locate where the change came from, rather than fitting a dynamic model | Goodhardt & Ehrenberg (1967), JMR 4, 155–161; Repeat-Buying §7.6 |
| Non-stationarity, single brand | NBD alone for the focal brand — single-brand predictions are unaffected by other brands moving | Repeat-Buying §4.10 |
| Any misfit | Dirichlet-as-benchmark — fit stationary norms, report the gap as the result | 1984 paper §3.1 |
| Category incidence fits badly | Empirical Dirichlet — substitute the observed \(P_n\) for the fitted NBD. Better base-period fit, but you lose time extrapolation | 1984 paper §2.4 |
| Repeat-buying across periods | NBD/LSD — one-parameter formulae in \(q\) where the full Dirichlet has no closed form | Repeat-Buying Ch. 8 |
| Just need brand-vs-not | Beta-binomial conditional on category purchases — cheaper and often sufficient | Chatfield & Goodhardt (1970), Appl. Statist. 19, 240–250 |
| Heavy-buyer tail understated | Poisson-generalised inverse Gaussian — heavier-tailed than the gamma mixture | Sichel (1982), Appl. Statist. 31, 193–204 |
| Buying more regular than Poisson | Erlang interpurchase times | Chatfield & Goodhardt (1973), JASA 68, 828–835 |
| Nested choice (brand→flavour, brand→pack, brand→store) | Hierarchical Dirichlet — fit within partition, then across | Sketched in 1984 paper §5.3; store-choice in Kau & Ehrenberg (1984) |
For long-cycle durables, contractual and B2B categories, the tradition largely abandons the compound-Poisson machinery and works with the empirical generalisations directly — Double Jeopardy, Natural Monopoly, Duplication of Purchase, the Law of Buyer Moderation (heavy buyers regress toward the mean in the next period, light buyers upward), Pareto Share (the top 20% of buyers account for roughly 50–60% of volume, not 80%) — validated on those categories rather than modelled through the Dirichlet. That evidence sits in How Brands Grow Part 2 (Sharp & Romaniuk 2016) rather than in the 1984 paper.
Note what that implies methodologically: the laws outlive the model. This data set is a small instance of the same thing — the Dirichlet misfits, but Double Jeopardy holds in the observed columns regardless. Ehrenberg’s order of business was always empirical generalisation first, model second.
8.3 Adjacency to the rest of this repository
The Pareto/NBD and BG/NBD family in models/purchasing/ is the same Poisson-gamma skeleton with the brand-choice layer removed and a dropout process bolted on. The Dirichlet asks “which of several brands, in a stationary market”; BTYD asks “is this customer still alive, at one firm”.
The bridge is Assumption C. The Dirichlet insists purchase rate and brand preference are independent; a customer-base view is precisely where you would relax it. This data set gives a concrete reason to want to — the small brands are disproportionately bought by heavy category buyers, which is a \(\mu\)–\(\mathbf p\) dependence that C forbids.
9 Time extrapolation
Only \(M\) scales; \(K\) and \(S\) are structural.
yr2 = buy_table(M0 * 2, K, S, brand_share, brand_name)
pd.concat({"year 1": buy, "years 1-2": yr2}, axis=1).round(2)| year 1 | years 1-2 | |||||
|---|---|---|---|---|---|---|
| pen.brand | pur.brand | pur.cat | pen.brand | pur.brand | pur.cat | |
| Alpha | 0.54 | 3.27 | 4.75 | 0.64 | 5.54 | 8.75 |
| Bravo | 0.50 | 3.20 | 4.78 | 0.60 | 5.38 | 8.79 |
| Charlie | 0.14 | 2.70 | 5.01 | 0.17 | 4.25 | 9.19 |
| Delta | 0.06 | 2.63 | 5.05 | 0.08 | 4.09 | 9.26 |
| Other | 0.03 | 2.60 | 5.07 | 0.04 | 4.02 | 9.29 |
Code
tseq = np.arange(0, 2.01, 0.1)
pen_g = np.zeros((nbrand, len(tseq)))
buy_g = np.ones((nbrand, len(tseq)))
for i, t in enumerate(tseq):
if t == 0:
continue
tbl = buy_table(M0 * t, K, S, brand_share, brand_name)
pen_g[:, i] = tbl["pen.brand"].values
buy_g[:, i] = tbl["pur.brand"].values
fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
colors = plt.cm.viridis(np.linspace(0, 0.85, nbrand))
for j in range(nbrand):
axes[0].plot(tseq, pen_g[j], color=colors[j], lw=1.8, label=brand_name[j])
axes[1].plot(tseq, buy_g[j], color=colors[j], lw=1.8)
# observed year-1 penetration, for reference
axes[0].scatter(np.ones(nbrand), obs["b (O)"].values, color=colors,
zorder=5, s=28, marker="D", edgecolor="white", linewidth=0.8)
axes[0].set(xlabel="Years", ylabel="Penetration",
title="Penetration growth (♦ = observed year 1)")
axes[1].set(xlabel="Years", ylabel="Purchases per buyer", title="Buying-rate growth")
for ax in axes:
ax.grid(alpha=0.25, lw=0.5)
axes[0].legend(fontsize=8, frameon=False, loc="upper left")
fig.tight_layout()
plt.show()
Penetration grows less than pro rata; the extra buyers in a longer window are lighter buyers, and \(w\) absorbs the rest. The observed year-1 markers sit close to the curves for Alpha and Bravo and above them for the small brands — the same under-prediction seen in the O-vs-T table, now visible as a vertical gap.
10 What this application shows
Running the model on a category it was not chosen to flatter is more instructive than another clean reproduction:
- The parameters are readable. \(K = 2.6\) (vs 0.78) says a near-universally bought category; \(S = 0.62\) (vs 1.30) says a polarised one. Both match what the descriptive analysis found independently.
- Double Jeopardy survives. \(w(1-b)\) is far flatter than \(b\) or \(w\) across a twenty-fold share range, in the observed data as well as the fit. The empirical generalisation outlives the model that predicts it — which is Ehrenberg’s whole methodological point.
- The misfit is a finding, not a failure. Four independent measures agree that the small brands are repertoire add-ons rather than scaled-down leaders, and that Alpha and Bravo are near-mutually-exclusive. Assumption A2 is what breaks: the consumer base is segmented into two camps, which no single Dirichlet can represent.
- Measure duplication on the symmetric scale. Against the Dirichlet’s own \(b_{j\mid k}\) the anomaly looked broad; against \(D\,b_j b_k\) with observed penetrations it collapses to one dominant residual. Choosing the homoscedastic form changed the conclusion.
nstaris not a set-and-forget default. The package ships 50; this category needs ~170. The check is one line and belongs in any application.
The caveats from the main essay all still apply — stationarity, short periods, the heavy-buyer tail, the instability of \(S\) across period lengths. To those, this data set adds the multi-brand-trip wrinkle: 2.2% of occasions here allocate to more than one brand, which the model’s \(\sum_j r_j = n\) arithmetic does not admit.
The constructive close is that the tradition already has the tool for this: three partitioned duplication coefficients cut the residuals by 61% and say something usable about the market — the leaders avoid each other, the followers cluster. That is a better outcome than a better-fitting black box.