NBD-Dirichlet Model of Consumer Buying Behavior for Marketing Research

Author

Abdullah Mahmood

Published

July 19, 2026

Sources:

This essay implements the Dirichlet model from scratch in NumPy/pandas and reproduces the worked example from the NBDdirichlet package — UK toothpaste purchasing, first quarter of 1973, from the AGB panel of 5,240 continuously-reporting households. That is the same data set used in §3 of the 1984 paper.

Chapter-by-chapter reading notes on the underlying monographs are in references/papers/brand-choice/.

1 The question this model answers

Every other model in this repository looks at one firm’s customer base and asks who is alive and what they will be worth. The Dirichlet asks something different, of different data: given a whole product category and the brands competing in it, what should each brand’s loyalty measures look like, purely as a consequence of its size?

That framing matters. A brand manager sees that 53% of last quarter’s buyers came back this quarter and wants to know whether that is good. The number alone cannot answer it. What is needed is an interpretative norm — what a 19%-share brand in this category would show if nothing at all were special about it. The Dirichlet supplies that norm, and it does so from remarkably little input: category penetration, category buying rate, and each brand’s market share.

The inputs are panel aggregates, not per-customer RFM histories. No transaction log is required.

2 Imports

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)

3 The data

Four numbers describe the category, and two vectors describe the brands.

cat_pen = 0.56       # category penetration: 56% bought toothpaste in the quarter
cat_buyrate = 2.6    # category buyers bought, on average, 2.6 times in the quarter

brand_name = np.array(["Colgate DC", "Macleans", "Close Up", "Signal",
                       "ultrabrite", "Gibbs SR", "Boots Priv. Label",
                       "Sainsbury Priv. Lab."])
brand_share = np.array([0.25, 0.19, 0.10, 0.10, 0.09, 0.08, 0.03, 0.02])
brand_pen_obs = np.array([0.20, 0.17, 0.09, 0.08, 0.08, 0.07, 0.03, 0.02])

nbrand = len(brand_name)

pd.DataFrame({"Market share": brand_share, "Penetration (obs)": brand_pen_obs},
             index=brand_name)
Market share Penetration (obs)
Colgate DC 0.25 0.20
Macleans 0.19 0.17
Close Up 0.10 0.09
Signal 0.10 0.08
ultrabrite 0.09 0.08
Gibbs SR 0.08 0.07
Boots Priv. Label 0.03 0.03
Sainsbury Priv. Lab. 0.02 0.02

Market share here is the share of category purchase occasions, not of revenue. Working in purchase occasions rather than units or money is one of Ehrenberg’s four foundational analysis choices — it lets multi-unit purchases and multiple pack-sizes be handled by a single theory, and it is justified empirically (it produces generalisable results) rather than a priori.

Note already that share ranges over a factor of 12.5 (25% down to 2%) while penetration ranges over a factor of 10 (20% down to 2%). The two move together almost exactly. That is the first hint of the result the model will formalise.

4 The model

The Dirichlet is a mixture of distributions at four levels, two describing how often people buy the category and two describing which brand they pick. It is a counting + choice model in the taxonomy used elsewhere in this repository.

Level Assumption Aggregate consequence
B1 Each consumer’s category purchases follow a Poisson process with rate \(\mu\) —
B2 Rates \(\mu\) vary across consumers as a gamma distribution Category counts are NBD\((M, K)\)
A1 Each consumer’s brand choices are multinomial with probability vector \(\mathbf{p}\) —
A2 Choice vectors \(\mathbf{p}\) vary across consumers as a Dirichlet\((\alpha_1 \ldots \alpha_g)\) Brand counts given \(n\) are Dirichlet-multinomial
C Purchase rates and choice probabilities are independent across consumers The two halves multiply

Written as a compound distribution, the number of purchases of each of the \(g\) brands in a period of length \(T\) is

\[ \big[\mathcal{M}(\mathbf{r} \mid \mathbf{p}, n) \underset{\mathbf{p}}{\wedge} \mathcal{D}(\mathbf{p} \mid \boldsymbol{\alpha})\big] \underset{n}{\wedge} \big[\mathcal{P}(n \mid \mu) \underset{\mu}{\wedge} \mathcal{G}(\mu \mid MT, K)\big] \]

where \(\mathcal{M}, \mathcal{D}, \mathcal{P}, \mathcal{G}\) are the multinomial, Dirichlet, Poisson and gamma distributions.

Three structural parameters have to be estimated:

  • \(M\) — the mean category purchase rate per capita;
  • \(K\) — the NBD exponent, measuring how much consumers differ in how often they buy the category (smaller \(K\) ⇒ more diversity);
  • \(S = \sum_j \alpha_j\) — measuring how much consumers differ in which brand they buy (smaller \(S\) ⇒ more polarised, more sole-brand loyalty; larger \(S\) ⇒ everyone has much the same propensities).

Brand \(j\)’s share is \(\mu_j = \alpha_j / S\), so once \(S\) is known the individual \(\alpha_j = S\mu_j\) follow directly from the observed market shares.

