lifetimesライブラリを使って購買回数の期待値を求めてみる

## 自己紹介
初めまして、リテールでEコマースマーケティングの担当をしています。
今回、自身がまだEC担当業務に慣れていない中、マーケティング方法として何ができるか模索していた時に勉強したlifetimesライブラリを紹介したいと思います。
参考書はこちらで勉強しました。

https://www.ohmsha.co.jp/book/9784274230615/

## そもそもlifetimesライブラリとは
lifetimesライブラリはRFM（最終購買日、購買頻度、購入金額の3つの指標）を使い、LTV予測を作るために使用されます。
今回LTV予測までは行わず、lifetimesの中のBetaGeoFitterというモデルを使って購買頻度の期待値を求めるところまでにとどめてモデルを作成していきます。

## まずはデータを作成してみる。
今回はGoogle Colabを使用していきます。
必要なライブラリはこちら。
※lifetimesはcolab内にデフォルトでインストールされていないので注意

```py
!pip install lifetimes

import pandas as pd
import numpy as np
import random
from datetime import datetime, timedelta

from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score

from lifetimes import BetaGeoFitter
from lifetimes.utils import summary_data_from_transaction_data
from lifetimes.utils import calibration_and_holdout_data

import warnings
warnings.simplefilter("ignore")

```



今回生のデータを使うことはできないので、Copilotを使ってデータを作成します。
お願いしたプロンプトはこちら。

>以下の条件に沿ってCSVを作成してほしい
>・注文番号：Tから始まり6桁の数字を持つ（例T000001）
>・注文者番号：Pから始まり5桁の数字を持つ（例P00001）
>・商品名：商品1～100 ・注文数：注文数
>・注文合計金額：1商品のこの注文番号における合計金額
>・注文日：2023/01/01～2024/12/31までの2年分
>この時、同じ注文番号で複数商品を購入することもある。
>注文者のうち約3割はリピートせずに1回きりの購入だが、そのほかの注文者は2～250回の注文を行っている。
>頻度はランダムであるが多い人ほど週1回程度の頻度で買っている。
>商品名は例として3つ挙げたが、おおむね100種類ほどあり、最も注文が多いものは商品077である。
>今回ユーザーはおおむね20000人ぐらいいた。
>ファイルはおおむね400000行ぐらいになるようにしてほしい。 この条件でネットショッピングの購買データを作成してほしい。

<details><summary>この条件で作成されたコードはこちら</summary>

```py
# 商品名のリスト
products = [f"商品{i:03d}" for i in range(1, 101)]

# 商品077の選ばれる確率を高く設定
weights = [10 if product == "商品077" else 1 for product in products]

# 注文者数と注文数の設定
num_customers = 20000
repeat_customers_ratio = 0.7
repeat_customers = int(num_customers * repeat_customers_ratio)
one_time_customers = num_customers - repeat_customers

# 注文データを生成
orders = []
order_id = 1
total_rows = 0

while total_rows < 400000:
    for customer_id in range(1, num_customers + 1):
        num_orders = 1 if customer_id <= one_time_customers else random.randint(2, 250)
        for _ in range(num_orders):
            order_date = datetime(2023, 1, 1) + timedelta(days=random.randint(0, 729))
            num_items = random.randint(1, 10)
            for _ in range(num_items):
                product = random.choices(products, weights=weights, k=1)[0]
                quantity = random.randint(1, 5)
                price_per_unit = random.randint(100, 1000)
                total_price = quantity * price_per_unit
                orders.append([
                    f"T{order_id:06d}",
                    f"P{customer_id:05d}",
                    product,
                    quantity,
                    total_price,
                    order_date.strftime("%Y/%m/%d")
                ])
                total_rows += 1
                if total_rows >= 400000:
                    break
            order_id += 1
            if total_rows >= 400000:
                break
        if total_rows >= 400000:
            break

# データフレームに変換
df = pd.DataFrame(orders, columns=["注文番号", "注文者番号", "商品名", "注文数", "注文合計金額", "注文日"])

# CSVファイルに書き出し
df.to_csv("shopping_data.csv", index=False, encoding="utf-8-sig")
```
</details>

以上でデータの作成が終わり、以下のようなデータが作成されました。

|  | 注文番号 | 注文者番号 | 商品名 | 注文数 | 注文合計金額 | 注文日 |
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|0|  T000001 | P00001  |  商品027 | 4  |3564 | 2023/10/17  |
|1|  T000001 | P00001  |  商品062 |  1 | 703  | 2023/10/17  |
|2| T000001  |  P00001 | 商品045  |   2|  1822 |  2023/10/17 |
|3| T000001  | P00001  |  商品077 |  1 |  446 | 2023/10/17  |
|4|  T000001 |  P00001 |  商品086 |  5 |  3595 |  2023/10/17 |

