MarTech Farm

#31│BG/NBDモデルをSQLで近似する ― 「この顧客はまだ生きているか」

RFMの限界と確率論的顧客分析

RFM分析は強力だが、根本的な限界がある。

企業によって定義は様々あるが「最終購買から180日経った顧客」を「休眠顧客」と呼ぶ。しかしこの顧客は本当に「まだいつか戻ってくる可能性がある顧客」なのか、それとも「もう二度と買わない顧客（チャーン済み）」なのか、RFMだけでは判断できない。

その問いに確率論的な答えを与えるのが BG/NBD モデル（Beta-Geometric / Negative Binomial Distribution model）だ。
2004年にFader・Hardie・Leeが提唱した手法で、非契約型ビジネス（= 退会手続きなしに黙って去るECや通販）における顧客の「生死」と「購買頻度」を同時にモデル化する。

このモデルが与えてくれる答えは2つだ。

  1. この顧客がまだアクティブである確率（P(alive)）
  2. 次のN期間でこの顧客が購買する期待回数（期待購買数）

Python の lifetimes ライブラリを使えば自動計算できるが、今回はその推定ロジックをSQLで近似実装する。なぜか。データが Treasure Data にある場合、わざわざ Python にエクスポートしなくてもSQLだけで P(alive) を計算し、そのままセグメント配信に使えるからだ。

BG/NBDモデルの直感的な理解

数式の前にモデルの考え方を整理する。

BG/NBD は2つのプロセスを仮定する。

プロセス1：購買のプロセス（Negative Binomial Distribution）

アクティブな顧客は「ランダムなタイミング」で購買する。顧客ごとに購買頻度（λ）が異なり、その分布がガンマ分布に従う、と仮定する。

プロセス2：離脱のプロセス（Beta-Geometric）

顧客は「購買のたびに一定確率 p で離脱する」と仮定する。
この離脱確率 p も顧客ごとに異なり、ベータ分布に従う、と仮定する。

つまり「たくさん買う顧客」と「あまり買わない顧客」が混在し、「すぐ離脱する顧客」と「長く続く顧客」が混在する現実を、2つの分布の組み合わせで表現する。

観察データから何が推定できるか

各顧客について観察できるのは次の3つだ。

  • x：観察期間中の購買回数
  • t_x：最終購買日から観察開始日までの経過時間（観察期間内での最終購買タイミング）
  • T：観察期間の長さ（最初の購買日から観察終了日まで）

この3つだけから、「この顧客がまだアクティブである確率」を推定できる。

使用するテーブル

-- orders テーブル
-- order_id     : 注文ID
-- customer_id  : 顧客ID
-- order_date   : 注文日（DATE型）
-- total_amount : 注文金額
-- status       : 'completed' / 'cancelled'

STEP 1 ― BG/NBD の入力パラメータを顧客ごとに計算する

まず各顧客の x（購買回数）、t_x（最終購買タイミング）、T（観察期間）を計算する。

観察期間は「初回購買日から基準日まで」とする。

-- STEP1: BG/NBD 入力パラメータの計算

WITH customer_rfm AS (
    SELECT
        customer_id,
        MIN(order_date)   AS first_order_date,   -- 初回購買日
        MAX(order_date)   AS last_order_date,    -- 最終購買日
        COUNT(order_id)   AS x,                  -- 総購買回数

        -- T: 観察期間（初回購買日から基準日まで、週単位）
        DATE_DIFF('day', MIN(order_date), DATE '2024-12-31') / 7.0  AS T,

        -- t_x: 最終購買タイミング（初回購買日から最終購買日まで、週単位）
        DATE_DIFF('day', MIN(order_date), MAX(order_date)) / 7.0    AS t_x
    FROM orders
    WHERE
        status     = 'completed'
        AND order_date >= DATE '2022-01-01'  -- 観察開始日
        AND order_date <= DATE '2024-12-31'  -- 観察終了日（基準日）
    GROUP BY customer_id
    HAVING COUNT(order_id) >= 1  -- 少なくとも1回購買
)
SELECT
    customer_id,
    first_order_date,
    last_order_date,
    x,
    ROUND(T, 2)    AS T_weeks,
    ROUND(t_x, 2)  AS t_x_weeks,
    -- 初回購買が最終購買の場合、t_x = 0（1回のみ購買）
    CASE WHEN x = 1 THEN 0 ELSE ROUND(t_x, 2) END  AS t_x_safe
