import arviz as az 
import matplotlib.pyplot as plt
import pandas as pd
import pymc as pm
import pymc_extras as pmx
import seaborn as sns

from pymc_marketing import clv

Cargar datos#

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

raw_data = pd.read_csv(url_cdnow)

Agregaciones RFM#

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

rfm_data.info()
rfm_data.head()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2357 entries, 0 to 2356
Data columns (total 4 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   customer_id  2357 non-null   int64  
 1   frequency    2357 non-null   float64
 2   recency      2357 non-null   float64
 3   T            2357 non-null   float64
dtypes: float64(3), int64(1)
memory usage: 73.8 KB

	customer_id	frequency	recency	T
0	1	3.0	49.0	78.0
1	2	1.0	2.0	78.0
2	3	0.0	0.0	78.0
3	4	0.0	0.0	78.0
4	5	0.0	0.0	78.0

ModeloParetoNBD#

ajuste del MAP#

pnbd_map = clv.ParetoNBDModel(data=rfm_data)
pnbd_map.build_model()  # required for prior predictive checks
pnbd_map

Pareto/NBD
            alpha ~ Weibull(2, 10)
             beta ~ Weibull(2, 10)
                r ~ Weibull(2, 1)
                s ~ Weibull(2, 1)
recency_frequency ~ ParetoNBD(r, alpha, s, beta, <constant>)

with pnbd_map.model:
    prior_idata = pm.sample_prior_predictive(random_seed=45, draws=1)

obs_freq = prior_idata.observed_data["recency_frequency"].sel(obs_var="frequency")
ppc_freq = prior_idata.prior_predictive["recency_frequency"].sel(obs_var="frequency")[
    0
][0]

Sampling: [alpha, beta, r, recency_frequency, s]

pnbd_map.fit()
map_fit = pnbd_map.fit_summary()  # save for plotting later

obs_freq = pnbd_map.idata.observed_data["recency_frequency"].sel(obs_var="frequency")
ppc_freq = pnbd_map.distribution_new_customer_recency_frequency(
    rfm_data,
    random_seed=42,
).sel(chain=0, draw=0, obs_var="frequency")

Sampling: [recency_frequency]

DEMZ se ajusta#

pnbd_full.fit_summary()

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	15.614	1.050	13.668	17.537	0.036	0.026	862.0	996.0	1.0
beta	12.750	3.691	6.362	19.946	0.135	0.096	739.0	1323.0	1.0
r	0.619	0.046	0.542	0.712	0.002	0.001	855.0	898.0	1.0
s	0.431	0.060	0.320	0.544	0.002	0.001	824.0	1251.0	1.0

Nutpie Fit#

pnbd_nutpie = clv.ParetoNBDModel(data=rfm_data)
bgnbd_nutpie = clv.BetaGeoModel(data=rfm_data)

idata_nutpie = pnbd_nutpie.fit(fit_method = 'mcmc', nuts_sampler="nutpie")

/Users/coltallen/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pymc/pytensorf.py:1066: FutureWarning: compile_pymc was renamed to compile. Old name will be removed in a future release of PyMC
  warnings.warn(

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[17], line 4
pnbd_nutpie = clv.ParetoNBDModel(data=rfm_data)
bgnbd_nutpie = clv.BetaGeoModel(data=rfm_data)
----> 4 idata_nutpie = pnbd_nutpie.fit(fit_method = 'mcmc', nuts_sampler="nutpie")

File ~/Projects/pymc-marketing/pymc_marketing/clv/models/pareto_nbd.py:350, in ParetoNBDModel.fit(self, fit_method, **kwargs)
with warnings.catch_warnings():
   warnings.simplefilter(
       action="ignore",
       category=UserWarning,
   )
--> 350     return super().fit(fit_method, **kwargs)

File ~/Projects/pymc-marketing/pymc_marketing/clv/models/basic.py:135, in CLVModel.fit(self, fit_method, **kwargs)
match fit_method:
   case "mcmc":
--> 135         idata = self._fit_mcmc(**kwargs)
   case "map":
       idata = self._fit_MAP(**kwargs)

File ~/Projects/pymc-marketing/pymc_marketing/clv/models/basic.py:164, in CLVModel._fit_mcmc(self, **kwargs)
   sampler_config = self.sampler_config.copy()
sampler_config.update(**kwargs)
--> 164 return pm.sample(**sampler_config, model=self.model)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pymc/sampling/mcmc.py:781, in sample(draws, tune, chains, cores, random_seed, progressbar, progressbar_theme, step, var_names, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, blas_cores, model, compile_kwargs, **kwargs)
       raise ValueError(
           "Model can not be sampled with NUTS alone. It either has discrete variables or a non-differentiable log-probability."
       )
   with joined_blas_limiter():
--> 781         return _sample_external_nuts(
           sampler=nuts_sampler,
           draws=draws,
           tune=tune,
           chains=chains,
           target_accept=kwargs.pop("nuts", {}).get("target_accept", 0.8),
           random_seed=random_seed,
           initvals=initvals,
           model=model,
           var_names=var_names,
           progressbar=progressbar,
           idata_kwargs=idata_kwargs,
           compute_convergence_checks=compute_convergence_checks,
           nuts_sampler_kwargs=nuts_sampler_kwargs,
           **kwargs,
       )
if exclusive_nuts and not provided_steps:
   # Special path for NUTS initialization
   if "nuts" in kwargs:

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pymc/sampling/mcmc.py:337, in _sample_external_nuts(sampler, draws, tune, chains, target_accept, random_seed, initvals, model, var_names, progressbar, idata_kwargs, compute_convergence_checks, nuts_sampler_kwargs, **kwargs)
   if kwarg in nuts_sampler_kwargs:
       compile_kwargs[kwarg] = nuts_sampler_kwargs.pop(kwarg)
--> 337 compiled_model = nutpie.compile_pymc_model(
   model,
   **compile_kwargs,
)
t_start = time.time()
idata = nutpie.sample(
   compiled_model,
   draws=draws,
   (...)
   **nuts_sampler_kwargs,
)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/nutpie/compile_pymc.py:391, in compile_pymc_model(model, backend, gradient_backend, **kwargs)
   backend = "numba"
if backend.lower() == "numba":
--> 391     return _compile_pymc_model_numba(model, **kwargs)
elif backend.lower() == "jax":
   return _compile_pymc_model_jax(
       model, gradient_backend=gradient_backend, **kwargs
   )

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/nutpie/compile_pymc.py:176, in _compile_pymc_model_numba(model, **kwargs)
   raise ImportError(
       "Numba is not installed in the current environment. "
       "Please install it with something like "
       "'mamba install -c conda-forge numba' "
       "and restart your kernel in case you are in an interactive session."
   )
import numba
(
   n_dim,
   n_expanded,
   logp_fn_pt,
   expand_fn_pt,
   shape_info,
--> 176 ) = _make_functions(model, mode="NUMBA", compute_grad=True, join_expanded=True)
expand_fn = expand_fn_pt.vm.jit_fn
logp_fn = logp_fn_pt.vm.jit_fn

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/nutpie/compile_pymc.py:497, in _make_functions(model, mode, compute_grad, join_expanded)
   (logp, grad) = pytensor.clone_replace([logp, grad], replacements)
   with model:
--> 497         logp_fn_pt = compile_pymc((joined,), (logp, grad), mode=mode)
else:
   (logp,) = pytensor.clone_replace([logp], replacements)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pymc/pytensorf.py:1070, in compile_pymc(*args, **kwargs)
def compile_pymc(*args, **kwargs):
   warnings.warn(
       "compile_pymc was renamed to compile. Old name will be removed in a future release of PyMC",
       FutureWarning,
   )
-> 1070     return compile(*args, **kwargs)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pymc/pytensorf.py:1055, in compile(inputs, outputs, random_seed, mode, **kwargs)
opt_qry = mode.provided_optimizer.including("random_make_inplace", check_parameter_opt)
mode = Mode(linker=mode.linker, optimizer=opt_qry)
-> 1055 pytensor_function = pytensor.function(
   inputs,
   outputs,
   updates={**rng_updates, **kwargs.pop("updates", {})},
   mode=mode,
   **kwargs,
)
return pytensor_function

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/compile/function/__init__.py:318, in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
   fn = orig_function(
       inputs, outputs, mode=mode, accept_inplace=accept_inplace, name=name
   )
else:
   # note: pfunc will also call orig_function -- orig_function is
   #      a choke point that all compilation must pass through
--> 318     fn = pfunc(
       params=inputs,
       outputs=outputs,
       mode=mode,
       updates=updates,
       givens=givens,
       no_default_updates=no_default_updates,
       accept_inplace=accept_inplace,
       name=name,
       rebuild_strict=rebuild_strict,
       allow_input_downcast=allow_input_downcast,
       on_unused_input=on_unused_input,
       profile=profile,
       output_keys=output_keys,
   )
return fn

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/compile/function/pfunc.py:465, in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys, fgraph)
   profile = ProfileStats(message=profile)
inputs, cloned_outputs = construct_pfunc_ins_and_outs(
   params,
   outputs,
   (...)
   fgraph=fgraph,
)
--> 465 return orig_function(
   inputs,
   cloned_outputs,
   mode,
   accept_inplace=accept_inplace,
   name=name,
   profile=profile,
   on_unused_input=on_unused_input,
   output_keys=output_keys,
   fgraph=fgraph,
)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/compile/function/types.py:1769, in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys, fgraph)
   m = Maker(
       inputs,
       outputs,
   (...)
       fgraph=fgraph,
   )
   with config.change_flags(compute_test_value="off"):