4.1 Why the Dirichlet, specifically

The choice of the Dirichlet for A2 is not arbitrary convenience. There is a characterisation theorem behind it. If a market is unsegmented — meaning the proportion of purchases a consumer devotes to any one brand is independent of how they split the rest among the others — then the mixing distribution over choice vectors must be Dirichlet. It is the unique distribution expressing “independence except for the constraint \(\sum_j p_j = 1\)” (Mosimann 1962).

So the Dirichlet’s close fit across 40+ product categories is not merely a good curve fit; it is the quantitative evidence that these markets are largely unsegmented.

5 Layer 1 — category purchase incidence (the NBD)

5.1 From Poisson × gamma to the NBD

An individual buying at Poisson rate \(\mu\) over a period of length \(T\) makes \(n\) purchases with probability \(e^{-\mu T}(\mu T)^n / n!\). Mixing over a gamma distribution of \(\mu\) with shape \(K\) and mean \(M\) gives the negative binomial:

\[ P_n = \left(1 + \frac{MT}{K}\right)^{-K} \frac{\Gamma(K+n)}{n!\,\Gamma(K)} \left(\frac{MT}{MT+K}\right)^{n} \]

This is the same NBD as in nbd-overview — Poisson counting with gamma heterogeneity — applied here to the category rather than to one firm’s transactions.

5.2 Estimating \(M\)

\(M\) is the per-capita category purchase rate, which is simply penetration times the buyers’ rate:

\[ M = B \times W = 0.56 \times 2.6 \]

M0 = cat_pen * cat_buyrate
print(f"M = {M0:.4f} category purchases per capita per quarter")
M = 1.4560 category purchases per capita per quarter

5.3 Estimating \(K\) by “mean and zeros”

\(K\) is fitted by matching the model’s proportion of non-buyers to the observed one. Setting \(n = 0\) and \(T = 1\) in the NBD gives \(P_0 = (1 + M/K)^{-K}\), and \(P_0 = 1 - B\), so we need the \(K\) solving

\[ 1 - B = \left(1 + \frac{M}{K}\right)^{-K} \qquad\Longleftrightarrow\qquad K\log\!\left(1 + \frac{M}{K}\right) + \log(1 - B) = 0 \]

There is no closed form, so we minimise the squared residual numerically — exactly what the R package does with optimize.

This estimator is worth pausing on. It is not maximum likelihood on the frequency counts, which is what the Fader–Hardie essays in this repository use. It is Ehrenberg’s “mean and zeros” method, and it is deliberate: it is ≥90% efficient for typical consumer purchase data (against under 50% for the method of moments), and it requires only the two numbers that panel tabulations actually report — the mean and the proportion of non-buyers. You do not need the full frequency distribution.

The existence condition is \(M > -\log P_0\); below that no NBD can be fitted.

def estimate_K(M, cat_pen, max_K=30.0):
    """Fit the NBD exponent K by 'mean and zeros': solve (1 + M/K)^-K = 1 - B."""
    cp = np.log(1.0 - cat_pen)
    obj = lambda K: (K * np.log(1.0 + M / K) + cp) ** 2
    return minimize_scalar(obj, bounds=(1e-4, max_K), method="bounded").x

K = estimate_K(M0, cat_pen)
print(f"K = {K:.4f}")
print(f"check: implied P(0) = {(1 + M0 / K) ** -K:.4f}  vs observed {1 - cat_pen:.4f}")
K = 0.7788
check: implied P(0) = 0.4400  vs observed 0.4400

A \(K\) well below 1 signals a highly skewed, reverse-J purchase distribution — a large mass of light buyers and a thin tail of heavy ones. That is the norm for frequently-bought consumer goods.

5.4 The category purchase distribution \(P_n\)

We compute \(P_n\) in logs via a cumulative sum, which is both faster and numerically safer than forming the gamma ratios directly:

\[ \log P_n = -K\log\!\left(1+\frac{M}{K}\right) + \sum_{a=0}^{n-1}\big[\log(K+a) - \log(1+a)\big] + n\log\!\left(\frac{M}{M+K}\right) \]

NSTAR = 50  # truncation point for the infinite sums over category purchases

def Pn(M, K, n):
    """NBD probability of exactly n category purchases. Vectorised over n."""
    n = np.atleast_1d(np.asarray(n, dtype=int))
    nmax = int(n.max())
    a = np.arange(nmax)
    # cumulative log-ratio term; index 0 is the empty sum
    logratio = np.concatenate(([0.0], np.cumsum(np.log(K + a) - np.log(1.0 + a))))
    log_p = (-K * np.log(1.0 + M / K)
             + logratio[n]
             + np.where(n > 0, n * np.log(M / (M + K)), 0.0))
    return np.exp(log_p)

ns = np.arange(NSTAR + 1)
P = Pn(M0, K, ns)

