Implementación del modelo Gamma-gamma en PyMC por el desarrollador#

Referencia: Fader, P. S., & Hardie, B. G. (2013). El modelo Gamma-Gamma del valor monetario. Febrero, 2, 1-9.

http://www.brucehardie.com/notes/025/gamma_gamma.pdf

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
import seaborn as sns

from pymc_marketing import clv

Simular datos#

rng = np.random.default_rng(42)

# Hyperparameters
p_true = 6.
q_true = 4.
v_true = 15.

# Number of subjects
N = 500  
# Subject level parameters
nu_true = pm.draw(pm.Gamma.dist(q_true, v_true, size=N), random_seed=rng)

# Number of observations per subject
x = rng.poisson(lam=2, size=N) + 1  
idx = np.repeat(np.arange(0, N), x)
# Observations
z = pm.draw(pm.Gamma.dist(p_true, nu_true[idx]), random_seed=rng)

print(sum(x))
assert len(nu_true[idx]) == sum(x)

plt.hist(z, bins=50, ec="w")
plt.xlabel("transaction value")
plt.ylabel("counts")
plt.title("Simulated data");

../../../_images/3061ffa4807d34bcb2812afa4f34cc89a0efedc9a37ae37c2d34b5af09e60e13.png

df = pd.DataFrame(data={"individual_transaction_value": z, "customer_id": idx})
z_mean = df.groupby("customer_id").mean()["individual_transaction_value"].values
z_mean[:10]

array([ 17.5597973 ,  41.05272046,  15.90609488,  83.95307047,
        20.36896009,  23.8572992 ,  46.09000842,  47.49876237,
       131.16095313,  16.42659393])

Implementación de PyMC#

Podemos utilizar la implementación preconstruida de PyMMMC del modelo Gamma-Gamma, que también ofrece métodos de graficado y predicción atractivos.

Implementaciones manuales de PyMC#

Mostramos cómo se puede implementar el modelo Gamma-Gamma manualmente utilizando PyMC. Esto aclara cómo se puede modificar o extender el modelo para incluir más información previa o estructura adicional.

Modelo Gamma-Gamma condicionado a transacciones individuales \(z\)#

with pm.Model() as m1:
    p = pm.HalfFlat("p")
    q = pm.HalfFlat("q")
    v = pm.HalfFlat("v")
    
    nu = pm.Gamma("nu", q, v, size=N)
    pm.Gamma("z", p, nu[idx], observed=z)

    pm.Deterministic("mean_spend", p / nu)
    
    trace1 = pm.sample(random_seed=rng)

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

100.00% [8000/8000 00:20<00:00 Sampling 4 chains, 0 divergences]

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

az.summary(trace1, var_names=["p", "q", "v"])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	6.055	0.260	5.578	6.564	0.020	0.014	178.0	382.0	1.02
q	3.914	0.295	3.358	4.453	0.007	0.005	1877.0	1959.0	1.00
v	14.879	1.486	12.194	17.717	0.071	0.050	438.0	1131.0	1.01

az.plot_posterior(trace1, var_names=["p", "q", "v"], ref_val=[p_true, q_true, v_true]);

../../../_images/27701817275b87fbd5dee64d0c0fb99db7e4a71e809946acc972155a158e56ec.png

Modelo gamma-gamma condicionado a las transacciones promedio por usuario \(\overline{z}\)#

Esto no puede muestrear porque el modelo contiene «casi» dos parámetros independientes por observación. Para más detalles, consulta este tema de Discourse

with pm.Model() as m2:
    p = pm.HalfFlat("p")
    q = pm.HalfFlat("q")
    v = pm.HalfFlat("v")

    nu = pm.Gamma("nu", q, v, size=N)
    # We use the convolution properties of the gamma distribution to model
    # the mean of multiple transaction using the parameters of individual
    # transactions
    pm.Gamma("z_mean", p*x, nu*x, observed=z_mean)
    
    trace2 = pm.sample(random_seed=rng)

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

100.00% [8000/8000 00:26<00:00 Sampling 4 chains, 0 divergences]

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

az.summary(trace2, var_names=["p", "q", "v"])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	14.454	11.819	4.299	40.762	5.738	4.375	5.0	11.0	2.17
q	3.629	0.376	2.939	4.294	0.129	0.094	8.0	41.0	1.41
v	9.154	5.062	1.340	16.853	2.264	1.705	5.0	11.0	2.09

az.plot_posterior(trace2, var_names=["p", "q", "v"], ref_val=[p_true, q_true, v_true]);

../../../_images/6c9e155c1ba21aaf01e564ccabcda2b201da92591a80fc016dd1323bb7652817.png

Modelo Gamma-Gamma condicionado en la transacción promedio por usuario con \(\nu\) marginalizado#

with pm.Model() as m3:
    p = pm.HalfFlat("p")
    q = pm.HalfFlat("q")
    v = pm.HalfFlat("v")

    # Likelihood of z_mean, marginalizing over nu
    likelihood = pm.Potential(
        "likelihood", 
        (
            pt.gammaln(p * x + q)
            - pt.gammaln(p * x)
            - pt.gammaln(q)
            + q * pt.log(v)
            + (p * x - 1) * pt.log(z_mean)
            + (p * x) * pt.log(x)
            - (p * x + q) * pt.log(x * z_mean + v)
        ),
    )

    # Closed form solution posterior individual nu
    nu = pm.Deterministic("nu", pm.Gamma.dist(p * x + q, v + x * z_mean))
    pm.Deterministic("mean_spend", p / nu)
    
    trace3 = pm.sample(random_seed=rng)

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

100.00% [8000/8000 00:32<00:00 Sampling 4 chains, 0 divergences]

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

az.summary(trace3, var_names=["p", "q", "v"])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	6.845	2.937	3.153	12.110	0.110	0.080	865.0	812.0	1.0
q	3.915	0.390	3.194	4.654	0.012	0.009	974.0	1118.0	1.0
v	15.375	6.510	4.445	26.879	0.221	0.156	816.0	738.0	1.0

az.plot_posterior(trace3, var_names=["p", "q", "v"], ref_val=[p_true, q_true, v_true]);

../../../_images/6d766e3b0f5cf530430be6303cb8bf28db990e8966f0d1c358bb1c925f311962.png