次に、このデータをRFMの形にまとめていきます。

lifetimesにはデータさえあればRFMを作成してくれる「summary_data_from_transaction_data」があるため、今回はこちらを利用できる形に加工していきます。
必要となるデータは「顧客IDごとの日別の合計購入実績」となるので、

```py
df_pre = df.groupby(["注文者番号","注文日"])["注文合計金額"].sum().reset_index()
```

こちらを実行します。
すると、おおむね400,000ほどあったデータが65,851データにまで集約されました。

そして、こちらを使ってRFMを導きだします。

```py
df_rfm = summary_data_from_transaction_data(
        df_pre,
        "注文者番号",
        "注文日",
        observation_period_end = "2024/12/31",
        monetary_value_col = "注文合計金額")        
        
```
これを実行することで、以下のようなデータができました。

| 注文者番号 | frequency | recency | T | monetary_value |
|:-:|:-:|:-:|:-:|:-:|
| P00001  |  0.0 | 0.0  |  441.0 |  0.000000 |
| P00002  |  0.0 |  0.0 | 441.0  | 0.000000  |
| ...  | ...  | ...  | ...  | ...  |
|  P06544	 | 85.0  |  729.0 |  730.0	 | 9171.835294  |

各データについて説明すると
- frequency : 期間中のリピート回数
- recency : 期間中のアクティブ期間（最後の購入日 - 最初の購入日）
- T : 最初に購入してから期末までの期間（会員日数ともいえる）
- monetary_value : リピート時の平均購入金額

となっています。
現在、この分布がどのようになっているか、リピート回数ごとにヒストグラムで確認してみると6,000人がリピートしないデータになってしまったので、1度しか注文しない人という条件を消してもう一度データを作成しなおしました。

さて、今回はモデルを作成するため、calibration_and_holdout_dataを使います。これはtrain期間のRFMを「cal」、test期間のRFMを「holdout」として先ほどと同様のカラムを抽出させることができます。

```py
df_rfm = calibration_and_holdout_data(
        df_pre,
        "注文者番号",
        "注文日",
        calibration_period_end = "2023/12/31",
        observation_period_end = "2024/12/31",
        monetary_value_col = "注文合計金額")    
```
これを実行することで、以下のようなデータができました。

| 注文者番号 | frequency_cal | recency_cal	 | T_cal | monetary_value_cal | frequency_holdout | monetary_value_holdout | duration_holdout |
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
| P00002  | 0  |  0 |  130 |  0 |  0 | 0  | 366  |
| P00010  |  0 |  0 |  347 |  0 |  0 | 0  |  366 |
|  P00012 |  0 |  0 |  75 |  0 |  0 | 0  |  366 |
|  ... |  ... | ...  | ...  | ...  | ...  | ...  |  ... |
|  P06512 |  2 |  213 |  228 | 1739.5  | 2  | 4346  |  366 |

- duration_holdout : テスト期間の日数

そのうえで、現在のデータ上学習期間中にfrequencyが0の人は予想が難しいため、一度除外します。
```py
df_rfm = df_rfm.loc[df_rfm["frequency_cal"] > 0]

```

## モデル作成
ここまでで、データの作成が終わりましたのでいよいよモデルを作成します。

```py
BGF_model = BetaGeoFitter()

BGF_model.fit(
            df_rfm["frequency_cal"],
            df_rfm["recency_cal"],
            df_rfm["T_cal"])

print(BGF_model.summary)

```

今回、このまま実行したところ、「ConvergenceError」のようなエラーが起きました。
これはデータが収束せず、モデルの作成ができなかったために起きます。
そこでモデル内にパラメータ「penalizer_coef」を設定することで、学習を複雑化させずに回答を見つけやすく設定することができます。
デフォルトでこの値は「０」になっているので、今回は「0.01」に設定しました。
また、データそのものに関しても改めて「リピート0の人が全体の3割いて、かつfrequencyの中央値は30ぐらいになるようにする」条件を追加してデータを再修正しました。

<details><summary>再修正したデータ作成コードはこちら</summary>

