CLV Quickstart

Customer Lifetime Value (CLV) is the measure of a customer’s contribution over time to a business. This metric is used to inform spending levels on new customer acquisition, retention, and other marketing and sales efforts, so reliable estimation is essential.

PyMC-Marketing provides tools to segment customers on their past behavior (see RFM Segmentation) as well as the following Buy Till You Die (BTYD) probabilistic models to predict future behavior:

• BG/NBD model for continuous time, non-contractual modeling

• Pareto/NBD model for continuous time, non-contractual modeling with covariates

• Shifted BG model for discrete time, contractual modeling with cohorts and covariates

• BG/BB model for discrete time, contractual modeling

• Gamma-Gamma model for expected monetary value

• Modified BG/NBD model, similar to the BG/NBD model, but assumes non-repeat customers are still active.

This table contains a breakdown of the four BTYD modeling domains, and examples for each:

	Non-contractual	Contractual
Continuous	online purchases	ad conversion time
Discrete	concerts & sports events	recurring subscriptions

In this notebook we will demonstrate how to estimate future purchasing activity and CLV with the CDNOW dataset, a popular benchmarking dataset in CLV and BTYD research. Data is available here, with additional details here.

import arviz as az
import arviz_plots as azp
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
from arviz_base.labels import MapLabeller

from pymc_marketing import clv

az.style.use("arviz-darkgrid")

%config InlineBackend.figure_format = "retina" # nice looking plots

1.1 Data Requirements

For all models, the following nomenclature is used:

• customer_id represents a unique identifier for each customer.

• frequency represents the number of repeat purchases that a customer has made, i.e. one less than the total number of purchases.

• T represents a customer’s “age”, i.e. the duration between a customer’s first purchase and the end of the period of study. In this example notebook, the units of time are in weeks.

• recency represents the time period when a customer made their most recent purchase. This is equal to the duration between a customer’s first and last purchase. If a customer has made only 1 purchase, their recency is 0.

• monetary_value represents the average value of a given customer’s repeat purchases. Customers who have only made a single purchase have monetary values of zero.

The rfm_summary function can be used to preprocess raw transaction data for modeling:

raw_trans = pd.read_csv(
    "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/cdnow_transactions.csv"
)

raw_trans.head(5)

	_id	id	date	cds_bought	spent
0	4	1	19970101	2	29.33
1	4	1	19970118	2	29.73
2	4	1	19970802	1	14.96
3	4	1	19971212	2	26.48
4	21	2	19970101	3	63.34

rfm_data = clv.utils.rfm_summary(
    raw_trans,
    customer_id_col="id",
    datetime_col="date",
    monetary_value_col="spent",
    datetime_format="%Y%m%d",
    time_unit="W",
)

rfm_data

	customer_id	frequency	recency	T	monetary_value
0	1	3.0	49.0	78.0	23.723333
1	2	1.0	2.0	78.0	11.770000
2	3	0.0	0.0	78.0	0.000000
3	4	0.0	0.0	78.0	0.000000
4	5	0.0	0.0	78.0	0.000000
...	...	...	...	...	...
2352	2353	2.0	53.0	66.0	19.775000
2353	2354	5.0	24.0	66.0	44.928000
2354	2355	1.0	44.0	66.0	24.600000
2355	2356	6.0	62.0	66.0	31.871667
2356	2357	0.0	0.0	66.0	0.000000

2357 rows × 5 columns

It is important to note these definitions differ from that used in RFM segmentation, where the first purchase is included, T is not used, and recency is the number of time periods since a customer’s most recent purchase.

To visualize data in RFM format, we can plot the recency and T of the customers with the plot_customer_exposure function. We see a large chunk (>60%) of customers haven’t made another purchase in a while.

fig, ax = plt.subplots(figsize=(10, 5))
(
    rfm_data.sample(n=100, random_state=42)
    .sort_values(["recency", "T"])
    .pipe(clv.plot_customer_exposure, ax=ax, linewidth=0.5, size=0.75)
);

Predicting Future Purchasing Behavior with the BG/NBD Model