FROM customer_rfm
ORDER BY x DESC, T DESC;

出力イメージ

customer_idfirst_order_datelast_order_datexT_weekst_x_weeks
C0042022-03-122024-11-1518143.4136.1
C0012022-08-052024-10-208125.3116.0
C0232023-01-202024-09-08497.781.3
C0892024-10-012024-10-01113.00.0

STEP 2 ― モデルパラメータを決定する

BG/NBD モデルには4つのパラメータがある。

α, β  ： 購買頻度の分布（ガンマ分布）のパラメータ
r, s  ： 離脱確率の分布（ベータ分布）のパラメータ

本来はこれらを最尤推定で求める（Python の lifetimes ライブラリが自動計算する）。SQLでの近似では、事前に Python で推定したパラメータ値を定数として使うことが現実的だ。

典型的な通販データでの参考値（実際のデータで必ずキャリブレーションする）

-- モデルパラメータ（事前にPythonで推定した値を定数として定義）
-- 実データで推定した値に必ず置き換えること

WITH model_params AS (
    SELECT
        0.243  AS r,    -- 購買頻度分布の形状パラメータ
        4.414  AS alpha, -- 購買頻度分布のスケールパラメータ
        0.793  AS a,    -- 離脱確率分布のパラメータ
        2.425  AS b     -- 離脱確率分布のパラメータ
)

Pythonでのパラメータ推定（SQLの外で実行する）

from lifetimes import BetaGeoFitter
bgf = BetaGeoFitter()
bgf.fit(df['frequency'], df['recency'], df['T'])
print(bgf.params_)
# {'r': 0.243, 'alpha': 4.414, 'a': 0.793, 'b': 2.425}

このパラメータを上記CTEに貼り付けてSQLで使う。

STEP 3 ― P(alive) を計算する（BG/NBDの核心）

BG/NBD モデルで「顧客がまだアクティブである確率」は次の式で求まる。

P(alive | x, t_x, T) = 1 / (1 + (a / (b + x - 1)) × A)

ここで：
A = (α + T)^(r + x) / (α + t_x)^(r + x)
  × Σ_{j=0}^{x-1} [ C(x-1, j) × (-1)^j × (b + x - 1) / (b + x - 1 - j) ]
  × [(α + T)/(α + t_x)]^(-j)  （超幾何関数の近似）

これを厳密に SQL で実装するのは複雑すぎるため、実務でよく使われる簡略近似式を使う。

顧客がアクティブな確率の直感的な近似として、最も単純かつ実用的なのは「最終購買日から現在までの沈黙期間が、その顧客の通常の購買間隔と比べてどれだけ異常か」に基づくアプローチだ。

-- STEP3: P(alive) の実用的近似計算