pd.DataFrame({"n": ns[:9], "P(n)": P[:9]}).set_index("n").T
n 0 1 2 3 4 5 6 7 8
P(n) 0.44 0.22 0.13 0.08 0.05 0.03 0.02 0.01 7.50e-03

5.5 Checking the truncation

The theoretical sums run to infinity; we truncate at NSTAR. The 1984 paper uses a careful truncation-plus-lump procedure to preserve the mean exactly; the R package skips that for simplicity and just takes nstar large enough. We should verify that choice, because a too-small NSTAR silently biases every measure downward.

Two checks: the probabilities should sum to ≈1, and the implied mean should recover \(M\).

def check_truncation(M, K, nstar=NSTAR, tol_mean=0.1):
    ns = np.arange(nstar + 1)
    P = Pn(M, K, ns)
    total, mean = P.sum(), (ns * P).sum()
    ok = (total >= 0.99) and (abs(mean - M) <= tol_mean)
    return {"sum P(n)": total, "implied mean": mean, "target M": M, "OK": ok}

pd.DataFrame([check_truncation(M0 * t, K) | {"period t": t} for t in (1, 4, 8)]
             ).set_index("period t")
sum P(n) implied mean target M OK
period t
1 1.00 1.46 1.46 True
4 1.00 5.77 5.82 True
8 0.98 10.16 11.65 False

NSTAR = 50 is comfortable at one and four quarters. It would need raising for a much longer horizon or a more frequently-bought category — the diagnostic above is the way to tell.

6 Layer 2 — brand choice (Dirichlet-multinomial → beta-binomial)

6.1 The additivity property

The Dirichlet’s decisive technical property is additivity: any set of brands can be merged into a single “super-brand” whose parameter is the sum of theirs,

\[ \alpha_{(j+k)} = \alpha_j + \alpha_k \]

and nothing else about the model changes. This is not shared by competing models (e.g. Hendry) and it does two jobs here. It justifies lumping minor brands into an “all other brands” category, and — crucially for computation — it means that for a single brand \(j\) versus everything else, the \(g\)-dimensional problem collapses to just two categories.

A Dirichlet-multinomial with two categories is a beta-binomial. So conditional on \(n\) category purchases, the number \(r_j\) of purchases of brand \(j\) is

\[ p(r_j \mid n) = \binom{n}{r_j}\, \frac{B(\alpha_j + r_j,\; S - \alpha_j + n - r_j)}{B(\alpha_j,\; S - \alpha_j)} \]

with \(\alpha_j = S\mu_j\). Every brand performance measure below is a sum of the form \(\sum_n P_n \times (\text{something from this beta-binomial})\).

The special case \(r_j = 0\) telescopes into a product that avoids the beta function entirely, which is what the R package uses inside the \(S\) optimisation:

\[ p(0 \mid n) = \frac{(S-\alpha_j)(S-\alpha_j+1)\cdots(S-\alpha_j+n-1)}{S(S+1)\cdots(S+n-1)}, \qquad p(0 \mid 0) = 1 \]

def p_r_given_n(r, n, alpha, S):
    """Beta-binomial P(r purchases of the brand | n category purchases).

    Computed in log space for numerical stability at large n:
      log C(n,r) + log B(a+r, S-a+n-r) - log B(a, S-a)
    """
    r = np.asarray(r, dtype=float)
    return np.exp(
        gammaln(n + 1) - gammaln(r + 1) - gammaln(n - r + 1)
        + betaln(alpha + r, S - alpha + n - r)
        - betaln(alpha, S - alpha)
    )

def p_zero_given_n(n, alpha, S):
    """P(0 purchases of the brand | n category purchases), telescoped product form."""
    if n == 0:
        return 1.0
    a = np.arange(n)
    return np.exp(np.sum(np.log(S - alpha + a) - np.log(S + a)))

6.2 Estimating \(S\)

\(S\) is fitted from the observed brand penetrations. For a trial value of \(S\), the model’s penetration for brand \(j\) is

\[ b_j(S) = 1 - \sum_{n=0}^{n^*} P_n \, p(0 \mid n; \alpha_j = S\mu_j) \]

and we choose the \(S\) that reproduces the observed \(b_j\). Because this can be done brand by brand, we get \(g\) separate estimates \(S_1 \ldots S_g\) — which is itself useful: if the model fits, they should all be similar. Wide disagreement is a diagnostic that the category is segmented or otherwise ill-described.

def S_for_brand(j, M, K, share, pen_obs, nstar=NSTAR, max_S=30.0):
    """Solve for the S that reproduces brand j's observed penetration."""
    ns = np.arange(nstar + 1)
    P = Pn(M, K, ns)

    def obj(S):
        theo_pen = 1.0 - sum(P[n] * p_zero_given_n(n, S * share[j], S) for n in ns)
        return (theo_pen - pen_obs[j]) ** 2

    return minimize_scalar(obj, bounds=(1e-9, max_S), method="bounded").x

