マーケターのための非契約型商材のライフタイムバリュー(CLV)予測 python lifetimesパッケージを使用した実践

from google.colab import drive
drive.mount('/content/drive')

!pip install lifetimes

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import timedelta
import seaborn as sns
import os
from sklearn.metrics import mean_absolute_percentage_error
from lifetimes import BetaGeoFitter, GammaGammaFitter
from lifetimes.utils import calibration_and_holdout_data, summary_data_from_transaction_data, calculate_alive_path
from lifetimes.plotting import plot_frequency_recency_matrix, plot_probability_alive_matrix, plot_period_transactions, \
    plot_history_alive, plot_cumulative_transactions, plot_calibration_purchases_vs_holdout_purchases, plot_transaction_rate_heterogeneity, plot_dropout_rate_heterogeneity
import warnings
warnings.filterwarnings('ignore')
os.chdir(‘/content/drive/My Drive/[先ほどのファイルを置いた場所]')
print (os.getcwd())

# データ読み込み
df0 = pd.read_excel("Online Retail.xlsx")
# 中身確認
display (df0.head())

# 欠損確認
df0.isnull().sum()

# CustomerIDがないユーザのレコードを落とす
df1 = df0.copy()
df1 = df1[pd.notnull(df1["CustomerID"])]

# 返品レコードを除く (Quantity<=0除く)
df1 = df1[df1["Quantity"] > 0]
# datetimeから日時抜き出す
df1["InvoiceDate"] = pd.to_datetime(df1["InvoiceDate"]).dt.date #normalize()
# object型に変換
df1["CustomerID"] = df1["CustomerID"].astype(np.int64).astype(object)
# 確認
df1.describe(include='object').T

# 収益変数作成 revenues = quantity * unitprice
df1["Revenues"] = df1["Quantity"] * df1["UnitPrice"]
df1.head()

# (例) RFの参考値
# recency, frequency, T
dfx = df1[df1["CustomerID"] == 14527]
xmax_date = dfx["InvoiceDate"].max()
xmin_date = dfx["InvoiceDate"].min()
# recency:
print("customer minimum date:", xmin_date)
print("customer maximum date:", xmax_date)
xrec = (xmax_date - xmin_date).days
print("recency:", xrec)                                    # recency = time span between first and last purchase
# age T:
xmaxall_date = df1["InvoiceDate"].max()
print("population maximum date:", xmaxall_date)
xage = (xmaxall_date - xmin_date).days                    # age T
print("T:", xage)
# frequency:
xfreq = len(dfx[dfx["Quantity"] > 0].groupby("InvoiceDate"))-1    # frequency: periods with repeat purchases
print("frequency:", xfreq)

# determine recency, frequency, T, monetary value for each customer
df_rft = summary_data_from_transaction_data(
    transactions = df1, 
    customer_id_col = "CustomerID", 
    datetime_col = "InvoiceDate", 
    monetary_value_col = "Revenues", 
    observation_period_end = max_date, 
    freq = "D")
pd.options.display.float_format = '{:,.0f}'.format
df_rft

# axis length
max_freq = df_rft["frequency"].quantile(0.98)
max_rec = df_rft["recency"].max()
max_T = df_rft["T"].max()
# frequency
fig = plt.figure(figsize=(8, 4)) 
ax = sns.distplot(df_rft["frequency"])
ax.set_xlim(0, max_freq)
ax.set_title("frequency (days): distribution of the customers");
# recency
fig = plt.figure(figsize=(8, 4)) 
ax = sns.distplot(df_rft["recency"])
ax.set_xlim(0, max_rec)
ax.set_title("recency (days): distribution of the customers")
# T
fig = plt.figure(figsize=(8, 4)) 
ax = sns.distplot(df_rft["T"])
ax.set_xlim(0, max_T)
ax.set_title("customer age T (days): distribution of the customers")

# BG/NBD model
bgf = BetaGeoFitter(penalizer_coef=1e-06)
bgf.fit(
        frequency = df_rft["frequency"], 
        recency = df_rft["recency"], 
        T = df_rft["T"],   
        weights = None,  
        verbose = True,   
        tol = 1e-06)
pd.options.display.float_format = '{:,.3f}'.format
bgf.summary

# 検証
# frequency of repeat transactions: predicted vs actual
fig = plt.figure(figsize=(12, 12))
plot_period_transactions(
            model = bgf,  
            max_frequency = 10);

#axis length
max_freq = int(df_rft["frequency"].max())  #quantile(0.95))          
max_T = int(df_rft["T"].max())                         
max_rec = int(df_rft["recency"].max())  
t = 10 
fig = plt.figure(figsize=(7, 7)) 
plot_frequency_recency_matrix(
        model = bgf,   
        T = t,  
        max_frequency = max_freq,  
        max_recency = max_rec)

# probability that a customer has not churned (= is alive), based on 
# customer's specific recency r and frequency f
fig = plt.figure(figsize=(7, 7)) 
plot_probability_alive_matrix(
            model = bgf,  
            max_frequency = max_freq,  
            max_recency = max_rec);

# select a customer
custID = 15107
df1C = df1[df1["CustomerID"] == custID]
span_days = 365*2
# history of the selected customer: probability over time of being alive
fig = plt.figure(figsize=(20,4))
print("customer",custID,": probability of being alive over time")
plot_history_alive(
                    model = bgf, 
                    t = span_days, 
                    transactions = df1C, 
                    datetime_col = "InvoiceDate");

# monetary value > 0 のものを抽出
df_rftv = df_rft[df_rft["monetary_value"] > 0]
pd.options.display.float_format = '{:,.2f}'.format
df_rftv.describe()

# Gamma-Gamma model は frequencyとmonetary valueの相関係数が0に近いことが必要
corr_matrix = df_rftv[["monetary_value", "frequency"]].corr()
corr = corr_matrix.iloc[1,0]
print("Pearson correlation: %.3f" % corr)

# fitting the Gamma-Gamma model
ggf = GammaGammaFitter(penalizer_coef = 1e-06)
ggf.fit(
        frequency = df_rftv["frequency"],
        monetary_value = df_rftv["monetary_value"],  
        weights = None,   
        verbose = True,  
        tol = 1e-06,  
        q_constraint = True)
pd.options.display.float_format = '{:,.3f}'.format
ggf.summary

# 予測購買回数付与
def predict_purch(df, t):
        df["predict_purch_" + str(t)] = \
                bgf.predict(
                    t, 
                    df["frequency"], 
                    df["recency"], 
                    df["T"])
                    
                    
# call helper function: predict each customer's purchases over multiple time periods
t_FC = [10, 30, 60, 90]
_ = [predict_purch(df_rft, t) for t in t_FC]
pd.options.display.float_format = '{:,.1f}'.format
print("predicted number of purchases for each customer over next t days:")
df_rft

# 予測継続率付与
prob_alive = bgf.conditional_probability_alive(
        frequency = df_rft["frequency"], 
        recency = df_rft["recency"], 
        T = df_rft["T"])
df_rft["prob_alive"] = prob_alive
pd.options.display.float_format = '{:,.2f}'.format
df_rft.describe()

# 予測購買額付与
exp_avg_rev = ggf.conditional_expected_average_profit(
        df_rftv["frequency"],
        df_rftv["monetary_value"])
df_rftv["exp_avg_rev"] = exp_avg_rev
df_rftv["avg_rev"] = df_rftv["monetary_value"]
df_rftv["error_rev"] = df_rftv["exp_avg_rev"] - df_rftv["avg_rev"]
mape = mean_absolute_percentage_error(exp_avg_rev, df_rftv["monetary_value"])
print("MAPE of predicted revenues:", f'{mape:.2f}')
pd.options.display.float_format = '{:,.3f}'.format
df_rftv.head()

# CLVを計算する
DISCOUNT_a = 0                # 年間割引率 ここでは0でよい
LIFE = 12                        # 12か月後のCLV
clv = ggf.customer_lifetime_value(
        transaction_prediction_model = bgf,  
        frequency = df_rftv["frequency"],  
        recency = df_rftv["recency"],  
        T = df_rftv["T"],  
        monetary_value = df_rftv["monetary_value"],  
        time = LIFE,    
        freq = "D",                          
        discount_rate = DISCOUNT_a)
df_rftv.insert(0, "CLV", clv)             # expected customer lifetime values
df_rftv.describe().T

plt.hist(df_rftv['CLV'],bins=1000)
plt.xlim(0,30000)
plt.grid(True)
plt.title('Predict CLV ')

# MとCLVの関係
fig, axes = plt.subplots(1,2,figsize=(14,7))
axes[0].scatter(df_rftv['monetary_value'],df_rftv['CLV'])
axes[0].set_title('M * CLV')
# FとCLVの関係
axes[1].scatter(df_rftv['frequency'],df_rftv['CLV'])
axes[1].set_title('F * CLV')