This dataset is an example of continuous time, non-contractual transactions because it comprises purchases from an online music store. PyMC-Marketing provides several models for this use case:

• Beta-Geometric/Negative Binomial Distribution (BG/NBD)

• Pareto/Negative Binomial Distribution (Pareto/NBD)

• Modified Beta-Geometric/Negative Binomial Distribution (MBG/NBD)

We will be using the BG/NBD model in this notebook because it works well for basic use cases. For more comprehensive modeling, the Pareto/NBD and MBG/NBD models have expanded functionality and fewer limitations.

Let’s create a BG/NBD model with the default model configuration:

bgm = clv.BetaGeoModel()
bgm.build_model(data=rfm_data)
bgm

BG/NBD
            alpha ~ Weibull(2, 10)
      phi_dropout ~ Uniform(0, 1)
    kappa_dropout ~ Pareto(1, 1)
                r ~ Weibull(2, 1)
                a = Deterministic(f(kappa_dropout, phi_dropout))
                b = Deterministic(f(kappa_dropout, phi_dropout))
recency_frequency ~ BetaGeoNBD(a, b, r, alpha, <constant>)

All parameters for the BG/NBD model must be continuous and positive-valued. Weibull distributions work very well for parameters with these requirements, and are used as the defaults for the alpha and r parameters.

The a and b dropout parameters are pooled with the hierarchical kappa_dropout and phi_dropout priors for improved performance. This heirarchical structure can be displayed visually:

bgm.graphviz()

Model Fitting with MAP

By default, fitted models generates full Bayesian posteriors via MCMC sampling. However, for extremely large datasets where uncertainty estimates are not needed and/or MCMC is too slow, Maximum a Posteriori can be used to estimate point estimates for model parameters.

Use rfm_train_test_split for train/test splits. Here we are training on 52 weeks of data, and withholding the remaining 26 weeks for testing:

rfm_train_test_data = clv.utils.rfm_train_test_split(
    raw_trans,
    customer_id_col="id",
    datetime_col="date",
    monetary_value_col="spent",
    train_period_end="19980101",
    datetime_format="%Y%m%d",
    time_unit="W",
)

rfm_train_test_data

	customer_id	frequency	recency	T	monetary_value	test_frequency	test_monetary_value	test_T
0	1	3.0	49.0	52.0	23.723333	0.0	0.000	26.0
1	2	1.0	2.0	52.0	11.770000	0.0	0.000	26.0
2	3	0.0	0.0	52.0	0.000000	0.0	0.000	26.0
3	4	0.0	0.0	52.0	0.000000	0.0	0.000	26.0
4	5	0.0	0.0	52.0	0.000000	0.0	0.000	26.0
...	...	...	...	...	...	...	...	...
2352	2353	0.0	0.0	40.0	0.000000	2.0	19.775	26.0
2353	2354	5.0	24.0	40.0	44.928000	0.0	0.000	26.0
2354	2355	0.0	0.0	40.0	0.000000	1.0	24.600	26.0
2355	2356	4.0	26.0	40.0	33.317500	2.0	28.980	26.0
2356	2357	0.0	0.0	40.0	0.000000	0.0	0.000	26.0

2357 rows × 8 columns

Note any customers who made their first purchase during the testing period will be excluded. Use rfm_summary to retain all customers for the final model.

bgm_map = clv.BetaGeoModel()

bgm_map.fit(
    data=rfm_train_test_data,
    method="map",
)
bgm_map.fit_summary()

MAP ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━   0% 0:00:02 logp = -11,946, ||grad|| = 3.026

alpha             6.8
phi_dropout      0.22
kappa_dropout     2.2
r                0.28
a                0.47
b                 1.7
Name: value, dtype: object

Let’s evaluate model performance by tracking predictions against historical purchases.

clv.plot_expected_purchases_over_time(
    model=bgm_map,
    purchase_history=raw_trans,
    datetime_col="date",
    customer_id_col="id",
    datetime_format="%Y%m%d",
    time_unit="W",
    t=78,
    set_index_date=True,
    t_start_eval=52,
    plot_cumulative=False,
);

