Gamma-Gamma Model of Monetary Value

Author

Abdullah Mahmood

Published

July 19, 2026

Source:

The Gamma-Gamma Model of Monetary Value

Code

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')

Code

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.

Code

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.

	P1X	zbar
Correlations Between Frequency & Monetary Value
P1X	1.000	0.114
zbar	0.114	1.000

Code

# 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.

shape: (9, 2)

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

Code

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

Code

# 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.

Code

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)

1 Parameter Estimation

Code

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.2503

\(q\) = 3.7440

\(\gamma\) = 15.4403

Log-Likelihood = -4055.9177

The distribution where the means have been computed across \(x → ∞\) transactions

Code

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) 
)

Code

# 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)

Code

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)

2 Computing Conditional Expectations

Code

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()

shape: (2_357, 5)

ID	P1X	zbar	Weight	E(Z\|x,zbar)
i32	u32	f64	f64	f64
1	2	2234.5	0.179996	24.653556
2	1	1177.0	0.305079	18.909032
3	0	0.0	1.0	35.170584
4	0	0.0	1.0	35.170584
5	0	0.0	1.0	35.170584
…	…	…	…	…
2353	0	0.0	1.0	35.170584
2354	5	4492.8	0.080716	44.140425
2355	0	0.0	1.0	35.170584
2356	4	3331.75	0.098899	33.500768
2357	0	0.0	1.0	35.170584