Modelo BG/NBD#

In this notebook we show how to fit a BG/NBD model in PyMC-Marketing. The model is presented in the paper: Fader, P. S., Hardie, B. G., & Lee, K. L. (2005). “Counting your customers” the easy way: An alternative to the Pareto/NBD model. Marketing science, 24(2), 275-284.

Preparar Notebook#

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
from fastprogress.fastprogress import progress_bar

from pymc_marketing import clv

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

%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "retina"

Leer datos#

We use the CDNOW dataset.

data_path = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/clv_quickstart.csv"

df = pd.read_csv(data_path)

df.head()

	frequency	recency	T	monetary_value
0	2	30.43	38.86	22.35
1	1	1.71	38.86	11.77
2	0	0.00	38.86	0.00
3	0	0.00	38.86	0.00
4	0	0.00	38.86	0.00

The following definitions apply:

frequency representa el número de compras repetidas que ha realizado el cliente. Esto significa que es uno menos que el número total de compras. En realidad, esto es ligeramente incorrecto. Es el conteo de los períodos de tiempo en los que el cliente realizó una compra. Así que, si se utilizan días como unidades, entonces es el conteo de días en los que el cliente realizó una compra.

T representa la edad del cliente en las unidades de tiempo elegidas (semanalmente, en el conjunto de datos anterior). Esto es igual a la duración entre la primera compra de un cliente y el final del período bajo estudio.

recency representa la antigüedad del cliente en el momento de sus compras más recientes. Esto es igual a la duración entre la primera compra de un cliente y su última compra. (Por lo tanto, si solo ha realizado 1 compra, la recency es 0.)

Truco

Renombramos la columna del índice a customer_id, ya que esto es requerido por el modelo.

data = (
    df.reset_index()
    .rename(columns={"index": "customer_id"})
    .drop(columns="monetary_value")
)

Especificación del Modelo#

El modelo BG/NBD es un modelo probabilístico que describe el comportamiento de compra de un cliente en un entorno no contractual. Se basa en las siguientes suposiciones para cada cliente:

Proceso de Frecuencia#

Mientras está activo, el tiempo entre transacciones se distribuye de manera exponencial con la tasa de transacciones, es decir,

\[f(t_{j}|t_{j-1}; \lambda) = \lambda \exp(-\lambda (t_{j} - t_{j - 1})), \quad t_{j} \geq t_{j - 1} \geq 0\]
La heterogeneidad en \(\lambda\) sigue una distribución gamma con pdf

\[f(\lambda|r, \alpha) = \frac{\alpha^{r}\lambda^{r - 1}\exp(-\lambda \alpha)}{\Gamma(r)}, \quad \lambda > 0\]

Proceso de abandono#

Después de cualquier transacción, un cliente se vuelve inactivo con una probabilidad de \(p\).
La heterogeneidad en \(p\) sigue una distribución beta con pdf

\[f(p|a, b) = \frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} p^{a - 1}(1 - p)^{b - 1}, \quad 0 \leq p \leq 1\]
La tasa de transacción \(\lambda\) y la probabilidad de abandono \(p\) varían de manera independiente entre los clientes.

En lugar de estimar \(\lambda\) y \(p\) para cada cliente específico, lo hacemos para un cliente elegido al azar, es decir, trabajamos con los valores esperados de los parámetros. Por lo tanto, estamos interesados en encontrar la distribución posterior de los parámetros \(r\), \(\alpha\), \(a\) y \(b\).

Ajuste del Modelo#

Estimar tales parámetros es muy fácil en PyMC-Marketing. Instanciamos el modelo de manera similar:

model = clv.BetaGeoModel()

Y construya el modelo para ver la configuración del modelo:

model.build_model(data=data)
model

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

Observe los priors adicionales phi_dropout y kappa_dropout. Estos se añadieron a la configuración predeterminada para mejorar el rendimiento, pero se pueden omitir al especificar una model_config personalizada con a y b.

La estructura del modelo especificado también se puede visualizar:

model.graphviz()