Customers purchasing behavior in this dataset shows a strong trend changepoint in March 1997, which the model was able to capture!

Model Fitting with MCMC

Fitting the same model with MCMC on the same data and the same priors lets us check the point estimate is consistent, and gives us the full posterior needed for credibility intervals on parameters and downstream predictions.

MCMC sampling produces full posteriors, which give credibility intervals around parameter estimates and predictions. Here we fit a second model with MCMC on the same rfm_train_test_data split and the same model_config used for MAP, so the comparison isolates the inference method.

bgm_mcmc = clv.BetaGeoModel()

bgm_mcmc.fit(data=rfm_train_test_data)
bgm_mcmc.fit_summary()

Computing ... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
alpha	6.85	0.53	6	7.7	1868	2215	1.00	0.012	0.0087
phi_dropout	0.2147	0.0191	0.18	0.25	1972	2284	1.00	0.00043	0.00031
kappa_dropout	2.48	0.75	1.6	3.8	2082	2170	1.00	0.016	0.023
r	0.2848	0.014	0.26	0.31	1811	1936	1.00	0.00033	0.00023
a	0.524	0.127	0.36	0.75	2683	2446	1.00	0.0024	0.0026
b	1.96	0.63	1.2	3.1	2002	2032	1.00	0.014	0.021

We can use ArviZ, a Python library tailored to produce visualizations for Bayesian models, to plot the posterior distribution of each parameter.

azp.plot_dist(bgm_mcmc.fit_result);

Mean values of the fitted posterior distributions strongly align with MAP point estimates. For the remainder of the notebook we will fit the BG/NBD model on the full RFM dataset so downstream predictions and CLV estimates use all available information.

bgm.fit(data=rfm_data)
bgm.fit_summary()

Computing ... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
alpha	7.11	0.51	6.3	7.9	2072	2087	1.00	0.011	0.0078
phi_dropout	0.1843	0.0151	0.16	0.21	1910	2430	1.00	0.00034	0.00025
kappa_dropout	3.25	0.86	2.1	4.8	1836	2000	1.00	0.02	0.023
r	0.277	0.0123	0.26	0.3	2122	1967	1.00	0.00027	0.00019
a	0.591	0.123	0.42	0.8	2261	2171	1.00	0.0026	0.0027
b	2.66	0.74	1.7	4	1790	1990	1.00	0.018	0.02

Visualizing Prediction Matrices

clv.plot_frequency_recency_matrix(bgm);

We can see our best customers have been active for over 60 weeks and have made over 20 purchases (bottom-right). Note the “tail” sweeping up towards the upper-left corner - these customers are infrequent and/or may not have purchased recently. What is the probability they are still active?

clv.plot_probability_alive_matrix(bgm)

<Axes: title={'center': 'Probability Customer is Alive,\nby Frequency and Recency of a Customer'}, xlabel="Customer's Historical Frequency", ylabel="Customer's Recency">

Note that all non-repeat customers have an alive probability of 1, which is one of the quirks of BetaGeoModel. In many use cases this is still a valid assumption, but if non-repeat customers are a key focus in your use case, you may want to try ParetoNBDModel instead.

Looking at the probability alive matrix, we can infer that customers who have made fewer purchases are less likely to return, and may be worth targeting for retention.

Ranking customers from best to worst

Having fit the model, we can ask what is the expected number of purchases for our customers over the next 10 time periods. Let’s look at the four most promising customers.

num_purchases = bgm.expected_purchases(future_t=10)

sdata = rfm_data.copy()
sdata["expected_purchases"] = num_purchases.mean(("chain", "draw")).values
sdata.sort_values(by="expected_purchases").tail(4)

	customer_id	frequency	recency	T	monetary_value	expected_purchases
812	813	30.0	72.0	74.0	35.654000	3.475635
1202	1203	32.0	71.0	72.0	47.172187	3.844154
156	157	36.0	74.0	77.0	30.603611	3.948038
1980	1981	35.0	66.0	68.0	46.748857	4.352503

We can also plot credibility intervals for the expected purchases:

