Non contractual Customer Lifetime Value estimated probabilistically with the Beta Geometric/Negative Binomial Distribution (BG/NBD) Model

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('darkgrid')
sns.set(rc={"figure.figsize":(10, 8)})

rng = np.random.default_rng(2022)

shape, scale = 4, 11
tenure = rng.gamma(shape, scale, 10000) *2
plt.hist(tenure, density=True, edgecolor='black', bins=50)
plt.title('Histogramm - Customer Tenure')
plt.xlabel('Tenure')
plt.ylabel('Count')
plt.show()

%%time
from faker import Faker
import pandas as pd
fake = Faker(['it_IT', 'en_UK', 'fr_FR', 'de_DE', 'uk_UA'])
newcols = fake.simple_profile().keys()

profiles = pd.DataFrame(columns=list(newcols))
for i in range(10000):
    profiles.loc[i] = fake.simple_profile()

profiles['customer_id'] = profiles.index   
first_column = profiles.pop('customer_id')
profiles.insert(0, 'customer_id', first_column)
profiles.head()

# model params rounded from CDNOW sample in the paper p. 281
params=dict()

params['r'] = 0.25, 
params['alpha'] = 4.5

params['a'] = 0.8
params['b'] = 2.4

params

observation_period_end='2021-12-31'

%%time
from lifetimes.generate_data import beta_geometric_nbd_model_transactional_data

transactions = beta_geometric_nbd_model_transactional_data(tenure, params['r'], params['alpha'], params['a'], params['b'], 
                                                  observation_period_end=observation_period_end, freq='D', size=10000)
transactions.shape

trans_df = transactions.merge(profiles, left_on='customer_id', right_on='customer_id')

trans_df.head()

from lifetimes.utils import summary_data_from_transaction_data

summary = summary_data_from_transaction_data(trans_df, 'customer_id', 'date', observation_period_end=observation_period_end)

summary = pd.concat([profiles, summary], axis=1)
summary

"""from lifetimes.generate_data import beta_geometric_nbd_model
#lifetimes.generate_data.beta_geometric_nbd_model(T, r, alpha, a, b, size=1)
#Generate artificial data according to the BG/NBD model.

df = beta_geometric_nbd_model(tenure, params['r'], params['alpha'], params['a'], params['b'], size=10000)
"""
"""Parameters:	
T (array_like) – The length of time observing new customers.
alpha, a, b (r,) – Parameters in the model. See [1]_
size (int, optional) – The number of customers to generate
Returns:	
DataFrame – With index as customer_ids and the following columns: ‘frequency’, ‘recency’, ‘T’, ‘lambda’, ‘p’, ‘alive’, ‘customer_id’"""

from lifetimes import BetaGeoFitter
bgf = BetaGeoFitter(penalizer_coef=0)
bgf.fit(summary['frequency'], summary['recency'], summary['T'])
bgf.summary

from lifetimes.plotting import plot_frequency_recency_matrix
sns.set(rc={"figure.figsize":(10, 8)})
sns.set_style('dark')
plot_frequency_recency_matrix(bgf, cmap='plasma')
plt.show()

from lifetimes.plotting import plot_probability_alive_matrix

plot_probability_alive_matrix(bgf, cmap='plasma')
plt.show()

t = 7
ppcolname = 'predicted_purchases_' + str(t)
summary[ppcolname] = bgf.conditional_expected_number_of_purchases_up_to_time(t, summary['frequency'],summary['recency'], summary['T'])
summary.sort_values(by=ppcolname, ascending=False).head(5)

t = 30
ppcolname = 'predicted_purchases_' + str(t)
summary[ppcolname] = bgf.conditional_expected_number_of_purchases_up_to_time(t, summary['frequency'],summary['recency'], summary['T'])
sorted_summary = summary.sort_values(by=ppcolname, ascending=False)
sorted_summary.head()

from lifetimes.plotting import plot_period_transactions
plot_period_transactions(bgf)
plt.show()

calibration_period_end='2021-10-10'

from lifetimes.utils import calibration_and_holdout_data

summary_cal_holdout = calibration_and_holdout_data(trans_df, 'customer_id', 'date',
                                        calibration_period_end=calibration_period_end,
                                        observation_period_end=observation_period_end )
summary_cal_holdout.info()

from lifetimes.plotting import plot_calibration_purchases_vs_holdout_purchases