S_all = np.array([S_for_brand(j, M0, K, brand_share, brand_pen_obs)
                  for j in range(nbrand)])

pd.DataFrame({"Market share": brand_share, "Penetration (obs)": brand_pen_obs,
              "S_j": S_all.round(3)}, index=brand_name)
Market share Penetration (obs) S_j
Colgate DC 0.25 0.20 1.30
Macleans 0.19 0.17 2.21
Close Up 0.10 0.09 1.60
Signal 0.10 0.08 0.84
ultrabrite 0.09 0.08 1.45
Gibbs SR 0.08 0.07 1.29
Boots Priv. Label 0.03 0.03 2.22
Sainsbury Priv. Lab. 0.02 0.02 2.14

6.3 Pooling the per-brand estimates

The eight estimates range from 0.84 to 2.22 — a fair spread. Goodhardt, Ehrenberg & Chatfield’s own instruction is to drop discrepant brands and take a share-weighted average of the rest:

\[ \hat S = \frac{\sum_{j \in g^*} S_j \mu_j}{\sum_{j \in g^*} \mu_j} \]

Share-weighting gives more influence to large brands, whose penetrations are estimated from bigger sub-samples and are therefore more reliable.

The R package operationalises “discrepant” with a two-part boxplot rule:

  1. the usual 1.5 × IQR fence, and
  2. anything above the upper notch, \(\text{median} + 1.58 \times \text{IQR}/\sqrt{n}\).

Worth knowing if you are porting this: neither test appears in the package source. Its entire outlier step is four lines that call boxplot() and read two fields off the result —

bp <- boxplot(Sall, plot=F)
outlier  <- bp$out                    # 1.5 x IQR fence
outlier2 <- Sall[Sall > bp$conf[2]]   # above the upper notch

— so the actual arithmetic lives three packages away, in grDevices::boxplot.stats, which computes the fence and \(\text{conf} = \text{median} \pm 1.58\,\text{IQR}/\sqrt{n}\) from quartiles supplied by stats::fivenum. Reading dirichlet.R alone does not tell you what the rule is.

fivenum uses Tukey’s hinges rather than interpolated quartiles, so it is not np.percentile. We reproduce it directly below. (On this particular data the two definitions happen to flag the same three brands, so np.percentile would give the same \(S\) here — but the hinge version is what the package actually does, and there is no reason to expect the agreement to hold on other data.)

def fivenum(x):
    """Tukey's five-number summary, a direct port of R's stats::fivenum:

        n4 <- floor((n + 3)/2)/2
        d  <- c(1, n4, (n + 1)/2, n + 1 - n4, n)
        0.5 * (x[floor(d)] + x[ceiling(d)])

    Uses hinges, not linearly-interpolated quartiles, so the values feeding the
    fence and the notch match R's boxplot.stats() exactly.
    """
    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)])
    lo = np.floor(d).astype(int) - 1
    hi = np.ceil(d).astype(int) - 1
    return 0.5 * (x[lo] + x[hi])


def pool_S(S_all, share):
    """Share-weighted mean of per-brand S, dropping boxplot outliers and
    anything above the upper notch (matches R's NBDdirichlet)."""
    _, hinge_lo, median, hinge_hi, _ = fivenum(S_all)
    iqr = hinge_hi - hinge_lo
    outside_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 = outside_fence | (S_all > notch_hi)
    keep = ~drop
    return np.average(S_all[keep], weights=share[keep]), drop, notch_hi

S, dropped, notch_hi = pool_S(S_all, brand_share)

display_markdown(
    f"Upper notch = **{notch_hi:.4f}** → dropped: "
    f"**{', '.join(brand_name[dropped]) if dropped.any() else 'none'}**\n\n"
    f"Pooled **S = {S:.4f}**  (an unfiltered share-weighted mean would give "
    f"{np.average(S_all, weights=brand_share):.4f})",
    raw=True,
)

Upper notch = 2.0181 → dropped: Macleans, Boots Priv. Label, Sainsbury Priv. Lab.

Pooled S = 1.2953 (an unfiltered share-weighted mean would give 1.5499)

The notch rule removes Macleans and both private labels, and it moves \(S\) from 1.55 to 1.30 — a 16% difference that propagates into every measure below. It is not a cosmetic detail. (The 1984 paper reports \(S = 1.2\) for this data, fitted with the authors’ truncation procedure rather than a plain cutoff.)

6.4 The three estimated parameters

display_markdown(
    f"""
| Parameter | Value | Interpretation |
|---|---|---|
| $M$ | {M0:.3f} | category purchases per capita per quarter |
| $K$ | {K:.3f} | diversity in *how often* people buy the category |
| $S$ | {S:.3f} | diversity in *which brand* they buy |
""",
    raw=True,
)
Parameter Value Interpretation
\(M\) 1.456 category purchases per capita per quarter
\(K\) 0.779 diversity in how often people buy the category
\(S\) 1.295 diversity in which brand they buy