-> 1769         fn = m.create(defaults)
finally:
   if profile and fn:

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/compile/function/types.py:1661, in FunctionMaker.create(self, input_storage, storage_map)
start_import_time = pytensor.link.c.cmodule.import_time
with config.change_flags(traceback__limit=config.traceback__compile_limit):
-> 1661     _fn, _i, _o = self.linker.make_thunk(
       input_storage=input_storage_lists, storage_map=storage_map
   )
end_linker = time.perf_counter()
linker_time = end_linker - start_linker

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/basic.py:245, in LocalLinker.make_thunk(self, input_storage, output_storage, storage_map, **kwargs)
def make_thunk(
   self,
   input_storage: Optional["InputStorageType"] = None,
   (...)
   **kwargs,
) -> tuple["BasicThunkType", "InputStorageType", "OutputStorageType"]:
--> 245     return self.make_all(
       input_storage=input_storage,
       output_storage=output_storage,
       storage_map=storage_map,
   )[:3]

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/basic.py:695, in JITLinker.make_all(self, input_storage, output_storage, storage_map)
for k in storage_map:
   compute_map[k] = [k.owner is None]
--> 695 thunks, nodes, jit_fn = self.create_jitable_thunk(
   compute_map, nodes, input_storage, output_storage, storage_map
)
[fn] = thunks
fn.jit_fn = jit_fn

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/basic.py:647, in JITLinker.create_jitable_thunk(self, compute_map, order, input_storage, output_storage, storage_map)
# This is a bit hackish, but we only return one of the output nodes
output_nodes = [o.owner for o in self.fgraph.outputs if o.owner is not None][:1]
--> 647 converted_fgraph = self.fgraph_convert(
   self.fgraph,
   order=order,
   input_storage=input_storage,
   output_storage=output_storage,
   storage_map=storage_map,
)
thunk_inputs = self.create_thunk_inputs(storage_map)
thunk_outputs = [storage_map[n] for n in self.fgraph.outputs]

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/numba/linker.py:10, in NumbaLinker.fgraph_convert(self, fgraph, **kwargs)
def fgraph_convert(self, fgraph, **kwargs):
   from pytensor.link.numba.dispatch import numba_funcify
