Modelo gamma-gamma#

In this notebook we show how to fit a Gamma-Gamma model in PyMC-Marketing. The model is presented in the paper: Fader, P. S., & Hardie, B. G. (2013). The Gamma-Gamma model of monetary value. February, 2, 1-9.

Preparar el cuaderno#

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd

from pymc_marketing import clv

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

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

Cargar Datos#

Comenzamos cargando el conjunto de datos CDNOW.

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

summary_with_money_value = pd.read_csv(data_path)
summary_with_money_value["customer_id"] = summary_with_money_value.index
summary_with_money_value.head()

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

Para el modelo Gamma-Gamma, necesitamos filtrar a los clientes que han realizado solo una compra.

returning_customers_summary = summary_with_money_value.query("frequency > 0")

returning_customers_summary.head()

	frequency	recency	T	monetary_value	customer_id
0	2	30.43	38.86	22.35	0
1	1	1.71	38.86	11.77	1
5	7	29.43	38.86	73.74	5
6	1	5.00	38.86	11.77	6
8	2	35.71	38.86	25.55	8

Especificación del modelo#

Aquí describimos brevemente las suposiciones y la parametrización del modelo Gamma-Gamma del artículo mencionado.

El modelo de gasto por transacción se basa en las siguientes tres suposiciones generales:

El valor monetario de una transacción dada por un cliente varía aleatoriamente en torno a su valor promedio de transacción.
Los valores promedio de transacción varían entre los clientes, pero no varían con el tiempo para ningún individuo en particular.
La distribución de los valores promedio de transacción entre los clientes es independiente del proceso de transacción.

Para un cliente con x transacciones, sea \(z_1, z_2, \ldots, z_x\) el valor de cada transacción. El valor promedio de transacción observado del cliente por

\[::\]

Ahora describamos la parametrización:

Asumimos que \(z_i \sim \text{Gamma}(p, ν)\), con \(E(Z_i| p, ν) = \xi = p/ν\).

– Dadas las propiedades de convolución de la gamma, se deduce que el gasto total en x transacciones se distribuye como \(\text{Gamma}(px, ν)\).

– Dada la propiedad de escalado de la distribución gamma, se deduce que \(\bar{z} \sim \text{Gamma}(px, νx)\).
Asumimos \(ν \sim \text{Gamma}(q, \gamma)\).

Estamos interesados en estimar los parámetros \(p\), \(q\) y \(ν\).

Nota

El modelo Gamma-Gamma asume que no existe relación entre el valor monetario y la frecuencia de compra. Podemos verificar esta suposición calculando la correlación entre el gasto promedio y la frecuencia de compras.

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

	monetary_value	frequency
monetary_value	1.000000	0.113884
frequency	0.113884	1.000000

El valor de esta correlación es cercano a \(0.11\), lo cual en la práctica se considera lo suficientemente bajo como para proceder con el modelo.

Implementación de PyMC-Marketing#

Podemos utilizar la implementación preconstruida de PyMC-Marketing del modelo Gamma-Gamma, que también ofrece métodos de representación gráfica y predicción agradables:

Podemos construir el modelo para que podamos ver la especificación del modelo:

model = clv.GammaGammaModel()
model.build_model(data=returning_customers_summary)
model

Gamma-Gamma Model (Mean Transactions)
         p ~ HalfFlat()
         q ~ HalfFlat()
         v ~ HalfFlat()
likelihood ~ Potential(f(q, p, v))

Nota

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

Usando MAP#

Para comenzar, utilicemos un optimizador numérico (L-BFGS-B) de scipy.optimize para encontrar la estimación de máxima a posteriori (MAP) de los parámetros.

idata_map = model.fit(method="map").posterior.to_dataframe()

idata_map

		p	q	v
chain	draw
0	0	6.248787	3.744591	15.447813

MCMC#

También podemos utilizar MCMC para muestrear de la distribución posterior de los parámetros. MCMC es un método más robusto que MAP y proporciona estimaciones de incertidumbre para los parámetros.

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

idata_mcmc = model.fit(**sampler_kwargs)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [p, q, v]

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

idata_mcmc