That is the whole model. Everything that follows uses only these three numbers plus the vector of market shares — the observed brand penetrations have done their job in fitting \(S\) and are no longer needed.

7 Brand performance measures

With \(\alpha_j = S\mu_j\) fixed, the three headline measures are:

\[ \begin{aligned} b_j &= 1 - \sum_{n\ge 0} P_n\, p(0 \mid n) && \text{penetration} \\[4pt] w_j &= \frac{1}{b_j}\sum_{n\ge 1} P_n \sum_{r=1}^{n} r\, p(r \mid n) && \text{purchases of the brand, per brand buyer} \\[4pt] w_{P,j} &= \frac{1}{b_j}\sum_{n\ge 1} n\, P_n \big[1 - p(0 \mid n)\big] && \text{purchases of the \emph{category}, per brand buyer} \end{aligned} \]

\(w_P\) is the one people find surprising, so it is worth stating plainly: it counts all category purchases made by someone who bought brand \(j\) at least once — including all the times they bought competitors.

def brand_measures(M, K, S, alpha, nstar=NSTAR, limit=None):
    """Return (penetration, E[brand purchases], E[category purchases]) for one brand.

    `limit` restricts the sum to a given set of category purchase frequencies,
    which is what the 'heavy buyer' tables need. The two expectations are
    returned *unnormalised* (i.e. not yet divided by penetration) so that the
    caller can divide by whichever base is appropriate.
    """
    ns = np.arange(nstar + 1) if limit is None else np.asarray(limit, dtype=int)
    P = Pn(M, K, ns)
    p0 = np.array([p_r_given_n(0, n, alpha, S) for n in ns])
    pen = 1.0 - np.sum(P * p0)

    pos = ns[ns >= 1]
    P_pos = Pn(M, K, pos)
    e_brand = np.array([np.sum(np.arange(1, n + 1)
                               * p_r_given_n(np.arange(1, n + 1), n, alpha, S))
                        for n in pos])
    p0_pos = np.array([p_r_given_n(0, n, alpha, S) for n in pos])

    return pen, np.sum(P_pos * e_brand), np.sum(pos * P_pos * (1.0 - p0_pos))