```py
# Constants
NUM_USERS = 20000
MIN_ORDERS_PER_USER = 0
MAX_ORDERS_PER_USER = 250
MEDIAN_FREQUENCY = 30
PERCENT_ZERO_FREQUENCY = 0.3
MIN_DATE = datetime(2023, 1, 1)
MAX_DATE = datetime(2024, 12, 31)
NUM_PRODUCTS = 100
MOST_POPULAR_PRODUCT = "商品077"
MOST_POPULAR_PRODUCT_PROB = 0.1  # 10% chance for the most popular product
TOTAL_ROWS = 400000
OUTPUT_PATH = "shopping_data.csv"

# Generate product names
products = [f"商品{str(i).zfill(3)}" for i in range(1, NUM_PRODUCTS + 1)]

# Adjust probabilities for products
product_probs = [MOST_POPULAR_PRODUCT_PROB if p == MOST_POPULAR_PRODUCT else 
                 (1 - MOST_POPULAR_PRODUCT_PROB) / (NUM_PRODUCTS - 1) for p in products]

# Generate random order dates within the date range
def random_date(start, end):
    delta = end - start
    random_days = random.randint(0, delta.days)
    return start + timedelta(days=random_days)

# Data generation
rows = []
user_ids = [f"P{str(i).zfill(5)}" for i in range(1, NUM_USERS + 1)]

# Ensure 30% of users have frequency = 0
zero_frequency_users = random.sample(user_ids, int(NUM_USERS * PERCENT_ZERO_FREQUENCY))
remaining_users = [user for user in user_ids if user not in zero_frequency_users]

# Allocate frequencies to remaining users to achieve desired median
remaining_frequencies = [random.randint(1, MAX_ORDERS_PER_USER) for _ in range(len(remaining_users))]
remaining_frequencies.sort()
current_median = remaining_frequencies[len(remaining_frequencies) // 2]
adjustment_factor = MEDIAN_FREQUENCY - current_median
remaining_frequencies = [max(1, f + adjustment_factor) for f in remaining_frequencies]

for user_id in zero_frequency_users:
    # Users with zero frequency will not place any orders
    pass

for user_id, num_orders in zip(remaining_users, remaining_frequencies):
    for _ in range(num_orders):
        order_id = f"T{str(random.randint(1, 999999)).zfill(6)}"
        num_items = random.randint(1, 10)  # Items per order
        order_date = random_date(MIN_DATE, MAX_DATE)

        for _ in range(num_items):
            product = random.choices(products, weights=product_probs, k=1)[0]
            quantity = random.randint(1, 5)  # Quantity per item
            price_per_item = random.randint(100, 10000)  # Price per item (in yen)
            total_price = price_per_item * quantity

            rows.append({
                "注文番号": order_id,
                "注文者番号": user_id,
                "商品名": product,
                "注文数": quantity,
                "注文合計金額": total_price,
                "注文日": order_date.strftime("%Y-%m-%d")
            })

            # Limit total rows to approximately 400,000
            if len(rows) >= TOTAL_ROWS:
                break
        if len(rows) >= TOTAL_ROWS:
            break
    if len(rows) >= TOTAL_ROWS:
        break

# Create DataFrame and save to CSV
df = pd.DataFrame(rows)
df.to_csv(OUTPUT_PATH, index=False, encoding="utf-8-sig")

print(f"Generated shopping data saved to {OUTPUT_PATH}")

```
</details>

再モデル

```py
BGF_model = BetaGeoFitter(penalizer_coef=0.01)

BGF_model.fit(
            df_rfm["frequency_cal"],
            df_rfm["recency_cal"],
            df_rfm["T_cal"])

print(BGF_model.summary)
```

この結果が以下になります。

```
              coef  se(coef)  lower 95% bound  upper 95% bound
r      2.533405e+00  0.069344         2.397492         2.669318
alpha  7.240144e+01  2.206135        68.077414        76.725461
a      1.470955e-20       NaN              NaN              NaN
b      4.453423e-10       NaN              NaN              NaN
```

ではこのモデルを実行し、予測結果まで見ていきたいと思います。

```py
#予測日数を設定
duration_holdout = 366

predict = BGF_model.predict(
                        duration_holdout,
                        df_rfm["frequency_cal"],
                        df_rfm["recency_cal"],
                        df_rfm["T_cal"])

```

```py
#R2を出力
act = df_rfm["frequency_holdout"]
pred = predict
print(r2_score(act,pred))

```
結果は
```
0.5232499601079412
```
とあまりフィットしませんでした。

## まとめ
今回モデルを作成する上で生データを使えないこともあり、かなり生成AIのプロントを書くことに苦労しましたが、モデルの作成は簡単にできることが分かりました。
コードは冒頭で紹介しました参考書に準拠しています。参考書ではlifetimesの中にほかにもGGモデルといった「1回あたりの購買金額」を予測するモデルもあり、組み合わせることでLTV予測を作成することができると紹介していました。
私も引き続きLTV予測の作成に挑んでいきたいと思っています。