ids = [813, 1203, 157, 1981]
pc = azp.plot_dist(
    xr.Dataset({"num_purchases": num_purchases.sel(customer_id=ids)}),
    figure_kwargs={"figsize": (14, 4)},
)
for axi, id in zip(pc.viz["plot"]["num_purchases"].values, ids, strict=True):
    axi.set_title(f"Customer: {id}", size=14)
pc.viz["figure"].item().suptitle(
    "Expected Number of Purchase over 10 Time Periods",
    fontsize=18,
    fontweight="bold",
    y=1.15,
);

Predicting purchase behavior of a new customer

We can use the fitted model to predict the number of purchases for a fresh new customer.

azp.plot_dist(
    xr.Dataset(
        {"num_purchases": bgm.expected_purchases_new_customer(t=10).sel(customer_id=1)}
    )
)
plt.title("Expected purchases of a new customer in the first 10 periods");

Customer Probability Histories

Given a customer transaction history, we can calculate their historical probability of being alive, according to our trained model.

Let’s look at active customer 1516 and assess the change in probability that the user will ever return if they do no other purchases in the next 9 time periods.

customer_1516 = rfm_data.loc[1515]
customer_1516

customer_id       1516.000000
frequency           27.000000
recency             67.000000
T                   70.000000
monetary_value      51.944074
Name: 1515, dtype: float64

customer_1516_history = pd.DataFrame(
    dict(
        customer_id=np.arange(10),
        frequency=np.full(10, customer_1516["frequency"], dtype="int"),
        recency=np.full(10, customer_1516["recency"]),
        T=(np.arange(0, 10) + customer_1516["recency"]).astype("int"),
    )
)
customer_1516_history

	customer_id	frequency	recency	T
0	0	27	67.0	67
1	1	27	67.0	68
2	2	27	67.0	69
3	3	27	67.0	70
4	4	27	67.0	71
5	5	27	67.0	72
6	6	27	67.0	73
7	7	27	67.0	74
8	8	27	67.0	75
9	9	27	67.0	76

p_alive = bgm.expected_probability_alive(data=customer_1516_history)

hdi = az.hdi(p_alive, prob=0.94)
plt.fill_between(
    customer_1516_history["T"],
    hdi.sel(ci_bound="lower"),
    hdi.sel(ci_bound="upper"),
    color="C0",
    alpha=0.3,
)
plt.plot(customer_1516_history["T"], p_alive.mean(("draw", "chain")), marker="o")
plt.axvline(
    customer_1516_history["recency"].iloc[0], c="black", ls="--", label="Purchase"
)

plt.title("Probability Customer 1516 will purchase again")
plt.xlabel("T")
plt.ylabel("p")
plt.legend();

We can see that, if no purchases are being made in the next 9 weeks, the model has low confidence that the customer will ever return. What if they had done one purchase in between?

customer_1516_history.loc[7:, "frequency"] += 1
customer_1516_history.loc[7:, "recency"] = customer_1516_history.loc[7, "T"] - 0.5
customer_1516_history

	customer_id	frequency	recency	T
0	0	27	67.0	67
1	1	27	67.0	68
2	2	27	67.0	69
3	3	27	67.0	70
4	4	27	67.0	71
5	5	27	67.0	72
6	6	27	67.0	73
7	7	28	73.5	74
8	8	28	73.5	75
9	9	28	73.5	76

p_alive = bgm.expected_probability_alive(data=customer_1516_history)

hdi = az.hdi(p_alive, prob=0.94)
plt.fill_between(
    customer_1516_history["T"],
    hdi.sel(ci_bound="lower"),
    hdi.sel(ci_bound="upper"),
    color="C0",
    alpha=0.3,
)
plt.plot(customer_1516_history["T"], p_alive.mean(("draw", "chain")), marker="o")
plt.axvline(
    customer_1516_history["recency"].iloc[0], c="black", ls="--", label="Purchase"
)
plt.axvline(customer_1516_history["recency"].iloc[-1], c="black", ls="--")

plt.title("Probability Customer 1516 will purchase again")
plt.xlabel("T")
plt.ylabel("p")
plt.legend();