def buy_table(M, K, S, share, names, nstar=NSTAR):
    rows = []
    for j in range(len(share)):
        b, num_w, num_wp = brand_measures(M, K, S, S * share[j], nstar)
        rows.append([b, num_w / b, num_wp / 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
Colgate DC 0.20 1.82 3.16
Macleans 0.16 1.76 3.22
Close Up 0.09 1.68 3.30
Signal 0.09 1.68 3.30
ultrabrite 0.08 1.67 3.31
Gibbs SR 0.07 1.66 3.32
Boots Priv. Label 0.03 1.62 3.37
Sainsbury Priv. Lab. 0.02 1.61 3.38

7.1 Observed versus theoretical

Ehrenberg’s presentational discipline is to set observed against theoretical side by side, at two significant figures, so the reader can see where the model succeeds and fails rather than being handed a single fit statistic. The observed brand rate is \(w_j^{\text{obs}} = M\mu_j / b_j^{\text{obs}}\).

Code 
w_obs = M0 * brand_share / brand_pen_obs
ot = pd.DataFrame({
    "share": brand_share,
    "b (O)": brand_pen_obs, "b (T)": buy["pen.brand"].values,
    "w (O)": w_obs,         "w (T)": buy["pur.brand"].values,
    "wp (T)": buy["pur.cat"].values,
}, index=brand_name)
ot.loc["Average"] = ot.mean()
ot.round(2)
share b (O) b (T) w (O) w (T) wp (T)
Colgate DC 0.25 0.20 0.20 1.82 1.82 3.16
Macleans 0.19 0.17 0.16 1.63 1.76 3.22
Close Up 0.10 0.09 0.09 1.62 1.68 3.30
Signal 0.10 0.08 0.09 1.82 1.68 3.30
ultrabrite 0.09 0.08 0.08 1.64 1.67 3.31
Gibbs SR 0.08 0.07 0.07 1.66 1.66 3.32
Boots Priv. Label 0.03 0.03 0.03 1.46 1.62 3.37
Sainsbury Priv. Lab. 0.02 0.02 0.02 1.46 1.61 3.38
Average 0.11 0.09 0.09 1.64 1.69 3.30

Penetration is recovered to within about a point across all eight brands, from a 25%-share market leader down to a 2%-share private label — using market shares alone.

8 Double Jeopardy, read off the table

Look down the w (T) column. As share falls from 25% to 2%, the theoretical purchase frequency falls from 1.8 to 1.6 — it barely moves. Meanwhile penetration falls from 0.20 to 0.02, a factor of ten.

This is Double Jeopardy: small brands are punished twice, having both fewer buyers and slightly less loyal ones. The crucial point is that this is a prediction of the model, not an extra assumption fed into it. It falls out of the mathematics of share.

It is also derivable without the Dirichlet at all. From three separate empirical regularities — (A) buyers of any brand buy the category at about the same rate, (B) duplication follows \(b_{XY} = D\,b_X b_Y\) (implying near-zero correlation between buying X and buying Y), and (C) duplicated buyers buy a brand at about its normal rate — a few lines of algebra give

\[ w_X(1 - b_X) = w_Y(1 - b_Y) = \text{constant} \]

Let us check that our fitted model obeys it.

Code 
dj = pd.DataFrame({
    "share": brand_share,
    "b": buy["pen.brand"].values,
    "w": buy["pur.brand"].values,
    "w(1-b)": buy["pur.brand"].values * (1 - buy["pen.brand"].values),
}, index=brand_name)
dj.round(2)
share b w w(1-b)
Colgate DC 0.25 0.20 1.82 1.46
Macleans 0.19 0.16 1.76 1.49
Close Up 0.10 0.09 1.68 1.54
Signal 0.10 0.09 1.68 1.54
ultrabrite 0.09 0.08 1.67 1.54
Gibbs SR 0.08 0.07 1.66 1.55
Boots Priv. Label 0.03 0.03 1.62 1.58
Sainsbury Priv. Lab. 0.02 0.02 1.61 1.58

The w(1-b) column is essentially constant at 1.45–1.58 while share varies twelve-fold. The single factor \((1-b)\) does double duty: over long periods penetrations are high and \((1-b)\) changes fast, absorbing the strong trend in \(w\); over short periods penetrations are low and \((1-b) \approx 1\) for every brand, introducing no spurious trend. Analysis period length is handled indirectly, without appearing in the equation.

The pur.cat column carries the other half of the story. Buyers of the smallest brands buy the category most often (3.38 vs 3.16 for Colgate). Nobody is a devoted Sainsbury’s-private-label loyalist; its buyers are heavy toothpaste buyers who occasionally pick it up. This is the natural monopoly trend, and it too is a statistical selection effect predicted by the model rather than an assumption.

The marketing implication is the one Ehrenberg spent a career pressing, and it is where How Brands Grow comes from. Since \(m = b \times w\) and \(w\) barely varies, share differences are almost entirely penetration differences. Planning to grow by making existing buyers buy more is, in his phrase, “aiming at something altogether unusual or unlikely, like making pigs fly.”

9 The four summary tables

These reproduce the theoretical halves of Tables 3–6 of the 1984 paper, and match the summary.dirichlet method of the R package.

9.1 1. Purchase frequency distribution

How the brand’s buyers split by number of purchases: \(f(r) = \sum_{n \ge r} P_n \, p(r \mid n)\).

def freq_table(M, K, S, share, names, cutoff=5, nstar=NSTAR):
    ns = np.arange(nstar + 1)
    P = Pn(M, K, ns)

    def prob_r(r, alpha):
        return np.sum([P[n] * p_r_given_n(r, n, alpha, S) for n in range(r, nstar + 1)])

    rows = []
    for j in range(len(share)):
        a = S * share[j]
        head = [prob_r(r, a) for r in range(cutoff + 1)]
        tail = np.sum([prob_r(r, a) for r in range(cutoff + 1, nstar + 1)])
        rows.append(head + [tail])

    cols = [str(i) for i in range(cutoff + 1)] + [f"{cutoff + 1}+"]
    return pd.DataFrame(rows, index=names, columns=cols)

freq_table(M0, K, S, brand_share, brand_name).round(2)
0 1 2 3 4 5 6+
Colgate DC 0.80 0.12 0.04 0.02 0.01 0.0 0.01
Macleans 0.84 0.10 0.03 0.01 0.01 0.0 0.00
Close Up 0.91 0.06 0.02 0.01 0.00 0.0 0.00
Signal 0.91 0.06 0.02 0.01 0.00 0.0 0.00
ultrabrite 0.92 0.05 0.02 0.01 0.00 0.0 0.00
Gibbs SR 0.93 0.05 0.01 0.01 0.00 0.0 0.00
Boots Priv. Label 0.97 0.02 0.01 0.00 0.00 0.0 0.00
Sainsbury Priv. Lab. 0.98 0.01 0.00 0.00 0.00 0.0 0.00

Read row 1: of the whole population, 80% buy no Colgate in the quarter, 12% buy it once, 4% twice. The reverse-J shape is universal — most buyers of any brand are light buyers of it.

9.2 2. Heavy versus light category buyers

Restricting the sums to a set \(R\) of category purchase frequencies gives penetration and rate among consumers whose category buying falls in that range:

\[ b_{j \mid R} = 1 - \frac{\sum_{n \in R} P_n\, p(0 \mid n)}{\sum_{n \in R} P_n}, \qquad w_{j \mid R} = \frac{\sum_{n \in R} P_n \sum_{r=1}^{n} r\,p(r\mid n)}{\sum_{n \in R} P_n\,[1 - p(0 \mid n)]} \]

This is a fairly direct test of assumption (C) — that brand choice is independent of category purchase rate.

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, num_w, _ = 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, num_w / (P_sum - not_buying)])
    return pd.DataFrame(rows, index=names,
                        columns=["Penetration", "Avg Purchase Freq"])

