Source:
import polars as pl
import numpy as np
from sklearn.neighbors import KernelDensity
from scipy.optimize import minimize
from scipy.special import gammaln
from utils import CDNOW, modified_silverman
import altair as alt
from IPython.display import display_markdown
alt.renderers.enable("html")RendererRegistry.enable('html')cdnow = CDNOW(master=False, calib_p=273)
# For the Gamma-Gamma model, we need to filter out customers who have made only one purchase.
rfm_data = cdnow.rfm_summary().filter(pl.col('P1X') > 0)
rfm_data_array = rfm_data.select('P1X', 't_x', 'T', 'zbar').collect().to_numpy()
x = rfm_data_array[:,0] # frequency
zbar = rfm_data_array[:,3] / 100 # monetary value
t_x = rfm_data_array[:,1]
T = rfm_data_array[:,2]The Gamma-Gamma model assumes that there is no relationship between the monetary value and the purchase frequency. We can check this assumption by calculating the correlation between the average spend and the frequency of purchases.
corr_data = rfm_data.select('P1X', 'zbar').collect()
(
corr_data.corr()
.with_columns(pl.Series(corr_data.columns).alias("index"))
.style.tab_header(title="Correlations Between Frequency & Monetary Value")
.tab_stub(rowname_col="index")
.fmt_number(decimals=3)
)
# The value of this correlation is close to 0.11, which in practice is considered low enough to proceed with the model.| Correlations Between Frequency & Monetary Value | ||
|---|---|---|
| P1X | zbar | |
| P1X | 1.000 | 0.114 |
| zbar | 0.114 | 1.000 |
# Descriptive statistics of the average spend per repeat transaction
summary = rfm_data.select('zbar').with_columns(pl.col('zbar') / 100).describe()
summary
# We note that the distribution of observed individual means is highly skewed to the right.| statistic | zbar |
|---|---|
| str | f64 |
| "count" | 946.0 |
| "null_count" | 0.0 |
| "mean" | 35.077848 |
| "std" | 30.283506 |
| "min" | 2.99 |
| "25%" | 15.76 |
| "50%" | 27.54 |
| "75%" | 41.79 |
| "max" | 299.63381 |
Probability density estimate of the sample
m = np.arange(2.5, 301, 2.5) # Average transaction value range
# Apply log transformation for boundary correction
m_log = np.log(m)
zbar_log = np.log(zbar)
bw = modified_silverman(zbar_log)
print('Kernel Smoothing Bandwidth:', bw) Kernel Smoothing Bandwidth: 0.18800626075684287# Estimate the probability density function
# Method 1 - Using sklearn
kde = KernelDensity(kernel='gaussian', bandwidth=bw).fit(zbar_log.reshape(-1,1))
log_density = kde.score_samples(m_log.reshape(-1,1))
f = np.exp(log_density) / m # Transform the density back to the original scale
# Method 2 - Using statsmodels
# import statsmodels.api as sm
# kde = sm.nonparametric.KDEUnivariate(zbar_log)
# kde.fit(kernel='gau', bw=bw)
# f_log = kde.evaluate(m_log)
# f = f_log / m The distribution of average spend per (repeat) transaction across the 946 individuals who made a repeat transaction in the calibration period. Each customer’s average is computed across a (typically very) small number of transactions.
act_dist_plot = (
alt.Chart(pl.DataFrame({'Average Transaction Value (z)': m, 'f(z)': f}))
.mark_line().encode(
x=alt.X(
'Average Transaction Value (z)',
axis=alt.Axis(values=np.arange(0, 301, 50),
labelExpr='"$"+datum.value')),
y=alt.Y('f(z)', scale=alt.Scale(domain=[0, 0.04]))
)
)
act_dist_plot.properties(
width=500,
height=400,
title='Observed distribution of average transaction values across customers'
).configure_view(stroke=None).configure_axisY(grid=False).configure_axisX(grid=False) def gammagamma(x, zbar, guess={'p': 0.01, 'q': 0.01, 'gamma': 0.01}):
def log_likelihood(param):
p, q, gamma = param[0], param[1], param[2]
ll = (
gammaln(p*x+q) -
gammaln(p*x) -
gammaln(q) +
q*np.log(gamma) +
(p*x-1)*np.log(zbar) +
(p*x)*np.log(x) -
(p*x+q)*np.log(gamma+x*zbar)
)
return -np.sum(ll)
bnds = [(1e-6, np.inf) for _ in range(3)]
return minimize(log_likelihood, x0=list(guess.values()), bounds=bnds, method='L-BFGS-B')
res = gammagamma(x=x, zbar=zbar)
p, q, gamma = res.x
ll = res.fun
# Sample Parameters
# p = 6.24983547654959
# q = 3.7441106896737
# gamma = 15.4423198312514
display_markdown(f'''$p$ = {p:0.4f}
$q$ = {q:0.4f}
$\\gamma$ = {gamma:0.4f}
Log-Likelihood = {-ll:0.4f}''', raw=True)\(p\) = 6.2480
\(q\) = 3.7444
\(\gamma\) = 15.4483
Log-Likelihood = -4055.9177
The distribution where the means have been computed across \(x → ∞\) transactions
zeta = np.arange(300) + 1
f_zeta = (p * gamma)**q * zeta**(-q-1) * np.exp(-p*gamma/zeta) / np.exp(gammaln(q))
(
alt.Chart(pl.DataFrame({'Unobserved mean transaction value (ζ)': zeta, 'f(ζ)': f_zeta}))
.mark_line().encode(
x=alt.X(
'Unobserved mean transaction value (ζ)',
axis=alt.Axis(values=np.arange(0, 301, 50),
labelExpr='"$"+datum.value')),
y=alt.Y('f(ζ)', scale=alt.Scale(domain=[0, 0.04]))
).properties(
width=500,
height=400,
title='Distribution of the (unobserved) mean transaction value (ζ)'
).configure_view(stroke=None).configure_axisY(grid=False).configure_axisX(grid=False)
)# compute the density of average transaction value
# how many people with each level of x?
repeat_trans_dist = (
cdnow.rfm_summary()
.group_by('P1X')
.agg(pl.len().alias('Count'))
.sort('P1X')
.collect()
.to_numpy()
)
nx = repeat_trans_dist[1:, 1]
x_trans = repeat_trans_dist[1:, 0]
# compute the density of zbar for each x
y = np.arange(300) + 1
x_trans, y = np.meshgrid(x_trans,y)
a1 = gammaln(p*x_trans+q)-gammaln(p*x_trans)-gammaln(q)
a2 = q*np.log(gamma)
a3 = (p*x_trans-1)*np.log(y)
a4 = (p*x_trans)*np.log(x_trans)
a5 = (p*x_trans+q)*np.log(gamma+y*x_trans)
g1 = np.exp(a1+a2+a3+a4-a5)
# compute the weighted average
g = np.dot(nx, g1.T) / np.sum(nx)est_dist_plot = (
alt.Chart(pl.DataFrame({'Average Transaction Value (z)': np.arange(300) + 1, 'f(z)': g}))
.mark_line(strokeDash=[4,4]).encode(
x=alt.X(
'Average Transaction Value (z)',
axis=alt.Axis(values=np.arange(0, 301, 50),
labelExpr='"$"+datum.value')),
y=alt.Y('f(z)', scale=alt.Scale(domain=[0, 0.04]))
)
)
chart = act_dist_plot + est_dist_plot
chart.properties(
width=500,
height=400,
title='Observed versus theoretical distribution of average transaction value across customers'
).configure_view(stroke=None).configure_axisY(grid=False).configure_axisX(grid=False) E_Z = p*gamma/(q-1)
ce = (
cdnow.rfm_summary()
.select('ID', 'P1X', 'zbar')
.with_columns(((q - 1)/(p * pl.col('P1X') + q - 1)).alias('Weight'))
.with_columns(
(pl.col('Weight')*E_Z+(1-pl.col('Weight'))*pl.col('zbar')/100)
.alias('E(Z|x,zbar)')
)
)
ce.collect()| ID | P1X | zbar | Weight | E(Z|x,zbar) |
|---|---|---|---|---|
| i32 | u32 | f64 | f64 | f64 |
| 1 | 2 | 2234.5 | 0.180074 | 24.654499 |
| 2 | 1 | 1177.0 | 0.305191 | 18.911551 |
| 3 | 0 | 0.0 | 1.0 | 35.170302 |
| 4 | 0 | 0.0 | 1.0 | 35.170302 |
| 5 | 0 | 0.0 | 1.0 | 35.170302 |
| … | … | … | … | … |
| 2353 | 0 | 0.0 | 1.0 | 35.170302 |
| 2354 | 5 | 4492.8 | 0.080755 | 44.140022 |
| 2355 | 0 | 0.0 | 1.0 | 35.170302 |
| 2356 | 4 | 3331.75 | 0.098946 | 33.500827 |
| 2357 | 0 | 0.0 | 1.0 | 35.170302 |