../../_images/69248deb786fb470404348e4991c4d9f47d0d080b33ce90bfe16d649bf0a89c9.svg

Ahora podemos ajustar el modelo. El muestreador predeterminado en PyMC-Marketing es el Muestreador No-U-Turn (NUTS). Utilizamos las \(4\) cadenas predeterminadas y \(1000\) muestras por cadena.

Nota

No es necesario construir el modelo antes de ajustarlo. Podemos ajustar el modelo directamente.

sample_kwargs = {
    "draws": 2_000,
    "chains": 4,
    "target_accept": 0.9,
    "random_seed": 42,
}

idata_mcmc = model.fit(**sample_kwargs)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, r, phi_dropout, kappa_dropout]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 7 seconds.

Podemos consultar la tabla de resumen:

model.fit_summary()

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	4.504	0.374	3.810	5.201	0.006	0.004	4154.0	5019.0	1.0
r	0.245	0.012	0.222	0.268	0.000	0.000	4217.0	4680.0	1.0
phi_dropout	0.248	0.020	0.212	0.287	0.000	0.000	5004.0	4984.0	1.0
kappa_dropout	3.198	0.932	1.734	4.918	0.013	0.010	5262.0	4991.0	1.0
a	0.782	0.192	0.470	1.148	0.002	0.002	6141.0	5313.0	1.0
b	2.416	0.750	1.237	3.775	0.011	0.008	5096.0	4821.0	1.0

Vemos que los valores de r_hat están cerca de \(1\), lo que indica convergencia.

También podemos trazar distribuciones posteriores de los parámetros y los gráficos de rango:

axes = az.plot_trace(
    data=model.idata,
    compact=True,
    kind="rank_bars",
    backend_kwargs={"figsize": (12, 9), "layout": "constrained"},
)
plt.gcf().suptitle("BG/NBD Model Trace", fontsize=18, fontweight="bold");

../../_images/10f1c2afbd27dd655e7965556369cac63085f776e48cd6a1be0b48be79d7c500.png

Comprobaciones Predictivas Anteriores y Posteriores#

Los PPCs nos permiten verificar la eficacia de nuestras distribuciones a priori y el rendimiento de los posteriors ajustados.

Veamos cómo se desempeña el modelo en un chequeo predictivo previo, donde muestreamos de las distribuciones a priori predeterminadas antes de ajustar el modelo:

# PPC histogram plot
clv.plot_expected_purchases_ppc(model, ppc="prior");

Sampling: [alpha, kappa_dropout, phi_dropout, r, recency_frequency]

../../_images/7b4ba9f38c1ad25d61a45c9deebad1e05f5cf2d39b9a9be7461b989dd6efb1ee.png

clv.plot_expected_purchases_ppc(model, ppc="posterior");

Sampling: [recency_frequency]

../../_images/f4fdcb315b8aa832f8b3300d4e42b0157a13e5032f27ff4b41483b0a2a38895b.png

Algunas Aplicaciones#

Ahora que ha ajustado el modelo, podemos utilizarlo para hacer predicciones. Por ejemplo, podemos predecir la probabilidad esperada de que un cliente esté vivo en función del tiempo (pasos). Aquí hay un fragmento de código para hacerlo:

Número Esperado de Compras#

Tomemos una muestra de usuarios:

example_customer_ids = [1, 6, 10, 18, 45, 1412]

data_small = data.query("customer_id.isin(@example_customer_ids)")

data_small.head(6)

	customer_id	frequency	recency	T
1	1	1	1.71	38.86
6	6	1	5.00	38.86
10	10	5	24.43	38.86
18	18	3	28.29	38.71
45	45	12	34.43	38.57
1412	1412	14	30.29	31.57

Observe que los últimos dos clientes son compradores frecuentes en comparación con los demás.

steps = 90

expected_num_purchases_steps = xr.concat(
    objs=[
        model.expected_purchases(
            data=data_small,
            future_t=t,
        )
        for t in progress_bar(range(steps))
    ],
    dim="t",
).transpose(..., "t")

100.00% [90/90 00:01<00:00]