---> 10     return numba_funcify(fgraph, **kwargs)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
if not args:
   raise TypeError(f'{funcname} requires at least '
                   '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/numba/dispatch/basic.py:463, in numba_funcify_FunctionGraph(fgraph, node, fgraph_name, **kwargs)
@numba_funcify.register(FunctionGraph)
def numba_funcify_FunctionGraph(
   fgraph,
   (...)
   **kwargs,
):
--> 463     return fgraph_to_python(
       fgraph,
       numba_funcify,
       type_conversion_fn=numba_typify,
       fgraph_name=fgraph_name,
       **kwargs,
   )

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/utils.py:736, in fgraph_to_python(fgraph, op_conversion_fn, type_conversion_fn, order, storage_map, fgraph_name, global_env, local_env, get_name_for_object, squeeze_output, unique_name, **kwargs)
body_assigns = []
for node in order:
--> 736     compiled_func = op_conversion_fn(
       node.op, node=node, storage_map=storage_map, **kwargs
   )
   # Create a local alias with a unique name
   local_compiled_func_name = unique_name(compiled_func)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
if not args:
   raise TypeError(f'{funcname} requires at least '
                   '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/link/numba/dispatch/elemwise.py:357, in numba_funcify_Elemwise(op, node, **kwargs)
if not isinstance(op.scalar_op, Composite):
   scalar_inputs = [scalar(dtype=input.dtype) for input in node.inputs]
--> 357     scalar_node = op.scalar_op.make_node(*scalar_inputs)
scalar_op_fn = numba_funcify(
   op.scalar_op,
   node=scalar_node,
   (...)
   **kwargs,
)
nin = len(node.inputs)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/scalar/loop.py:180, in ScalarLoop.make_node(self, n_steps, *inputs)
cloned_constant = cloned_inputs[len(cloned_update) :]
# This will fail if the cloned init have a different dtype than the cloned_update
--> 180 op = ScalarLoop(
   init=cloned_init,
   update=cloned_update,
   constant=cloned_constant,
   until=cloned_until,
   name=self.name,
)
node = op.make_node(n_steps, *inputs)
return node

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/scalar/loop.py:69, in ScalarLoop.__init__(self, init, update, constant, until, name)
   inputs, outputs = clone([*init, *constant], update)
self.is_while = bool(until)
---> 69 self.inputs, self.outputs = self._cleanup_graph(inputs, outputs)
self._validate_updates(self.inputs, self.outputs)
self.inputs_type = tuple(input.type for input in self.inputs)

File ~/miniconda3/envs/pymc-marketing-dev/lib/python3.10/site-packages/pytensor/scalar/basic.py:3992, in ScalarInnerGraphOp._cleanup_graph(self, inputs, outputs)
for node in fgraph.apply_nodes:
   if not isinstance(node.op, ScalarOp):
-> 3992         raise TypeError(
           f"The fgraph of {self.__class__.__name__} must be exclusively "
           "composed of scalar operations."
       )
# Run MergeOptimization to avoid duplicated nodes
MergeOptimizer().rewrite(fgraph)

TypeError: The fgraph of ScalarLoop must be exclusively composed of scalar operations.

ADVI fit#

pnbd_advi.fit_summary()

arviz - WARNING - Shape validation failed: input_shape: (1, 500), minimum_shape: (chains=2, draws=4)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	15.547	0.567	14.514	16.599	0.030	0.021	366.0	453.0	NaN
beta	12.531	1.160	10.681	14.933	0.054	0.038	458.0	372.0	NaN
r	0.612	0.019	0.575	0.646	0.001	0.001	448.0	353.0	NaN
s	0.429	0.021	0.396	0.472	0.001	0.001	502.0	446.0	NaN

sns.lineplot(pnbd_advi.approx.hist, )
plt.yscale('log')

../../../_images/91c7b4cf415cae8c376314382f24087f0dfb480b5c9e4a2235c8cdc631e2d54e.png

Ajuste de rango completo#

pnbd_fullrank.fit_summary()

arviz - WARNING - Shape validation failed: input_shape: (1, 500), minimum_shape: (chains=2, draws=4)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	15.207	5.127	6.969	26.097	0.220	0.159	552.0	428.0	NaN
beta	14.941	10.196	2.099	33.381	0.503	0.394	467.0	363.0	NaN
r	0.592	0.209	0.253	0.973	0.010	0.007	391.0	409.0	NaN
s	0.417	0.192	0.143	0.786	0.009	0.007	458.0	462.0	NaN

sns.lineplot(pnbd_advi.approx.hist)

<Axes: >

../../../_images/c7e634de704f6c9fe294ce0a9f95fa8a01b0d7067c07f876e85207bcb3b9a910.png

Comparación visual#

_, axes = plt.subplots(
    nrows=2, ncols=2, figsize=(12, 7), sharex=False, sharey=False, layout="constrained"
)

axes = axes.flatten()

for i, var_name in enumerate(["r", "alpha", "s", "beta"]):
    ax = axes[i]
    az.plot_posterior(
        pnbd_full.idata.posterior[var_name].values.flatten(),
        color="C0",
        point_estimate="mean",
        ax=ax,
        label="DEMZ",
    )
    az.plot_posterior(
        pnbd_advi.idata.posterior[var_name].values.flatten(),
        color="C1",
        point_estimate="mean",
        ax=ax,
        label="ADVI",
    )
    # az.plot_posterior(
    #     pnbd_fullrank.idata.posterior[var_name].values.flatten(),
    #     color="C2",
    #     point_estimate="mean",
    #     ax=ax,
    #     label="FULLRANK_ADVI",
    # )
    ax.axvline(x=map_fit[var_name], color="C3", linestyle="--", label="MAP")
    ax.legend(loc="upper right")
    ax.set_title(var_name)

plt.gcf().suptitle("Pareto/NBD Model Parameters - DEMZ vs ADVI fits", fontsize=18, fontweight="bold");

../../../_images/eab1b0dfb159a79a01ec987bcfd926ccc70a9a257b4a9e44e6b7ba183949f475.png

_, axes = plt.subplots(
    nrows=2, ncols=2, figsize=(12, 7), sharex=False, sharey=False, layout="constrained"
)

axes = axes.flatten()

for i, var_name in enumerate(["r", "alpha", "s", "beta"]):
    ax = axes[i]
    az.plot_posterior(
        pnbd_full.idata.posterior[var_name].values.flatten(),
        color="C0",
        point_estimate="mean",
        ax=ax,
        label="DEMZ",
    )
    # az.plot_posterior(
    #     pnbd_advi.idata.posterior[var_name].values.flatten(),
    #     color="C1",
    #     point_estimate="mean",
    #     ax=ax,
    #     label="ADVI",
    # )
    az.plot_posterior(
        pnbd_fullrank.idata.posterior[var_name].values.flatten(),
        color="C2",
        point_estimate="mean",
        ax=ax,
        label="FULLRANK_ADVI",
    )
    ax.axvline(x=map_fit[var_name], color="C3", linestyle="--", label="MAP")
    ax.legend(loc="upper right")
    ax.set_title(var_name)

plt.gcf().suptitle("Pareto/NBD Model Parameters - DEMZ vs FULLRANK fits", fontsize=18, fontweight="bold");

../../../_images/eeb21844e208ebf25fa50ad1b73941e49d7db0508e4335a4623f5e4e825405e8.png

Observaciones:

Fullrank ofrece un ajuste bastante deficiente.
el ajuste de advi coincide con mcmc, aunque proporciona estimaciones más estrechas

Diferencias relativas en las estimaciones de parámetros#

(
    100*(pnbd_full.fit_summary()[['mean', 'sd']] - pnbd_advi.fit_summary()[['mean', 'sd']]) / 
    pnbd_full.fit_summary()[['mean', 'sd']]
    ).rename(
    columns={
        "mean": "pcnt_relative_diff_param_mean",
        "sd": "pcnt_relative_diff_param_sd",
    }
)

arviz - WARNING - Shape validation failed: input_shape: (1, 500), minimum_shape: (chains=2, draws=4)

	pcnt_relative_diff_param_mean	pcnt_relative_diff_param_sd
alpha	0.429102	46.000000
beta	1.717647	68.572203
r	1.130856	58.695652
s	0.464037	65.000000