From the plot above, say that customer 1516 makes a purchase at week 73.5, just over 6 weeks after they have made their last recorded purchase. We can see that the probability of the customer returning quickly goes back up!

Estimating Customer Lifetime Value Using the Gamma-Gamma Model

Until now we’ve focused mainly on transaction frequencies and probabilities, but to estimate economic value we can use the Gamma-Gamma model.

The Gamma-Gamma model assumes at least 1 repeat transaction has been observed per customer. As such we filter out those with zero repeat purchases.

nonzero_data = rfm_data.query("frequency>0")
nonzero_data

	customer_id	frequency	recency	T	monetary_value
0	1	3.0	49.0	78.0	23.723333
1	2	1.0	2.0	78.0	11.770000
5	6	14.0	76.0	78.0	76.503571
6	7	1.0	5.0	78.0	11.770000
7	8	1.0	61.0	78.0	26.760000
...	...	...	...	...	...
2351	2352	1.0	47.0	66.0	14.490000
2352	2353	2.0	53.0	66.0	19.775000
2353	2354	5.0	24.0	66.0	44.928000
2354	2355	1.0	44.0	66.0	24.600000
2355	2356	6.0	62.0	66.0	31.871667

1126 rows × 5 columns

If computing the monetary value from your own data, note that it is the mean of a given customer’s value, not the sum. monetary_value can be used to represent profit, or revenue, or any value as long as it is consistently calculated for each customer.

The Gamma-Gamma model relies on the important assumption that there is no relationship between the monetary value and the purchase frequency. In practice we need to check whether the Pearson correlation is less than 0.3:

nonzero_data[["monetary_value", "frequency"]].corr()

	monetary_value	frequency
monetary_value	1.000000	0.052819
frequency	0.052819	1.000000

Transaction frequencies and monetary values are uncorrelated; we can now fit our Gamma-Gamma model to predict average spend and expected lifetime values of our customers

The Gamma-Gamma model takes in a ‘data’ parameter, a pandas DataFrame with 3 columns representing Customer ID, average spend of repeat purchases, and number of repeat purchase for that customer. As with the BG/NBD model, these parameters are given HalfFlat priors which can be too diffuse for small datasets. For this example, we will use the default priors, but other priors can be specified just like with the BG/NBD example above.

gg = clv.GammaGammaModel()
gg.build_model(data=nonzero_data)
gg

Gamma-Gamma Model (Mean Transactions)
         p ~ Weibull(2, 1)
         q ~ Weibull(2, 1)
         v ~ Weibull(2, 10)
likelihood ~ Potential(f(q, p, v))

gg.fit(data=nonzero_data);

                                                                                                                   
                                                             Grad                                                  
  Progress               Draw        Divergen…   Step size   evals       Speed                Elapsed    Remaini…  
 ───────────────────────────────────────────────────────────────────────────────────────────────────────────────── 
  ━━━━━━━━━━━━━━━━━━━━   1400        0           0.375       3           1326.55 draws/s      0:00:01    0:00:00   
  ━━━━━━━━━━━━━━━━━━━━   1400        0           0.415       3           1255.07 draws/s      0:00:01    0:00:00   
  ━━━━━━━━━━━━━━━━━━━━   1400        0           0.439       7           1264.28 draws/s      0:00:01    0:00:00   
  ━━━━━━━━━━━━━━━━━━━━   1400        0           0.363       31          1314.31 draws/s      0:00:01    0:00:00   
                                                                                                                   

gg.fit_summary()

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
p	4.05	0.31	3.6	4.5	728	1002	1.00	0.011	0.0077
q	3.65	0.183	3.4	3.9	1134	1401	1.00	0.0053	0.0038
v	25.4	2.8	21	30	703	861	1.00	0.1	0.073

azp.plot_dist(gg.fit_result);

Predicting spend value of customers

Having fit our model, we can now use it to predict the conditional, expected average lifetime value of our customers, including those with zero repeat purchases.

expected_spend = gg.expected_customer_spend(data=rfm_data)