heavy_table(M0, K, S, brand_share, brand_name, limit=range(1, 7)).round(2)
Penetration Avg Purchase Freq
Colgate DC 0.34 1.61
Macleans 0.27 1.57
Close Up 0.15 1.51
Signal 0.15 1.51
ultrabrite 0.13 1.50
Gibbs SR 0.12 1.49
Boots Priv. Label 0.05 1.46
Sainsbury Priv. Lab. 0.03 1.45

9.3 3. Brand duplication

The proportion of brand \(k\)’s buyers who also bought brand \(j\). Additivity again does the work — form the composite brand \((j+k)\), get its penetration, and recover the overlap by inclusion–exclusion:

\[ b_{(j+k)} = 1 - \sum_n P_n\, p_j(0\mid n)\, p_k(0 \mid n), \qquad b_{jk} = b_j + b_k - b_{(j+k)}, \qquad b_{j \mid k} = \frac{b_{jk}}{b_k} \]

def dup_table(M, K, S, share, names, focal=0, nstar=NSTAR):
    ns = np.arange(nstar + 1)
    P = Pn(M, K, ns)
    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
        alpha_comp = S * (share[focal] + share[j])   # additivity: composite brand
        b_comp = 1.0 - np.sum([P[n] * p_r_given_n(0, n, alpha_comp, S) for n in ns])
        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=f"also buy | bought {names[focal]}")

dup = pd.concat([dup_table(M0, K, S, brand_share, brand_name, focal=f)
                 for f in (0, 1)], axis=1)
dup.round(2)
also buy | bought Colgate DC also buy | bought Macleans
Colgate DC 1.00 0.24
Macleans 0.19 1.00
Close Up 0.10 0.10
Signal 0.10 0.10
ultrabrite 0.09 0.09
Gibbs SR 0.08 0.08
Boots Priv. Label 0.03 0.03
Sainsbury Priv. Lab. 0.02 0.02

Now compare each column against the penetration vector. Colgate’s buyers buy Macleans at 0.19; Macleans’ buyers buy Colgate at 0.24. Against penetrations of 0.16 and 0.20 respectively, both are duplication ratios of about 1.2.

That is the Duplication of Purchase Law: \(b_{j\mid k} \approx D \cdot b_j\), with a single \(D\) for all brand pairs. Brands share customers in proportion to penetration, not according to positioning, price tier or manufacturer. Colgate does not have a special rivalry with Macleans; it shares buyers with every brand in proportion to how big that brand is.

Code 
D = dup.iloc[:, 0].values / buy["pen.brand"].values
pd.DataFrame({"b_j|Colgate": dup.iloc[:, 0].values,
              "b_j (penetration)": buy["pen.brand"].values,
              "ratio D": D}, index=brand_name).drop(index="Colgate DC").round(2)
b_j|Colgate b_j (penetration) ratio D
Macleans 0.19 0.16 1.19
Close Up 0.10 0.09 1.19
Signal 0.10 0.09 1.19
ultrabrite 0.09 0.08 1.19
Gibbs SR 0.08 0.07 1.19
Boots Priv. Label 0.03 0.03 1.19
Sainsbury Priv. Lab. 0.02 0.02 1.19

10 Time extrapolation: only \(M\) scales

For a period \(T\) times the base period, every formula is unchanged except that \(M\) becomes \(MT\). \(K\) and \(S\) are structural constants describing how consumers differ from one another, and those differences do not depend on how long you observe them.

This is the same \(k\)-invariance that underpins NBD penetration-growth prediction in the single-brand case, and it is empirically verified over 4/8/12/24-week periods in Repeat-Buying Ch. 7. It is what lets a single quarter of panel data predict the annual picture.

def set_period(M0, t):
    """Only M scales with the analysis period; K and S are time-invariant."""
    return M0 * t

annual = buy_table(set_period(M0, 4), K, S, brand_share, brand_name)
compare = pd.concat([buy.add_suffix(" (Q)"), annual.add_suffix(" (Yr)")], axis=1)
compare[["pen.brand (Q)", "pen.brand (Yr)", "pur.brand (Q)", "pur.brand (Yr)",
         "pur.cat (Q)", "pur.cat (Yr)"]].round(2)
pen.brand (Q) pen.brand (Yr) pur.brand (Q) pur.brand (Yr) pur.cat (Q) pur.cat (Yr)
Colgate DC 0.20 0.39 1.82 3.74 3.16 8.78
Macleans 0.16 0.31 1.76 3.51 3.22 8.98
Close Up 0.09 0.18 1.68 3.18 3.30 9.29
Signal 0.09 0.18 1.68 3.18 3.30 9.29
ultrabrite 0.08 0.17 1.67 3.15 3.31 9.32
Gibbs SR 0.07 0.15 1.66 3.11 3.32 9.35
Boots Priv. Label 0.03 0.06 1.62 2.91 3.37 9.46
Sainsbury Priv. Lab. 0.02 0.04 1.61 2.86 3.38 9.43