bgf = BetaGeoFitter(penalizer_coef=0.001)
bgf.fit(summary_cal_holdout['frequency_cal'], summary_cal_holdout['recency_cal'], summary_cal_holdout['T_cal'])
plot_calibration_purchases_vs_holdout_purchases(bgf, summary_cal_holdout)
plt.show()

bgf.fit(summary_cal_holdout['frequency_cal'], summary_cal_holdout['recency_cal'], summary_cal_holdout['T_cal'])
plot_cumulative_transactions(bgf, trans_df, 'date', 'customer_id', 280, 140);

plot_incremental_transactions(bgf, trans_df, 'date', 'customer_id',280, 140);

customer_X = sorted_summary[200:201]
customer_X

t = 30 # predict number of purchases in next t periods

individual = summary.iloc[customer_X.customer_id]
bgf.predict(t, individual['frequency'], individual['recency'], individual['T'])
# 0.0576511

from lifetimes.plotting import plot_history_alive

days_since_birth = 61 # of this customer X

sp_trans = trans_df.loc[trans_df['customer_id'] == int(customer_X.customer_id)]
plot_history_alive(bgf, days_since_birth, sp_trans, 'date', title='History of P_alive of customer_X='+str(customer_X.username) )
plt.show()

sp_trans

no_transactions = trans_df.groupby('customer_id').count().sort_values('date', ascending=False)['date']
no_transactions 

df3 = summary[summary['frequency']>0]
df3['prob_alive'] = bgf.conditional_probability_alive(df3['frequency'],df3['recency'],df3['T'])
sns.distplot(df3['prob_alive']);
plt.xlim(0,1)
plt.suptitle('Distribution of Alive Probabilities')
plt.show()

df3['churn'] = ['churned' if p < .1 else 'not churned' for p in df3['prob_alive']]
sns.countplot(df3['churn']);
plt.suptitle('Grouping of alive probabilities')
plt.show()

print('Grouping of alive probabilities')
df3['churn'][(df3['prob_alive']>=.1) & (df3['prob_alive']<.2)] = "high risk"
df3['churn'].value_counts()

## monetray value data gen

shape, scale = 3, 3.  # mean=4, std=2*sqrt(2)
monetary_value_average = rng.gamma(shape, scale, 10000) * 2
plt.hist(monetary_value_average*2, density=True, edgecolor='black', bins=50)
plt.show()

summary['monetary_value_average'] = transaction_amount

returning_customers_summary = summary[summary['frequency']>0]

returning_customers_summary.head()

returning_customers_summary[['monetary_value_average', 'frequency']].corr()

from lifetimes import GammaGammaFitter

ggf = GammaGammaFitter(penalizer_coef = 0.002)
ggf.fit(returning_customers_summary['frequency'], returning_customers_summary['monetary_value_average'])

ggf.summary

summary['Expected_cond_average_revenue'] = ggf.conditional_expected_average_profit(summary['frequency'], summary['monetary_value_average'])

summary['Expected_cond_average_revenue'].describe()

# MAPE
from sklearn.metrics import mean_absolute_percentage_error

summary["error_rev"] = summary['Expected_cond_average_revenue'] - summary['monetary_value_average']

mape = mean_absolute_percentage_error(summary['Expected_cond_average_revenue'], summary["monetary_value_average"])
print("MAPE of predicted revenues:", f'{mape:.3f}')

summary['Expected_cond_average_revenue'].head(20)

print("Expected conditional average profit: %s vs. Average profit: %s" % (
    ggf.conditional_expected_average_profit(
        summary['frequency'],
        summary['monetary_value_average']
    ).mean(),
    summary[summary['frequency']>0]['monetary_value_average'].mean()
))

# refit the BG model 
bgf.fit(summary['frequency'], summary['recency'], summary['T'])

# modelling CLV
summary['clv'] = ggf.customer_lifetime_value(
    bgf,
    summary['frequency'],
    summary['recency'],
    summary['T'],
    summary['monetary_value_average'],
    time=12, # lifetime in months
    discount_rate=0.006 
)
summary['clv'].head(10)

# describe the distribution
pd.options.display.float_format = '{:.8f}'.format
summary['clv'].describe()

ax = sns.distplot(summary['clv'], kde=False, rug=True)
ax.set_yscale('log')
plt.suptitle('Distribution of Customer Lifetime Value vs. logarithmic density')

plt.show()

summary.columns