az.summary(expected_spend.isel(customer_id=range(10)), kind="stats")

	mean	sd	eti89_lb	eti89_ub
x[1]	26	0.26	26	27
x[2]	22	0.71	21	24
x[3]	39	0.94	37	40
x[4]	39	0.94	37	40
x[5]	39	0.94	37	40
x[6]	75	0.2	74	75
x[7]	22	0.71	21	24
x[8]	31	0.41	31	32
x[9]	37	0.16	36	37
x[10]	39	0.94	37	40

labeller = MapLabeller(var_name_map={"x": "customer"})
azp.plot_forest(
    xr.Dataset({"expected_spend": expected_spend.isel(customer_id=range(10))}),
    combined=True,
    labeller=labeller,
)
plt.xlabel("Expected mean spend");

We can also look at the average expected mean spend across all customers

az.summary(expected_spend.mean("customer_id"), kind="stats")

	mean	sd	eti89_lb	eti89_ub
x	39	0.59	38	40

azp.plot_dist(xr.Dataset({"expected_spend": expected_spend.mean("customer_id")}))
plt.axvline(expected_spend.mean(), color="k", ls="--")
plt.title("Expected mean spend of all customers");

Predicting spend value of a new customer

azp.plot_dist(xr.Dataset({"expected_spend": gg.expected_new_customer_spend()}))
plt.title("Expected mean spend of a new customer");

Estimating CLV

Finally, we can combine the GG with the BG/NBD model to obtain an estimate of the customer lifetime value. This relies on the discounted cash flow model, adjusting for cost of capital.

If computational issues are encountered, use the thin_fit_result method prior to estimating CLV.

bgm.thin_fit_result(keep_every=2)

BG/NBD
            alpha ~ Weibull(2, 10)
      phi_dropout ~ Uniform(0, 1)
    kappa_dropout ~ Pareto(1, 1)
                r ~ Weibull(2, 1)
                a = Deterministic(f(kappa_dropout, phi_dropout))
                b = Deterministic(f(kappa_dropout, phi_dropout))
recency_frequency ~ BetaGeoNBD(a, b, r, alpha, <constant>)

clv_estimate = gg.expected_customer_lifetime_value(
    transaction_model=bgm,
    data=rfm_data,
    future_t=12,  # months
    discount_rate=0.01,  # monthly discount rate ~ 12.7% annually
    time_unit="W",  # original data is in weeks
)

az.summary(clv_estimate.isel(customer_id=range(10)), kind="stats")

	mean	sd	eti89_lb	eti89_ub
x[1]	30	1.1	28	32
x[2]	3.1	0.3	2.6	3.6
x[3]	5.8	0.23	5.4	6.2
x[4]	5.8	0.23	5.4	6.2
x[5]	5.8	0.23	5.4	6.2
x[6]	510	15	490	530
x[7]	4.1	0.34	3.6	4.6
x[8]	16	0.43	16	17
x[9]	48	1.3	46	50
x[10]	5.8	0.23	5.4	6.2

azp.plot_forest(
    xr.Dataset({"clv": clv_estimate.isel(customer_id=range(10))}),
    combined=True,
    labeller=labeller,
)
plt.xlabel("Expected CLV");

According to our models, customer[6] has a much higher expected CLV. There is also a large variability in this estimate that arises solely from uncertainty in the parameters of the BG/NBD and GG models.

When there is considerable ambiguity and risk in marketing business decisions, the Bayesian modeling paradigm is very useful for elucidating the uncertainty of financial impacts and retention rates over time.

%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pymc,pytensor

Last updated: Wed, 05 Aug 2026

Python implementation: CPython
Python version       : 3.12.13
IPython version      : 9.15.0

pymc_marketing: 1.0.0.dev0
pymc          : 6.0.1
pytensor      : 3.0.7

arviz         : 1.2.0
arviz_base    : 1.2.0
arviz_plots   : 1.2.0
matplotlib    : 3.10.9
numpy         : 2.4.6
pandas        : 2.3.3
pymc_marketing: 1.0.0.dev0
xarray        : 2026.4.0

Watermark: 2.6.0