Penetration grows less than pro rata — Colgate goes from 20% to 37%, not to 80% — because the extra buyers picked up in a longer window are progressively lighter buyers. Purchase frequency absorbs the rest of the growth.

Note too that the Double Jeopardy spread widens over the year: annual \(w\) runs from 3.8 down to 1.9, a much bigger spread than the quarterly 1.8 to 1.6. This is exactly what \(w = c/(1-b)\) predicts once penetrations get large enough for \((1-b)\) to bite.

Code 
tseq = np.arange(0, 8.01, 0.25)
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(set_period(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.5))
colors = plt.cm.viridis(np.linspace(0, 0.9, 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)

axes[0].set(xlabel="Quarters", ylabel="Penetration",
            title="Theoretical penetration growth")
axes[1].set(xlabel="Quarters", ylabel="Purchases per buyer",
            title="Theoretical buying-rate growth")
for ax in axes:
    ax.grid(alpha=0.25, lw=0.5)
axes[0].legend(fontsize=7, frameon=False, loc="upper left")
fig.tight_layout()
plt.show()

Theoretical penetration and purchase-frequency growth over eight quarters.

The curves never cross. A brand that starts smaller stays smaller on both measures, at every horizon — Double Jeopardy holding across time as well as across brands.

11 Validation against the R package

Every table above was checked against NBDdirichlet v1.4 run on the same inputs. The agreement is exact at the two decimal places the package reports.

r_buy = pd.DataFrame(
    {"pen.brand": [0.20, 0.16, 0.09, 0.09, 0.08, 0.07, 0.03, 0.02],
     "pur.brand": [1.82, 1.76, 1.68, 1.68, 1.67, 1.66, 1.62, 1.61],
     "pur.cat":   [3.16, 3.22, 3.30, 3.30, 3.31, 3.32, 3.37, 3.38]},
    index=brand_name)

diff = (buy.round(2) - r_buy).abs()
display_markdown(
    f"Parameters — R: `M = 1.46, K = 0.78, S = 1.3` · "
    f"here: `M = {M0:.2f}, K = {K:.2f}, S = {S:.2f}`\n\n"
    f"Max absolute difference across the `buy` table: **{diff.values.max():.4f}**",
    raw=True,
)

Parameters — R: M = 1.46, K = 0.78, S = 1.3 · here: M = 1.46, K = 0.78, S = 1.30

Max absolute difference across the buy table: 0.0000

Two notes on reproducing this package, both of which cost real time to track down:

  • The S outlier rule is two-part. Implementing only the 1.5 × IQR fence — which is what a natural reading of “boxplot outliers” suggests — removes nothing from this data set and yields \(S = 1.55\). The upper-notch test is what removes three brands and gives \(S = 1.30\). An existing Python port of this package (references/implementations/NBDdirichlet-main/) omits the notch and so reports \(S = 1.55\), with knock-on errors of ~0.05 in every purchase-rate figure.
  • The package’s own docs disagree with the package. man/print.dirichlet.Rd shows S = 1.55 in its example output, but running v1.4 gives S = 1.3, matching the vignette. The Rd example block is stale. Running the code is the only reliable arbiter, which is why the numbers above were checked against a live R session rather than against the documentation.

12 Where this model breaks down

Ehrenberg is scrupulous about the boundaries, and an essay that reported only the successes would misrepresent him. Carry these:

  • Stationarity is required. The model describes a market in equilibrium — no trend in aggregate sales. It does not explain change, and the authors are explicit that “nobody yet knows much about change, so we cannot try to explain it.” Its use in non-stationary settings is as the counterfactual: what would have happened without the promotion, the seasonal peak, the launch.
  • Very short periods fail. At or near a product’s minimum inter-purchase interval (about a week for most groceries) behaviour is qualitatively different, dominated by shopping-trip and “dead period” effects rather than brand choice.
  • The heavy-buyer tail is over-predicted — the “variance discrepancy” \(\sigma - s \approx m\), whose real cause is shelving: purchases bunch at the number of weeks in the period, and there are consistently too few buyers above that.
  • Saturated categories with no direct substitute (toothpaste is itself an example) show a small shortfall of once-only category buyers, so the NBD layer fits less well than the brand layer. Where this happens, the authors substitute the observed \(P_n\) for the fitted NBD and call it the Empirical Dirichlet — a better base-period fit, at the cost of losing the ability to extrapolate across time.
  • \(S\) is not as stable as the theory requires. It should be invariant to period length; in the 1984 paper’s own data it moved from 1.0 (4-week) to 2.2 (12-week) to 1.8 (annual). The authors flag this as unresolved.

The honest summary is Ehrenberg’s own: the justification for the model is not that it is true, but that it works — over 40-plus product fields, two continents, and three decades — and that where it fails, it fails systematically enough to be informative.