WITH model_params AS (
    SELECT 0.243 AS r, 4.414 AS alpha, 0.793 AS a, 2.425 AS b
),
customer_rfm AS (
    SELECT
        customer_id,
        COUNT(order_id)                                           AS x,
        DATE_DIFF('day', MIN(order_date), MAX(order_date)) / 7.0 AS t_x,
        DATE_DIFF('day', MIN(order_date), DATE '2024-12-31') / 7.0 AS T,
        DATE_DIFF('day', MAX(order_date), DATE '2024-12-31') / 7.0 AS recency_weeks  -- 最終購買からの経過週数
    FROM orders
    WHERE status = 'completed'
    AND order_date BETWEEN DATE '2022-01-01' AND DATE '2024-12-31'
    GROUP BY customer_id
),
palive_calc AS (
    SELECT
        c.customer_id,
        c.x,
        ROUND(c.t_x, 2)             AS t_x,
        ROUND(c.T, 2)               AS T,
        ROUND(c.recency_weeks, 2)   AS recency_weeks,

        -- BG/NBD P(alive) 近似式
        -- 参考：Fader et al. 2005 の式(7)を数値近似
        -- δ = a / (b + x - 1) ← 離脱確率のオッズ
        -- λ = (r + x) / (α + T) ← 期待購買率
        ROUND(
            1.0 / (
                1.0 + (
                    (p.a / NULLIF(p.b + c.x - 1, 0))
                    * POWER((p.alpha + c.T) / NULLIF(p.alpha + c.t_x, 0.001), p.r + c.x)
                )
            )
        , 4)  AS p_alive,

        -- 期待購買頻度（週あたり）
        ROUND((p.r + c.x) / (p.alpha + c.T), 4)  AS lambda_hat

    FROM customer_rfm   c
    CROSS JOIN model_params  p
    WHERE c.x >= 1  -- 1回以上購買
)
SELECT
    customer_id,
    x           AS purchase_count,
    t_x         AS recency_at_last_weeks,
    T           AS observation_weeks,
    recency_weeks AS weeks_since_last_order,
    p_alive,
    lambda_hat  AS weekly_purchase_rate,

    -- P(alive) のセグメント分類
    CASE
        WHEN p_alive >= 0.80 THEN 'アクティブ（確実）'
        WHEN p_alive >= 0.60 THEN 'アクティブ（高確率）'
        WHEN p_alive >= 0.40 THEN 'グレーゾーン'
        WHEN p_alive >= 0.20 THEN '離脱リスク高'
        ELSE                      'チャーン済みの可能性大'
    END  AS alive_segment,

    -- 次の12週間での期待購買数
    ROUND(p_alive * lambda_hat * 12, 2)  AS expected_orders_next_12w
FROM palive_calc
ORDER BY p_alive DESC, lambda_hat DESC;

出力イメージ

customer_idpurchase_countweeks_since_lastp_alivealive_segmentexpected_orders_next_12w
C004186.70.9821アクティブ（確実）1.82
C001810.40.9134アクティブ（確実）0.87
C023424.10.6412アクティブ（高確率）0.31
C089113.00.4821グレーゾーン0.12
C112378.20.0841チャーン済みの可能性大0.01

C004 は18回購買・最終購買から6.7週間と直近で、P(alive)=0.98。ほぼ確実にまだアクティブだ。C112 は3回購買したが78週間（約1.5年）沈黙しており、P(alive)=0.08。ほぼチャーンと見てよい。

STEP 4 ― RFM×P(alive) でセグメントを精緻化する（完成版）

RFMスコアと P(alive) を組み合わせることで、従来のRFMだけでは区別できなかった顧客を分離できる。

-- STEP4: RFM + P(alive) による精緻なセグメント分類