Podemos trazar el número esperado de compras para los próximos \(90\) períodos:

../../_images/77dcc58b90ac4f7a711b39f00d0d5c28c8d49e0277e84a78accc4a4531c46c8c.png

Tenga en cuenta que se espera que los compradores frecuentes realicen más compras en el futuro.

Probabilidad de que un Cliente esté Vivo#

Ahora examinamos la probabilidad de que un cliente esté vivo durante los próximos \(90\) períodos:

steps = 90

future_alive_all = []

for t in progress_bar(range(steps)):
    future_data = data_small.copy()
    future_data["T"] = future_data["T"] + t
    future_alive = model.expected_probability_alive(data=future_data)
    future_alive_all.append(future_alive)

expected_probability_alive_steps = xr.concat(
    objs=future_alive_all,
    dim="t",
).transpose(..., "t")

100.00% [90/90 00:00<00:00]

../../_images/d1439ffc49816e90a41d2915efcab05ea9731fe6f5fd050108eef5a82629eb66.png

Truco

Aquí hay algunas observaciones generales:

Estos gráficos suponen que no habrá compras futuras.
La probabilidad de descomposición no es la misma, ya que depende del historial de compras del cliente.
La probabilidad de estar vivo siempre está disminuyendo, ya que asumimos que no hay cambios en los otros parámetros.
Estas probabilidades son siempre no negativas, como se esperaba.

Advertencia

Para los compradores frecuentes, la probabilidad de estar vivos disminuye muy rápidamente ya que asumimos que no habrá compras futuras. Es muy importante tener esto en cuenta al interpretar los resultados.

Probabilidad de que un cliente no realice ninguna compra en un rango de tiempo#

Ahora analizamos la probabilidad de que un cliente realice 0 compras entre ahora y los próximos \(t\) períodos, que van de 0 a 30.

steps = 30
expected_probability_zero_purchases = xr.concat(
    objs=[
        model.expected_probability_no_purchase(
            data=data_small,
            t=t,
        )
        for t in progress_bar(range(steps))
    ],
    dim="t",
).transpose(..., "t")

100.00% [30/30 00:00<00:00]

fig, axes = plt.subplots(
    nrows=len(example_customer_ids),
    ncols=1,
    figsize=(12, 15),
    sharex=True,
    sharey=True,
    layout="constrained",
)

axes = axes.flatten()

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i]
    customer_expected_probability_zero_purchases = (
        expected_probability_zero_purchases.sel(customer_id=customer_id)
    )
    az.plot_hdi(
        range(steps),
        customer_expected_probability_zero_purchases,
        hdi_prob=0.94,
        color="C1",
        fill_kwargs={"alpha": 0.3, "label": "$94 \\%$ HDI"},
        ax=ax,
    )
    az.plot_hdi(
        range(steps),
        customer_expected_probability_zero_purchases,
        hdi_prob=0.5,
        color="C1",
        fill_kwargs={"alpha": 0.5, "label": "$50 \\%$ HDI"},
        ax=ax,
    )
    ax.plot(
        range(steps),
        customer_expected_probability_zero_purchases.mean(dim=("chain", "draw")),
        color="C1",
        label="posterior mean",
    )
    ax.legend(loc="upper right")
    ax.set(title=f"Customer {customer_id}", ylabel="Probability", ylim=(0, 1))

axes[-1].set(xlabel="steps")
plt.gcf().suptitle(
    "Expected Probability Zero Purchases between $(T, T+t]$.",
    fontsize=18,
    fontweight="bold",
);

../../_images/1681a3544a4c9e97fdad538ed97a04c403980f565bfe1c59e0a96c5b86175a9c.png

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

Last updated: Wed Feb 05 2025

Python implementation: CPython
Python version       : 3.10.16
IPython version      : 8.31.0

pymc    : 5.19.1
pytensor: 2.26.4

arviz         : 0.20.0
pandas        : 2.2.3
pymc_marketing: 0.11.0
xarray        : 2025.1.1
matplotlib    : 3.10.0
fastprogress  : 1.0.3

Watermark: 2.5.0