【BTYDモデル】阿部(2020)の解説
0. この記事の趣旨
この記事では BTYD モデルの一種として提案された 周期的な購買行動に対応した顧客の 生涯価値の導出と顧客維持介入戦略への応用(阿部 2020) について、その数理構造を解説することを目的としています。一応 python を用いた推定の実装についても触れるつもりではありますが、私は数値最適化や計算の高速化については十分な知見がないので、ご利用の際はご注意ください。ソースコードは以下にあります。
1. モデルの概要
本章では、以下の 2 点について解説したいと思います。
- シンプルな BTYD モデルが何をしようとしているのか
- 阿部(2020) ではどのような提案がなされているのか
本章を通じて、まずはこれらのモデルのアイデアについて理解していただければ幸いです。
1.1. BTYD モデルとは
おそらくこの記事を見ている人に BTYD モデルを全く知らない方はいらっしゃらないと思いますが、一応説明しておこうと思います。BTYD モデルとは、Buy-Till-You-Die の頭文字を取った名称で、顧客の購買行動をモデル化するために利用されるものです。最もシンプルなモデルは以下のように記述されます。
- 顧客には生存状態と死亡状態の 2 つがある。全ての顧客は最初の購買時点で生存状態となる。
- 生存状態である限り、顧客は毎日
p - 生存状態である限り、顧客は毎日
r
ただし、上記の説明は離散時間的な表現であり、実際はこれを連続時間でモデル化することが多いです。その場合、以下のように表現されます。
- 顧客には生存状態と死亡状態の 2 つがある。全ての顧客は最初の購買時点で生存状態となる。
- 生存状態において、顧客の購買の発生回数はパラメータを
\lambda - 生存状態が継続する期間はパラメータを
\mu
顧客が生存状態にあるか死亡状態にあるかは、契約を伴うデータであれば明示的に観察することができます。例えばコストコを考えてみると非常にわかりやすいです。コストコは会員のみが利用できるスーパーです。会員であれば購入が発生する可能性がありますが、退会してしまった場合は購入が発生しません。そのため、会員登録されているかどうかを生存状態か死亡状態の判断に利用できます。
しかし契約を伴う購買データはそれほど一般的ではありません。そうした状況に対応するため、BTYD モデルは非契約のデータに利用することができるように設計されています。BTYD モデルでは、その時々に顧客が生存状態にあるか死亡状態にあるかを確率的に評価しながら推論を行います。
1.2. 阿部(2020) モデルについて
上記で提案したシンプルなモデルでは顧客の購買は完全にランダムに発生しています。そのため周期的に購買が発生することが予想されるデータに対しては適切に購買行動がモデル化できない可能性があります。例えばヘアサロンの利用者は完全にランダムにヘアサロンへ行くかどうかを決めているというより、おそらく髪が伸びてくる時期に合わせて、例えば概ね 1 カ月に 1 回利用するということが考えられます。そのためこうしたデータに対してこのシンプルなモデルを当てることは少し良くないと言えるでしょう。
そこで、阿部(2020)ではこうした購買の周期性を捉えるためのモデルが提案されました。このモデルでは、前回の購買発生から
このようにモデル化することで、時間が経過するにしたがって購買が発生する確率が高まるモデルを構築することができます。以上が阿部モデルのシンプルなアイデアです。それではこの阿部モデルについて、より詳細にその数学構造を示します。
2. 阿部モデルの仮定
それではモデルの仮定を確認します。
2.1. 生存時間
生存時間
2.2. 購買金額
顧客が購入を意思決定した場合の購買金額
2.3. 購買確率
阿部モデルの概要で説明したように、前回の購買発生から
各個人のモデルの数理構造としては以上となります。また重要な前提として、これらの 3 つの確率変数は全て独立であると仮定します。
2.4. 各顧客の異質性について
阿部モデルではこの 3 つの確率的構造を前提として、各個人のパラメータを各個人のデータから推定します。これは自然なアプローチに思われると思いますが、実は BTYD モデルに関する先行研究では、各顧客のパラメータは全顧客に共通の何らかの確率密度関数に従って現れると仮定することも多いです。例えば、BG/NBD モデル では、各顧客の離脱率
こうした「個人のパラメータが確率分布に従って決定される」という考え方は直感的に感じられないかもしれませんが、このように設定する一つの理由は「ある顧客のパラメータの推定を、他の顧客のデータも参考にしながら行うことができる」ことにあります。実際阿部モデルのような顧客単位でデータを分けて各個人のパラメータを完全に別々に推定するアプローチは、顧客単位で十分なデータが得られない場合も多いため、注意が必要です。
3. 阿部モデルにおける尤度関数の導出
3.1. 観察変数の定義
さて、次に尤度関数の導出を行います。直前で述べたように阿部モデルでは個人のパラメータに全顧客に共通する確率分布を仮定しないため、ここで求めるべき尤度は各個人の購買履歴データからのみ計算できます。ある顧客が観察期間の打ち切りまでに
また、各個人は指数分布に従って死亡状態に移行します。この変数は観察できませんが、
ここからはここで用いた変数を使って、確率密度関数を計算します。詳細な導出に興味のない方は 3.4.まとめ まで飛ばしてください。
3.2. 購入金額の分離
尤度関数を記述するために、これらのデータ
と表現できます。さてここで、すべての確率変数は独立であったことを思い出しましょう。したがって購入金額とそれ以外の情報(=購買時点の情報)は独立に決定されるため
したがって、以下のように
この分離により、購買金額に関するパラメータ
購入金額については対数正規分布に従うことが仮定されているため、
3.3. 購買時点、頻度について
したがって一旦
つまり、
第 1 項
最後まで生存しているケース
まず
購買発生間隔について
また最後の購買が発生してから観測期間の打ち切りまでの期間である
期間中の状態推移について
ここでは期間の最後まで生存していることを仮定しているので、この状態遷移が発生する確率を考えます。生存時間の確率密度関数
以上のことから、
と記述できます。
観測期間の打ち切りまでのどこかで死亡状態に移行するケース
それでは次に
購買発生間隔について
先ほどと同じように
最終購買から死亡状態に移行する時点(
状態遷移について
ちょうど
これらの全ての結果を組み合わせると、
ここで、
となります。最後に
したがって
最終的な導出
ここまでで、
ただし、
3.4. まとめ
以上の議論で得られた結果を簡単にまとめておきます。確率密度関数を尤度関数を適宜読み替えていることに注意してください。
4. 推定方法
阿部モデルではこの尤度関数を用いて、（一部例外もありますが）最尤推定を行います。さてまずは扱いが簡単な購買金額に関する推定を考えましょう。
4.1. 購入金額に関するパラメータの最尤推定
この関数に対する最尤推定量は解析的に解を導くことができ、以下のようになります。
対数正規分布とは確率変数に対数を取った値が正規分布に従う分布なので、確率変数に対数を取って正規分布の最尤推定量と同じ形で表現できます。
ただし、分散の最尤推定量には一致性があるものの、不偏性がないことに注意しましょう。このことはサンプル数が十分に多い場合には良い推定量として機能するものの、サンプル数が少ない場合には推定に平均的なズレが発生してしまうことを意味します。今回は個人ごとにモデルを推定するためデータの量に限りがあるので、分散の推定量については最尤推定量ではなく不偏推定量を用いることにしましょう。
4.2. 残りのパラメータ (a,b,\mu)
先ほどのシンプルな対数正規分布を仮定して記述された尤度関数における最尤推定とは異なり、解析的に解を解くことができません。したがって数値解析の手法を用いて近似的に解を導く必要があります。そのための準備として、まずは対数尤度を計算しましょう。
論文では具体的な数値解析の手法には触れられていないため、どのような手法で実施されたかを知ることはできませんでした。そのためこの記事では、良く知られた数値解析の手法を用いて推定することにします。
5. Python を用いた実装
以下の実装は Google Colaboratory を用いて行うことを想定しています。ソースコードは以下にありますので、適宜ご利用ください。
5.1. データのダウンロード
論文ではヘアサロンのデータを用いて推定が実施されていましたが、このデータは取得することができなかったためここでは CDNOW データセットを利用します。CDNOW データセットは BTYD モデルの代表的な研究である、Fader et al (2005) にて利用されたデータセットです。
以下のコマンドでデータを取得・解凍しましょう。
!wget https://www.brucehardie.com/datasets/CDNOW_master.zip
!unzip CDNOW_master.zip
5.2. データの確認と前処理
続いてデータフレームの形式で読み込み、どのような情報が含まれているか確認します。
import pandas as pd
import warnings
warnings.simplefilter('ignore')
# テキストファイルを読み込む
# 区切り文字をスペースとして指定し、列名を明示的に設定
raw_df = pd.read_csv('CDNOW_master.txt', sep='\s+', names=['customer_id', 'date', 'quantity', 'price'])
df = raw_df.copy()
# 日付データを datetime 型に変更
df['date'] = pd.to_datetime(df['date'], format='%Y%m%d')
# 売上 (total_amount) 変数を作成
df['total_amount'] = df['quantity'] * df['price']
print('各データの型')
print(df.dtypes)
print('\nサンプル数')
print(len(df))
print('\n例')
print(df.head())
各データの型
customer_id int64
date datetime64[ns]
quantity int64
price float64
total_amount float64
dtype: object
サンプル数
69659
例
customer_id date quantity price total_amount
0 1 1997-01-01 1 11.77 11.77
1 2 1997-01-12 1 12.00 12.00
2 2 1997-01-12 5 77.00 385.00
3 3 1997-01-02 2 20.76 41.52
4 3 1997-03-30 2 20.76 41.52
顧客 ID, 日付、購入数量、購入価格、記述されているデータであることが分かります。データを確認すると、同一の顧客で同一の日付の購買データが複数レコード登場しているものが見受けられます。これはおそらく違う商品が購入された場合にレコードを分けて記述しているのだと推測されますが、私たちの分析で必要なのは各日付に合計いくら支払われたかという情報です。その情報を取得するための操作を行います。
agg_df = df.groupby(['customer_id', 'date'])['total_amount'].sum().reset_index()
print('各データの型')
print(agg_df.dtypes)
print('\nサンプル数')
print(len(agg_df))
print('\n例')
print(agg_df.head())
各データの型
customer_id int64
date datetime64[ns]
total_amount float64
dtype: object
サンプル数
67591
例
customer_id date total_amount
0 1 1997-01-01 11.77
1 2 1997-01-12 397.00
2 3 1997-01-02 41.52
3 3 1997-03-30 41.52
4 3 1997-04-02 39.08
5.3. データの選別
阿部モデルと同じようにここでは「5 回以上の購買がみられる顧客」にデータを限定します。実装は以下です。
customer_count_df = agg_df.groupby('customer_id').size().reset_index(name='count')
candidate_customers = customer_count_df[customer_count_df['count'] >= 5]['customer_id'].tolist()
print('元々のユーザー数', len(customer_count_df))
print('5回以上の購買ユーザー数', len(candidate_customers))
filtered_df = agg_df[agg_df['customer_id'].isin(candidate_customers)]
元々のユーザー数 23570
5回以上の購買ユーザー数 3835
つまり今回の推定には 23570 人のデータのうち 3835 人のデータが用いられます。
5.4. 尤度関数の記述に向けた変数の抽出
尤度関数の記述に向けて必要な変数は以下です。データセットの記述にもありますが、今回の CDNOW データセットにおける観測打ち切り時点は 1998 年 6 月 30 日です。
- 購買間隔日数
(t_1, t_2,...,t_{n-1}) - ただし、
t_1 \equiv y_2 - y_1 t_2 \equiv y_3 - y_2 t_{n-1} \equiv y_{n} - y_{n-1}
- ただし、
- 最終購買時点から観測打ち切り時点までの経過日数(
r -
1998年6月30日 - y_n
-
- 各購買における支払い金額
(s_1, s_2,..., s_n) - 観測期間開始から打ち切りまでの経過時間:
T
従って各顧客のデータを与えた際に、これらのデータを自動抽出するような関数を定義しましょう。
from typing import List, TypedDict
class Features(TypedDict):
'''尤度関数の記述に必要な変数を格納する辞書の型を定義
Args:
intervals(List[float]): 購買間隔日数のリスト
recency(float): 最終購買時点から観測打ち切り時点までの経過日数
expenditures(List[float]) 各購買における支払い金額
T(float): 観測期間開始から打ち切りまでの経過時間
'''
intervals: List[int]
recency: float
expenditures: List[float]
T: float
class CustomerFeatures(TypedDict):
customer_id: int
features: Features
def extract_features(customer_id: int, customer_df: pd.DataFrame) -> CustomerFeatures:
'''分析対象のユーザーのデータから、尤度の記述に必要な変数を取得します
Args:
df(pd.DataFrame): 以下の列が含まれるものです。
customer_id: int64
date: datetime64[ns]
total_amount: float64
'''
def _extract_intervals(customer_df: pd.DataFrame) -> List[int]:
'''各購買の間隔を計算'''
intervals = customer_df['date'].diff().dt.days.tolist()
# 先頭が nan になるので取り除く
return intervals[1:]
def _extract_recency(customer_df: pd.DataFrame) -> int:
'''最終購買時点から観測打ち切り時点までの経過日数を計算'''
observation_end = pd.to_datetime('1998-06-30')
last_purchase = customer_df['date'].max()
return (observation_end - last_purchase).days
def _extract_expenditures(customer_df: pd.DataFrame) -> List[float]:
'''各時点の購買金額を取得'''
return customer_df['total_amount'].tolist()
def _extract_T(customer_df: pd.DataFrame) -> int:
'''最初の購買が発生してから打ち切りまでの経過時間'''
observation_end = pd.to_datetime('1998-06-30')
observation_start = customer_df['date'].min()
return (observation_end - observation_start).days
intervals = _extract_intervals(customer_df)
recency = _extract_recency(customer_df)
expenditures = _extract_expenditures(customer_df)
T = _extract_T(customer_df)
result: CustomerFeatures = {
'customer_id': customer_id,
'features': {
'intervals': intervals ,
'recency': recency,
'expenditures': expenditures,
'T': T
}
}
return result
上記で定義した関数を利用して、顧客一人一人の特徴量を計算し、リスト形式で格納します。
from tqdm import tqdm
features_of_customers: List[CustomerFeatures] = []
for customer_id, customer_df in tqdm(filtered_df.groupby('customer_id')):
features = extract_features(customer_id, customer_df)
features_of_customers.append(features)
念のため、どのようなデータが取れているかいくつか確認しておきましょう。
from pprint import pprint
pprint(features_of_customers[:3])
[{'customer_id': 3,
'features': {'T': 544,
'expenditures': [41.52, 41.52, 39.08, 287.25, 83.84, 16.99],
'intervals': [87.0, 3.0, 227.0, 10.0, 184.0],
'recency': 33}},
{'customer_id': 5,
'features': {'T': 545,
'expenditures': [58.66,
13.97,
116.69999999999999,
136.64999999999998,
116.13,
52.28,
56.28,
121.41,
185.84,
121.41,
112.41],
'intervals': [13.0,
21.0,
66.0,
50.0,
16.0,
36.0,
55.0,
84.0,
4.0,
22.0],
'recency': 178}},
{'customer_id': 8,
'features': {'T': 545,
'expenditures': [9.77,
13.97,
135.87,
73.52,
356.0,
13.99,
48.92],
'intervals': [43.0, 124.0, 16.0, 136.0, 39.0, 94.0],
'recency': 93}}]
5.5. 推定の実施
さて、それでは推定を行うためのクラスを定義します。(
L-BFGS-B アルゴリズムを簡単に説明
おそらく数値最適化として最もなじみのある方法は勾配降下法かと思います。勾配降下法は一階微分の値に対し、分析者が与える学習率を掛け合わせて得られたルールに基づきパラメータを更新し続け、最適化問題を解くアプローチです。しかしこの手法では学習率の選択が非常に重要であり、不適切な設定は収束速度の低下や収束しないことにつながります。
この問題を解決するために開発された手法が ニュートン法 です。ニュートン法では一階微分の値に加えて 2 階微分（ヘッセ行列）の値まで用いることで、勾配降下法に比べてより適切にパラメータの更新が決定されます。しかしヘッセ行列の計算は計算コストが高くなるため、近似的に計算することにより計算コストを削減するアプローチが 準ニュートン法 といいます。この準ニュートン法の一つに BFGS アルゴリズム があります。
しかし、この BFGS アルゴリズムにおいて必要な近似ヘッセ行列は更新したいパラメータの数が増えれば増えるほど、メモリ要求量が大きくなることが問題になります。その問題を解決するためにヘッセ行列の近似の仕方を工夫したものが L-BFGS アルゴリズム です。
最後になりますがこの L-BFGS アルゴリズム をパラメータに制約のある最適化問題に適用できる形にしたものが L-BFGS-B アルゴリズム です。
import numpy as np
from scipy.integrate import quad
from typing import Union
from scipy.optimize import minimize, Bounds
class InitialGuess(TypedDict):
a: float
b: float
mu: float
class Estimation(TypedDict):
customer_id: int
log_m: float
omega_sq: float
a: Optional[float]
b: Optional[float]
mu: Optional[float]
success: bool
message: Optional[str]
class AbeModel:
def __init__(self):
'''各パラメータを保存するためのプロパティを作成'''
self.a: float = 1.0
self.b: float = 1.0
self.mu: float = 0.1
self.log_m: float = 1.0
self.omega_sq: float = 1.0
def _fit_log_m(self, expenditures: np.ndarray) -> float:
'''log_m の最尤推定を実施'''
return np.mean(np.log(expenditures))
def _fit_omega_sq(self, expenditures: np.ndarray, log_m: float) -> float:
'''omega^2 の不偏推定量を計算'''
degree_of_freedom = len(expenditures) - 1
return np.sum(np.square(np.log(expenditures) - log_m)) / degree_of_freedom
def _calc_G(self, intervals: Union[np.ndarray, float]) -> Union[np.ndarray, float]:
'''Gを計算'''
return np.exp(-self.a + self.b * intervals) / (1 + np.exp(-self.a + self.b * intervals))
def _calc_g(self, intervals: Union[np.ndarray, float]) -> Union[np.ndarray, float]:
'''gを計算'''
G = self._calc_G(intervals)
# gの計算で0になるのを防ぐためにepsilonを追加。ここが0になると log(g) の計算が不適切になる。
epsilon = 1e-10
g = self.b * (1 - G) * G
# gが非常に小さい値（epsilon以下）の場合は、epsilonを返す
g_with_min_value = np.maximum(g, epsilon)
return g_with_min_value
def _llh_first_term(self, intervals: np.ndarray) -> float:
'''$\sum_{i=1}^{n-1}\log{(g(t_i))}$ を計算'''
log_gs = np.log(self._calc_g(intervals))
return np.sum(log_gs)
def _llh_second_term(self, T: float) -> float:
return - self.mu * T
def _integral_term(self, r: float) -> float:
'''$\int_{0}^{r}\frac{\mu e^{-\mu x}}{1 + e^{-a+bx}}dx$ を近似計算'''
def integrand(x: float):
return (self.mu * np.exp(-self.mu * x)) / (1 + np.exp(-self.a + self.b * x))
approx_value, _ = quad(integrand, 0, r)
return approx_value
def _llh_third_term(self, r: float) -> float:
'''対数尤度の第三項を計算'''
value = 1 - self._calc_G(r) + np.exp(self.mu * r) * self._integral_term(r)
return np.log(value)
def _calc_llh(self, intervals: np.ndarray, T: float, r: float) -> float:
'''(a, b, mu) に関する対数尤度を計算'''
log_likelihood = (
self._llh_first_term(intervals) + self._llh_second_term(T) + self._llh_third_term(r)
)
return log_likelihood
def _neg_llh(self, params: np.ndarray, intervals: np.ndarray, T: float, r: float) -> float:
'''負の対数尤度関数を計算する'''
self.a, self.b, self.mu = params # パラメータを更新
return -self._calc_llh(intervals, T, r) # 対数尤度の負の値
def fit(self, intervals: np.ndarray, expenditures: np.ndarray, T: float, r: float, initial_guess: InitialGuess) -> Estimation:
'''最尤推定によるパラメータの推定'''
# まずは解析解が得られるものを計算
self.log_m = self._fit_log_m(expenditures)
self.omega_sq = self._fit_omega_sq(expenditures, self.log_m)
# 続いて数値解析でパラメータを推定
initial_guess_list = [initial_guess['a'], initial_guess['b'], initial_guess['mu']] # 初期値
bounds = Bounds([-np.inf, 1e-10, 0], [np.inf, np.inf, np.inf]) # パラメータの最小値と最大値を指定
result = minimize(self._neg_llh, initial_guess_list, args=(intervals, T, r), method='L-BFGS-B', bounds=bounds)
self.a, self.b, self.mu = result.x
return {
'log_m': self.log_m,
'omega_sq': self.omega_sq,
'a': self.a,
'b': self.b,
'mu': self.mu,
'success': result.success,
'message': result.message
}
このクラスを用いて各顧客のパラメータの推定を行いましょう。
from typing import Optional
estimations: List[Estimation] = []
for i in tqdm(range(len(features_of_customers))):
sample = features_of_customers[i]
customer_id = sample['customer_id']
intervals = np.array(sample['features']['intervals'])
expenditures = np.array(sample['features']['expenditures'])
T = sample['features']['T']
r = sample['features']['recency']
initial_guess = {'a': 1, 'b': 1, 'mu': 0.1}
model = AbeModel()
estimation = model.fit(intervals, expenditures, T, r, initial_guess)
estimation['customer_id'] = customer_id
estimations.append(estimation)
estimation_df = pd.DataFrame(estimations)
これにて推定が完了しました。お疲れさまでした。
6. まとめ
ここまで非常に長かったですが、最後までお読みくださった方は有難うございます。記述に誤りがある点などございましたら教えていただけますと幸いです。
Discussion