Modelo geométrico beta desplazado con datos individuales de clientes#
En este cuaderno replicamos los principales resultados y figuras de
Fader, P. S., & Hardie, B. G. (2007). Cómo proyectar la retención de clientes. Journal of Interactive Marketing, 21(1), 76-90. https://journals.sagepub.com/doi/pdf/10.1002/dir.20074
Los autores describen el modelo Beta Geométrico Desplazado (sBG) para el comportamiento del cliente en un entorno contractual discreto. Se asume que:
Al final de cada período, un cliente tiene una probabilidad
thetade renovar el contrato y1-thetade cancelar.La probabilidad
thetano cambia con el tiempo para un cliente dadoLa probabilidad
thetavaría entre los clientes de acuerdo con una distribución a priori Beta con hiperparámetrosalphaybeta.
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
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"
seed = sum(map(ord, "Individual sBG Model"))
rng = np.random.default_rng(seed)
Recreando el conjunto de datos#
El conjunto de datos contiene el porcentaje de clientes que aún están inscritos en el servicio después de cada período de tiempo. Hay dos grupos distintos de usuarios: regular y de alta gama, que se cree que tienen diferentes tasas de deserción. Cada grupo tenía 1000 clientes al inicio.
Al igual que en el documento original, ajustaremos los datos solo de los primeros 7 períodos de tiempo y utilizaremos los siguientes períodos para la validación.
df = pd.DataFrame(
{
"regular": [
100.0,
63.1,
46.8,
38.2,
32.6,
28.9,
26.2,
24.1,
22.3,
20.7,
19.4,
18.3,
17.3,
],
"highend": [
100.0,
86.9,
74.3,
65.3,
59.3,
55.1,
51.7,
49.1,
46.8,
44.5,
42.7,
40.9,
39.4,
],
}
)
df
| regular | highend | |
|---|---|---|
| 0 | 100.0 | 100.0 |
| 1 | 63.1 | 86.9 |
| 2 | 46.8 | 74.3 |
| 3 | 38.2 | 65.3 |
| 4 | 32.6 | 59.3 |
| 5 | 28.9 | 55.1 |
| 6 | 26.2 | 51.7 |
| 7 | 24.1 | 49.1 |
| 8 | 22.3 | 46.8 |
| 9 | 20.7 | 44.5 |
| 10 | 19.4 | 42.7 |
| 11 | 18.3 | 40.9 |
| 12 | 17.3 | 39.4 |
El modelo pymc-marketing implementado ajusta cada theta individual, en contraste con el enfoque del artículo donde los parámetros individuales son marginalizados.
Para nuestro análisis, tendremos que crear un conjunto de datos sintético con datos individuales, siguiendo las tendencias del grupo agregado.
def individual_data_from_percentage_alive(percentage_alive, initial_customers):
n_alive = np.asarray(percentage_alive / 100 * initial_customers, dtype=int)
died_at = np.zeros((initial_customers,), dtype=int)
counter = 0
for t, diff in enumerate((n_alive[:-1] - n_alive[1:]), start=1):
died_at[counter : counter + diff] = t
counter += diff
censoring_t = t + 1
died_at[counter:] = censoring_t
return died_at
T = 8 # In the paper only the first 7 years are included in the model
truncated_df = df[:T]
initial_customers = 1000 # Seems to be what F&H use, according to Appendix B
churn_highend = individual_data_from_percentage_alive(
truncated_df["highend"], initial_customers
)
churn_regular = individual_data_from_percentage_alive(
truncated_df["regular"], initial_customers
)
churn_regular[::10]
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
4, 4, 5, 5, 5, 5, 6, 6, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8])
Los clientes que no han abandonado tienen t_churn == T == 8. El modelo manejará estos de manera diferente, ya que sus vidas son censuradas.
Ajustando los modelos sBG para usuarios de alta gama y usuarios regulares#
dataset_highend = pd.DataFrame(
{
"customer_id": np.arange(initial_customers) + 1,
"t_churn": churn_highend,
"T": T,
"cohorts": "highend",
}
)
highend_users = clv.ShiftedBetaGeoModelIndividual()
highend_users.build_model(data=dataset_highend)
highend_users
Shifted-Beta-Geometric Model (Individual Customers)
alpha ~ HalfFlat()
beta ~ HalfFlat()
theta ~ Beta(alpha, beta)
churn_censored ~ Censored(Geometric(theta), -inf, <constant>)
dataset_regular = pd.DataFrame(
{
"customer_id": np.arange(initial_customers) + 1001,
"t_churn": churn_regular,
"T": T,
"cohorts": "regular",
}
)
regular_users = clv.ShiftedBetaGeoModelIndividual()
regular_users.build_model(data=dataset_regular)
regular_users
Shifted-Beta-Geometric Model (Individual Customers)
alpha ~ HalfFlat()
beta ~ HalfFlat()
theta ~ Beta(alpha, beta)
churn_censored ~ Censored(Geometric(theta), -inf, <constant>)
highend_users.fit(random_seed=rng, nuts_sampler="nutpie")
regular_users.fit(random_seed=rng, nuts_sampler="nutpie");
Sampler Progress
Total Chains: 4
Active Chains: 0
Finished Chains: 4
Sampling for now
Estimated Time to Completion: now
| Progress | Draws | Divergences | Step Size | Gradients/Draw |
|---|---|---|---|---|
| 2000 | 0 | 0.23 | 15 | |
| 2000 | 0 | 0.24 | 15 | |
| 2000 | 0 | 0.23 | 15 | |
| 2000 | 0 | 0.23 | 15 |
Sampler Progress
Total Chains: 4
Active Chains: 0
Finished Chains: 4
Sampling for now
Estimated Time to Completion: now
| Progress | Draws | Divergences | Step Size | Gradients/Draw |
|---|---|---|---|---|
| 2000 | 0 | 0.26 | 15 | |
| 2000 | 0 | 0.25 | 15 | |
| 2000 | 0 | 0.26 | 15 | |
| 2000 | 0 | 0.26 | 15 |
Contrastando inferencias posteriores con las estimaciones MLE del repositorio#
El modelo sBG tiene 2 parámetros poblacionales de interés: alpha y beta. Estos parámetros definen la distribución poblacional de las tasas de abandono theta de cada cliente individual. Cuanto mayores sean los valores de alpha y beta, más homogéneas serán las tasas de abandono entre diferentes clientes.
La relación de alpha a beta nos indica las tasas de cancelación esperadas. Si alpha/beta == 0.1, esperamos que el cliente promedio tenga una probabilidad de 0.1 de cancelar entre cada período de tiempo.
# MLE estimates from the paper
ref_val = {
"highend": [0.704, 3.806],
"regular": [0.688, 1.182],
}
highend_users.fit_summary(var_names=["alpha", "beta"])
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 0.742 | 0.139 | 0.495 | 0.985 | 0.017 | 0.011 | 71.0 | 152.0 | 1.06 |
| beta | 4.400 | 1.087 | 2.530 | 6.304 | 0.128 | 0.087 | 73.0 | 141.0 | 1.06 |
az.plot_posterior(
highend_users.fit_result, var_names=["alpha", "beta"], ref_val=ref_val["highend"]
);
regular_users.fit_summary(var_names=["alpha", "beta"])
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 0.723 | 0.067 | 0.608 | 0.854 | 0.004 | 0.002 | 225.0 | 589.0 | 1.02 |
| beta | 1.228 | 0.157 | 0.945 | 1.507 | 0.011 | 0.007 | 214.0 | 577.0 | 1.02 |
az.plot_posterior(
regular_users.fit_result, var_names=["alpha", "beta"], ref_val=ref_val["regular"]
);
El ajuste del modelo coincide con las estimaciones de Máxima Verosimilitud descritas en el documento original. Además, el muestreo MCMC nos proporciona información útil sobre la incertidumbre de los ajustes.
Recreando los resultados y figuras principales del artículo#
La figura 6 del artículo muestra la distribución grupal de las tasas de abandono de theta.
Las distribuciones inferidas son relativamente amplias en ambos grupos, lo que se deriva directamente del hecho de que tanto alpha como beta se inferieron como < 10.
Sigue siendo claro que la distribución de theta para los usuarios de alto nivel tiene mucha más masa cerca de valores más bajos, lo que sugiere que los usuarios de alto nivel tienen tasas de deserción más bajas en promedio, en comparación con los usuarios regulares.
new_highend_theta = highend_users.distribution_new_customer_theta(
n=100, random_seed=rng
)
new_regular_theta = regular_users.distribution_new_customer_theta(
n=100, random_seed=rng
);
Sampling: [theta]
Sampling: [theta]
ax = az.plot_dist(new_highend_theta, label="highend")
ax = az.plot_dist(new_regular_theta, label="regular", ax=ax, color="C1")
ax.set_xlabel("theta")
plt.ylim([0, 4])
ax.set_title("Figure 6");
Las figuras 4 y 5 muestran las tendencias promedio de abandono y retención predichas para los dos grupos. Podemos ver que las predicciones coinciden bien con los datos observados (línea negra), incluso al extrapolar en los períodos de tiempo que se excluyeron al ajustar el modelo.
Las gráficas también destacan una implicación interesante del modelo: se espera que las tasas de retención aumenten con el tiempo, a medida que los clientes más precarios se vayan retirando gradualmente. Esta es una consecuencia directa de modelar las tasas de abandono individuales como fijas a lo largo del tiempo.
churn_highend = highend_users.distribution_customer_churn_time(
customer_id=np.arange(initial_customers),
random_seed=rng,
)
churn_regular = regular_users.distribution_customer_churn_time(
customer_id=np.arange(initial_customers),
random_seed=rng,
)
Sampling: [churn, theta]
Sampling: [churn, theta]
weeks = xr.DataArray(np.arange(1, 14), dims=["weeks"])
survival_perc_highend = 100 - 100 * (churn_highend < weeks)
survival_perc_regular = 100 - 100 * (churn_regular < weeks)
az.plot_hdi(
weeks,
survival_perc_highend.mean("customer_id"),
hdi_prob=0.95,
color="C0",
fill_kwargs={"label": "high end"},
)
az.plot_hdi(
weeks,
survival_perc_regular.mean("customer_id"),
hdi_prob=0.95,
color="C1",
fill_kwargs={"label": "regular"},
)
plt.plot(weeks, df["highend"], color="k", label="observed")
plt.plot(weeks, df["regular"], color="k")
plt.axvline(7, ls=":", color="k")
plt.ylim([0, 105])
plt.ylabel("% Surviving")
plt.legend()
plt.title("Figure 4");
retention_rate_highend = survival_perc_highend.isel(
weeks=slice(1, None, None)
) / survival_perc_highend.isel(weeks=slice(None, -1, None))
retention_rate_regular = survival_perc_regular.isel(
weeks=slice(1, None, None)
) / survival_perc_regular.isel(weeks=slice(None, -1, None))
retention_rate_highend_obs = df["highend"][1:].values / df["highend"][:-1].values
retention_rate_regular_obs = df["regular"][1:].values / df["regular"][:-1].values
weeks_ = weeks[:-1]
az.plot_hdi(
weeks_,
retention_rate_highend.mean("customer_id"),
hdi_prob=0.95,
color="C0",
fill_kwargs={"label": "high end"},
)
az.plot_hdi(
weeks_,
retention_rate_regular.mean("customer_id"),
hdi_prob=0.95,
color="C1",
fill_kwargs={"label": "regular"},
)
plt.plot(weeks_, retention_rate_highend_obs, color="k", label="observed")
plt.plot(weeks_, retention_rate_regular_obs, color="k")
plt.axvline(7, ls="--", color="k")
plt.ylim([0.5, 1.05])
plt.ylabel("Retention Rate")
plt.legend()
plt.title("Figure 5");
%load_ext watermark
%watermark -n -u -v -iv -w -p pymc,pytensor
The watermark extension is already loaded. To reload it, use:
%reload_ext watermark
Last updated: Thu Oct 16 2025
Python implementation: CPython
Python version : 3.12.11
IPython version : 9.4.0
pymc : 5.25.1
pytensor: 2.31.7
arviz : 0.22.0
numpy : 2.2.6
pandas : 2.3.1
xarray : 2025.7.1
pymc_marketing: 0.15.1
pytensor : 2.31.7
matplotlib : 3.10.3
Watermark: 2.5.0