arviz.InferenceData

Podemos ver algunas estadísticas de la distribución posterior de los parámetros.

model.fit_summary()

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	6.403	1.303	4.217	8.710	0.032	0.023	1686.0	1783.0	1.0
q	3.781	0.290	3.238	4.338	0.007	0.005	1753.0	2187.0	1.0
v	16.051	4.243	8.050	23.907	0.107	0.076	1563.0	1746.0	1.0

Visualicemos las distribuciones posteriores y el gráfico de rangos:

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

../../_images/ab636a9863d07ec6eac028ef7619214f79ea1fb3363903d64fcaeab72f805308.png

We can compare the MCMC posterior with the MAP estimation.

fig, axes = plt.subplots(
    nrows=3, ncols=1, figsize=(12, 10), sharex=False, sharey=False, layout="constrained"
)

for i, var_name in enumerate(["p", "q", "v"]):
    ax = axes[i]
    az.plot_posterior(
        idata_mcmc.posterior[var_name].values.flatten(),
        color="C0",
        point_estimate="mean",
        hdi_prob=0.95,
        ax=ax,
        label="MCMC",
    )
    ax.axvline(x=idata_map[var_name].item(), color="C2", linestyle="-.", label="MAP")
    ax.legend(loc="upper right")

plt.gcf().suptitle("Gamma-Gamma Model Parameters", fontsize=18, fontweight="bold");

../../_images/db1c36f5c4ceec3ddfd59ee56c07de2cad9dce23bb9a76ff10dffbdc5bc357c1.png

We see that the MAP estimates are close to the mean of the posterior distribution obtained by MCMC.

Gasto Esperado del Cliente#

Una vez que tengamos la distribución posterior de los parámetros, podemos utilizar el método expected_average_profit para calcular la expectativa condicional del beneficio promedio por transacción para un grupo de uno o más clientes.

expected_spend = model.expected_customer_spend(data=summary_with_money_value)

Veamos cómo se ve para un subconjunto de clientes.

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

	mean	sd	hdi_3%	hdi_97%
x[0]	24.706	0.512	23.839	25.762
x[1]	18.994	1.311	16.641	21.547
x[2]	35.195	0.924	33.447	36.887
x[3]	35.195	0.924	33.447	36.887
x[4]	35.195	0.924	33.447	36.887
x[5]	71.387	0.599	70.165	72.409
x[6]	18.994	1.311	16.641	21.547
x[7]	35.195	0.924	33.447	36.887
x[8]	27.318	0.394	26.627	28.107
x[9]	35.195	0.924	33.447	36.887

ax, *_ = az.plot_forest(
    data=expected_spend.isel(customer_id=(range(10))), combined=True, figsize=(8, 7)
)
ax.set(xlabel="Expected Spend (10 Customers)", ylabel="Customer ID")
ax.set_title("Expected Spend", fontsize=18, fontweight="bold");

../../_images/61af525a1b22cb28eb96c5b2b2d9df17ce79f5a1c55b79a2a44c677f41412ab5.png

Finalmente, veamos algunas estadísticas y la distribución para todo el conjunto de datos.

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

	mean	sd	hdi_3%	hdi_97%
x	35.268	0.629	34.107	36.442

fig, ax = plt.subplots()
az.plot_dist(expected_spend.mean("customer_id"), label="Mean over Customer ID", ax=ax)
ax.axvline(x=expected_spend.mean(), color="black", ls="--", label="Overall Mean")
ax.legend(loc="upper right")
ax.set(xlabel="Expected Spend", ylabel="Density")
ax.set_title("Expected Spend", fontsize=18, fontweight="bold");

../../_images/b5a6e84270a0845c1558d4f6609653c9140b901681c95fadf97f1a897d25084d.png

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

Last updated: Mon Jul 01 2024

Python implementation: CPython
Python version       : 3.10.13
IPython version      : 8.22.2

pymc    : 5.14.0
pytensor: 2.20.0

arviz         : 0.17.1
pymc_marketing: 0.6.0
matplotlib    : 3.8.3
pandas        : 2.2.1

Watermark: 2.4.3