WITH model_params AS (
    SELECT 0.243 AS r, 4.414 AS alpha, 0.793 AS a, 2.425 AS b
),
customer_rfm AS (
    SELECT
        customer_id,
        COUNT(order_id)                                           AS x,
        DATE_DIFF('day', MIN(order_date), MAX(order_date)) / 7.0 AS t_x,
        DATE_DIFF('day', MIN(order_date), DATE '2024-12-31') / 7.0 AS T,
        DATE_DIFF('day', MAX(order_date), DATE '2024-12-31') / 7.0 AS recency_weeks,
        SUM(total_amount)                                         AS monetary
    FROM orders
    WHERE status = 'completed'
    AND order_date BETWEEN DATE '2022-01-01' AND DATE '2024-12-31'
    GROUP BY customer_id
),
with_palive AS (
    SELECT
        c.*,
        ROUND(
            1.0 / (
                1.0 + (
                    (p.a / NULLIF(p.b + c.x - 1, 0))
                    * POWER((p.alpha + c.T) / NULLIF(p.alpha + c.t_x, 0.001), p.r + c.x)
                )
            )
        , 4)  AS p_alive,
        ROUND((p.r + c.x) / (p.alpha + c.T), 4)  AS lambda_hat
    FROM customer_rfm   c
    CROSS JOIN model_params  p
)
SELECT
    customer_id,
    x,
    ROUND(recency_weeks, 1)  AS weeks_since_last,
    monetary,
    ROUND(p_alive, 3)        AS p_alive,
    ROUND(p_alive * lambda_hat * 52, 1)  AS predicted_orders_next_year,
    ROUND(p_alive * lambda_hat * 52 * (monetary / NULLIF(x, 0)), 0)
                             AS predicted_revenue_next_year,

    -- 4象限マトリクス（P(alive) × 購買金額）
    CASE
        WHEN p_alive >= 0.7 AND monetary >= 50000 THEN 'VIP顧客（高価値・高生存）'
        WHEN p_alive >= 0.7 AND monetary <  50000 THEN '育成顧客（低価値・高生存）'
        WHEN p_alive <  0.7 AND monetary >= 50000 THEN '緊急ウィンバック対象（高価値・低生存）'
        ELSE                                           '自然消滅待ち（低価値・低生存）'
    END  AS strategic_segment
FROM with_palive
ORDER BY predicted_revenue_next_year DESC NULLS LAST;

出力イメージ

customer_idxweeks_since_lastmonetaryp_alivepredicted_revenue_next_yearstrategic_segment
C004186.7284,6000.982308,000VIP顧客（高価値・高生存）
C001810.498,4000.913112,000VIP顧客（高価値・高生存）
C088542.1182,0000.41275,000緊急ウィンバック対象（高価値・低生存）
C019128.228,4000.95132,000育成顧客（低価値・高生存）
C112378.298,0000.0848,000自然消滅待ち（低価値・低生存）

「緊急ウィンバック対象」の C088 は累計18万円購買しているが、P(alive)=0.41 と離脱リスクが高い。RFMだけ見ると「高M・低R」で通常のウィンバック対象だが、P(alive) まで計算することで「今すぐ手を打たないとほぼ戻らない」という判断が数値に基づいてできる。

STEP 5 ― P(alive) を使ったコホート別の「生存率推移」

全顧客の P(alive) を月次バッチで計算し、コホート別の平均 P(alive) の推移を追うと、「どのコホートが速く死んでいくか」が定量化できる。

-- STEP5: 初回購買月（コホート）別の平均P(alive)推移

WITH model_params AS (SELECT 0.243 AS r, 4.414 AS alpha, 0.793 AS a, 2.425 AS b),
customer_rfm AS (
    SELECT
        customer_id,
        DATE_TRUNC('month', MIN(order_date))                       AS cohort_month,
        COUNT(order_id)                                            AS x,
        DATE_DIFF('day', MIN(order_date), MAX(order_date)) / 7.0   AS t_x,
        DATE_DIFF('day', MIN(order_date), DATE '2024-12-31') / 7.0 AS T
    FROM orders
    WHERE status = 'completed'
    AND order_date BETWEEN DATE '2022-01-01' AND DATE '2024-12-31'
    GROUP BY customer_id
),
with_palive AS (
    SELECT
        c.cohort_month,
        c.customer_id,
        ROUND(
            1.0 / (
                1.0 + (
                    (p.a / NULLIF(p.b + c.x - 1, 0))
                    * POWER((p.alpha + c.T) / NULLIF(p.alpha + c.t_x, 0.001), p.r + c.x)
                )
            )
        , 4)  AS p_alive
    FROM customer_rfm c CROSS JOIN model_params p
)
SELECT
    cohort_month,
    COUNT(customer_id)             AS cohort_size,
    ROUND(AVG(p_alive), 3)         AS avg_p_alive,
    ROUND(APPROX_PERCENTILE(p_alive, 0.5), 3)  AS median_p_alive,
    SUM(CASE WHEN p_alive >= 0.7 THEN 1 ELSE 0 END)  AS high_survival_count,
    ROUND(
        SUM(CASE WHEN p_alive >= 0.7 THEN 1 ELSE 0 END) * 100.0
        / NULLIF(COUNT(customer_id), 0)
    , 1)  AS high_survival_pct
