<no title>

import numpy as np
import pandas as pd
import pymc as pm

from pymc_marketing.clv.distributions import continuous_contractual, continuous_non_contractual

from scipy.special import expit, hyp2f1

import xarray as xr

import matplotlib.pyplot as plt

a = 0.8
b = 2.5
alpha = 3
r = 4

rng = np.random.default_rng(seed=34)

df = pd.read_csv("data/clv_quickstart.csv")
T = df["T"].values

def rng_fn(rng, a, b, r, alpha, T, T0, size):
    p = rng.beta(a, b, size=size)
    lam = rng.gamma(r, 1 / alpha, size=size)

    return continuous_contractual.rng_fn(rng, lam, p, T, T0, size=size)

data = rng_fn(rng, a, b, r, alpha, T, 0, size=len(T))

recency = data[..., 0]
frequency = data[..., 1]
alive = 1 - data[..., 2]

def conditional_probability_alive_reference(
    frequency, 
    recency, 
    T
):
    log_div = (r + frequency) * np.log((alpha + T) / (alpha + recency)) + np.log(
        a / (b + np.maximum(frequency, 1) - 1)
    )

    return np.atleast_1d(np.where(frequency == 0, 1.0, expit(-log_div)))

plt.hist(conditional_probability_alive_reference(frequency, recency, T), bins=40);

../../../_images/4647b7752b50ca457542ec590df4ef04ec75f07399ba164d76f9aa05f46b258e.png

with pm.Model() as model:
    pm.Normal("a", a, sigma=0.1)
    pm.Normal("b", b, sigma=0.1)
    pm.Normal("alpha", alpha, sigma=0.1)
    pm.Normal("r", r, sigma=0.1)
    
    trace = pm.sample(chains=2)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 4 jobs)
NUTS: [a, b, alpha, r]

100.00% [4000/4000 00:02<00:00 Sampling 2 chains, 0 divergences]

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

dims = ("customer_id",)
coords = {"customer_id": range(len(frequency))}

frequency = xr.DataArray(
    frequency, 
    dims=dims,
    coords=coords,
)

recency = xr.DataArray(
    recency, 
    dims=dims,
    coords=coords,
)

T = xr.DataArray(
    T, 
    dims=dims,
    coords=coords,
)

def conditional_probability_alive(frequency, recency, T, trace):

    dims = ("customer_id",)
    coords = {"customer_id": range(len(frequency))}
    
    to_xarray = lambda array: xr.DataArray(data=array, coords=coords, dims=dims)

    frequency = to_xarray(frequency)
    recency = to_xarray(recency)
    T = to_xarray(T)
    
    a = trace.posterior["a"]
    b = trace.posterior["b"]
    alpha = trace.posterior["alpha"]
    r = trace.posterior["r"]
    
    log_div = (r + frequency) * np.log((alpha + T) / (alpha + recency)) + np.log(
        a / (b + np.maximum(frequency, 1) - 1)
    )

    return xr.where(frequency == 0, 1.0, expit(-log_div))

plt.hist(conditional_probability_alive(frequency, recency, T, trace).mean(("draw", "chain")), bins=40);

../../../_images/1900ae23f3d07812130dd40a9951a36aa20fdfd6d36c0e15d5b50921cdd7922f.png

trace.posterior["a"].mean()

<xarray.DataArray 'a' ()>
array(0.80128725)

def expected_number_of_purchases(t, frequency, recency, T):
    numerator = 1 - ((alpha + T) / (alpha + T + t)) ** (r + frequency) * hyp2f1(
        r + frequency,
        b + frequency,
        a + b + frequency - 1,
        t / (alpha + T + t),
    )
    numerator *= (a + b + frequency - 1) / (a - 1)
    denominator = 1 + (frequency > 0) * (a / (b + frequency - 1)) * (
        (alpha + T) / (alpha + recency)
    ) ** (r + frequency)
    
    return numerator/denominator

to_xarray = lambda array: xr.DataArray(array,
    coords={"customer_id": range(len(array))},
    dims=("customer_id",),
)
t = to_xarray(range(20, 40, 2))
frequency = to_xarray([1, 3, 5, 7, 9]*2)
recency = to_xarray([20, 30]*5)
T = to_xarray([25, 35]*5)

expected_number_of_purchases(t, frequency, recency, T)

<xarray.DataArray (customer_id: 10)>
array([1.42554638, 2.17059763, 3.29959668, 4.06284244, 4.72380438,
       1.73223875, 3.1645582 , 4.08454744, 5.20603817, 6.34220761])
Coordinates:
  * customer_id  (customer_id) int64 0 1 2 3 4 5 6 7 8 9