FROM with_palive
GROUP BY cohort_month
ORDER BY cohort_month;

出力イメージ

cohort_monthcohort_sizeavg_p_alivehigh_survival_pct
2022-01-018420.51241.2%
2022-06-019180.49839.8%
2023-01-011,2410.62152.4%
2023-06-011,1840.63854.1%
2024-01-011,8920.78468.3%

2022年初頭のコホートは avg P(alive) が 0.51 と低下しており、2年以上経過して離脱が進んでいる。2023年以降のコホートは P(alive) が高く、近年獲得した顧客の生存率が改善していることが分かる。

実務での運用ヒント

① パラメータは必ず自社データで推定する

参考値として示したパラメータ（r=0.243 など）は Fader らの論文でのサンプル値だ。業種・商材・購買サイクルによって大きく変わる。Python の lifetimes ライブラリを使い、自社の過去データ（最低でも12〜24ヶ月）でパラメータを推定してから SQL に組み込む。

# Pythonでのパラメータ推定
from lifetimes import BetaGeoFitter
bgf = BetaGeoFitter(penalizer_coef=0.0)
bgf.fit(rfm_df['frequency'], rfm_df['recency_weeks'], rfm_df['T_weeks'])
print(bgf.params_)

② 購買1回のみの顧客の扱い

x = 1（1回購買のみ）の場合、t_x = 0 となり分母にゼロが生じるリスクがある。NULLIF(p.alpha + c.t_x, 0.001) で小さな値でのゼロ除算を防ぐか、初回のみ購買者は別途「新規顧客リスト」として扱う設計が安全だ。

③ 毎月の再計算サイクルを組む

P(alive) は新しい購買があるたびに更新される。少なくとも月次バッチで再計算し、Treasure Data の monthly_palive テーブルとして保存する。このテーブルをセグメント配信の基礎データとして使うことで、常に最新の「生存確率ベースのリスト」が維持できる。

④ 厳密な実装ではなく「BG/NBDの精神」を使う

学術的に厳密な BG/NBD の P(alive) 計算式は超幾何関数を含む複雑なものだ。今回の実装は近似だが、「沈黙期間が長いほど P(alive) が下がる」という本質的な動作は正確に再現している。精度を求めるなら Python の lifetimes ライブラリで計算し、結果をTDにロードする方法が最も正確だ。

まとめ

BG/NBD モデルの SQL 近似実装をまとめる。

  1. 顧客ごとに x（購買回数）・t_x（最終購買タイミング）・T（観察期間）を週単位で計算する
  2. Python で推定した4パラメータ（r, α, a, b）を定数として CTE に定義する
  3. P(alive) の近似式を POWER と NULLIF で実装する
  4. RFM の M（購買金額）と P(alive) を組み合わせて「高価値・高生存」「高価値・低生存」の4象限に分類する
  5. 「次の1年の期待購買数」= P(alive) × λ̂ × 52 で将来価値を推計する

RFM が「現在の状態のスナップショット」なら、BG/NBD は「将来の行動の確率的予測」だ。この2つを組み合わせることで、「今すぐ手を打つべき高価値・低生存顧客」を精確に特定し、施策のROIを最大化できる。

契約内容によるが、treasure dataでもdigdagでpythonを起動させることもできるのでまずはSQLが面倒くさいならばPythonでやるのもいいかもしれない。

おすすめ記事

#19│ポイント残高をトランザクションログから再計算するSQL

MarTech Farmをもっと見る

今すぐ購読し、続きを読んで、すべてのアーカイブにアクセスしましょう。

続きを読む