はじめに

購買回数/来店回数の分布は負の二項分布で表すことができるとよく言われます。
この事について、何故そうなるのかを書いていきたいと思います。

使用データ

共立出版 - データ分析プロセスに付属しているデータTafeng.csvを使用します。

データを見てみる

2000年11月のデータを使います。購買回数(来店回数)を出したいので、同じ日付かつ同じユーザーの重複を削除し、購買回数を出します。

python

df = pd.read_csv('Tafeng/Tafeng_dataset/Tafeng.csv')
df['dt'] = df['Time'].str[:10]
df['month'] = df['Time'].str[:7]
# 2000-11のみ使用
df_target = df[df['month'] == '2000-11']
# DB全体の人数
target_users = pd.DataFrame({
    'CustID': df['CustID'].unique()
})

# 同じ日付同じユーザーの重複を削除し、来店回数を出す
df_target_agg = df_target.drop_duplicates(['dt','CustID']).groupby(['CustID'], as_index=False).size()

# DB全体の人数とマージ
dfm = pd.merge(target_users, df_target_agg, how='left', on='CustID').fillna(0)

# プロット
plt.figure(figsize=(13, 8))
sns.countplot(x=dfm['size'])
plt.title(f'2000年11月 来店回数 n = {len(dfm)}')
plt.xticks(rotation=270)
plt.xlabel('購買回数')
plt.ylabel('人数')
plt.show()

購買回数	人数
0回	15506
1回	9997
2回	3526
3回	1491
4回	720
5回	361
...
22回	3
23回	1
26回	1

購買回数の分布は以上のようになりました。
0回が一番人数が多く、そこから購買回数が増えるとガクッと人数が減っていくような形です。

負の二項分布になるまでの理由

以下の2つを掛け合わせると、負の二項分布になります。

**単位期間あたりの購入回数が$\lambda$**の人が購入を単位期間中に起こす回数はポアソン分布に従う

P(X = x | \lambda) = \frac{e^{-\lambda}\lambda^x}{x!} \quad (ポアソン分布)

単位期間あたりの購入回数$\lambda$というのは、ひとによってばらつきがある。そこで**$\lambda$の分布としてガンマ分布を仮定する**
- ガンマ分布は$\alpha = 1$の時、指数分布となる
- $\beta$は尺度母数というもので、負の二項分布の導出では$\beta$は$\frac{1}{\beta}$と逆数として使用される

g(\lambda | \alpha, \beta) = \frac{\beta^{\alpha}\lambda^{\alpha - 1}e^{-\beta \lambda}}{\Gamma{(\alpha)}} \quad(ガンマ分布) \\
g(\lambda | 1, \beta) = \beta e^{-\beta \lambda} \quad (指数分布)

「単位期間あたりの購入回数がλ」とは

「単位期間」というのは自由に指定した期間の長さを1とした単位
- 今回の例では単位期間は1ヵ月
$\lambda$は1ヵ月の単位期間で購入するであろう回数の期待値

$\lambda$の値を色々と設定し、分布とランダムサンプルの結果を見てみます。

python

# ポアソン分布をプロットする関数
def poisson_exsample(lambda_, random_sample):
    x = np.arange(poisson.ppf(0.01, lambda_), poisson.ppf(0.99, lambda_))
    plt.figure(figsize=(14, 7))
    plt.subplot(121)
    plt.plot(x, poisson.pmf(x, lambda_), 'bo', ms=8, label=f'λ={lambda_}のポアソン分布')
    plt.vlines(x, 0, poisson.pmf(x, lambda_), colors='b', lw=5, alpha=0.5)
    plt.title(f'λ={lambda_}のポアソン分布 確率質量関数')
    poisson_random = poisson.rvs(lambda_, size=random_sample)
    plt.subplot(122)
    sns.countplot(poisson_random)
    plt.title(f'サンプル数={random_sample} λ={lambda_}のポアソン分布ランダム出力')
    plt.tight_layout()
    plt.savefig(f'data_folder/サンプル数={random_sample} λ={lambda_}のポアソン分布結果.png', bbox_inches='tight', pad_inches=0.3)
    plt.show()

λ = 12の場合

Aさんがいたとします。
Aさんの期待値$\lambda$ = 12の時、単位期間中に12回購入する確率は11%強です。20回や5回購入する確率は1%強ですね。
右グラフはAさんが「単位期間」を500回繰り返した時の回数分布です。

λ = 1の場合

Bさんがいたとします。
Bさんの期待値$\lambda$ = 1の時、単位期間中に1回購入する確率は35%強です。3回購入する確率は5%強ですね。
右グラフはBさんが「単位期間」を500回繰り返した時の回数分布です。

「単位期間あたりの購入回数が$\lambda$の人が購入を単位期間中に起こす回数はポアソン分布に従う」とは、
人一人にフォーカスした時、その人が購入する回数の分布です。

「λの分布としてガンマ分布を仮定する」とは

$\lambda$を持つのはAさんやBさん以外にたくさんいる
- 人それぞれ持っている$\lambda$そのものの分布
右裾が長く、非負の値を取るガンマ分布を$\lambda$の分布とする

$\alpha,\beta$の値を設定し、分布とランダムサンプルの結果を見てみます。

python

# ガンマ分布をプロットする関数
def gamma_exsample(alpha, beta, random_sample):
    x = np.arange(gamma.ppf(0.01, alpha, scale=1 / beta), gamma.ppf(0.99, alpha, scale=1 / beta))
    plt.figure(figsize=(14, 7))
    plt.subplot(121)
    plt.plot(x, gamma.pdf(x, alpha, scale=1 / beta), label=f'$\alpha$={alpha}のポアソン分布')
    plt.title(fr'$\alpha$={alpha} scale=1/{beta}のガンマ分布 確率密度関数')
    plt.subplot(122)
    gamma_random = gamma.rvs(alpha, scale=1 / beta, size=random_sample)
    plt.hist(gamma_random, bins=50)
    plt.title(fr'サンプル数={random_sample} $\alpha$={alpha} scale=1/{beta}のガンマ分布ランダム出力')
    plt.tight_layout()
    plt.savefig(f'data_folder/サンプル数={random_sample} alpha={alpha} scale=1_{beta}のガンマ分布.png', bbox_inches='tight', pad_inches=0.3)
    plt.show()

α = 1 ,β = 0.33 の場合

α = 2 ,β = 0.5 の場合

対象となる人数が500人いたとします
500人それぞれの$\lambda$の分布は上記となり、右グラフが$n$=500でランダムサンプルした結果になります。

「以下の2つを掛け合わせると、負の二項分布になる」とは

購買回数をシミュレートする流れとして
「ある人が$\lambda$を持ち、その$\lambda$を以って購入回数を$x$回起こす確率」を求めることになります。
つまりポアソン分布とガンマ分布の積という事になります。

P(X = x | \lambda)g(\lambda | \alpha, \beta) = \frac{e^{-\lambda}\lambda^x}{x!}\frac{\beta^{\alpha}\lambda^{\alpha - 1}e^{-\beta \lambda}}{\Gamma{(\alpha)}} \\
= \frac{\beta^{\alpha}\lambda^{x + \alpha - 1}e^{-\lambda(1 + \beta)}}{x!\Gamma{(\alpha)}}

$\lambda$の値を問わず、ある人が出来事をx回起こす確率を求めるため、上式を$\lambda$について積分します。

P(X = x | \alpha, \beta) = \int_{0}^{\infty}\frac{\beta^{\alpha}\lambda^{x + \alpha - 1}e^{-\lambda(1 + \beta)}}{x!\Gamma{(\alpha)}} d\lambda \\

$\lambda$と関係ないものを積分記号の外に出します。

= \frac{\beta^{\alpha}}{x!\Gamma{(\alpha)}}\int_{0}^{\infty}\lambda^{x + \alpha - 1}e^{-\lambda(1+\beta)}d\lambda

積分記号の中に$\frac{(1 + \beta)^{x + \alpha}}{\Gamma{(x + \alpha)}}$を掛け、積分記号の外に$\frac{\Gamma{(x + \alpha)}}{(1 + \beta)^{x + \alpha}}$を掛けて割り戻します

= \frac{\beta^{\alpha}}{x!\Gamma{(\alpha)}}\frac{\Gamma{(x + \alpha)}}{(1 + \beta)^{x + \alpha}}\int_{0}^{\infty}\frac{(1 + \beta)^{x + \alpha}\lambda^{x + \alpha - 1}e^{-\lambda(1+\beta)}}{\Gamma{(x + \alpha)}}d\lambda

すると、積分記号の中身は$Ga(x + \alpha , \frac{1}{1 + \beta})$のガンマ分布となります。
ガンマ分布の全確率は1なので、積分の中身は1となります。

= \frac{\beta^{\alpha}}{x!\Gamma{(\alpha)}}\frac{\Gamma{(x + \alpha)}}{(1 + \beta)^{x + \alpha}} \cdot 1 \\

残りの部分をガンマ関数と階乗パートとそれ以外のパートに分けます。

= \frac{\Gamma{(x + \alpha)}}{x!\Gamma{(\alpha)}}\frac{\beta^{\alpha}}{(1 + \beta)^{x + \alpha}} \\
= \frac{\Gamma{(x + \alpha)}}{x!\Gamma{(\alpha)}}\beta^{\alpha}\frac{1}{(1 + \beta)^{x}}\frac{1}{(1 + \beta)^{\alpha}} \\
= \frac{\Gamma{(x + \alpha)}}{x!\Gamma{(\alpha)}}\Biggl(\frac{\beta}{(1 + \beta)}\Biggr)^{\alpha} \Biggl(\frac{1}{(1 + \beta)}\Biggr)^x

ガンマ関数$\Gamma{(n)}$は$n$ = 整数の時、$\Gamma{(n)} = (n - 1)!$が成り立ちます。
つまり$x$が整数の時、$x! = \Gamma{(x + 1)}$という事になります。
そして、$x$は購入回数を指すので、整数です。

\frac{\Gamma{(x + \alpha)}}{x!\Gamma{(\alpha)}} = \frac{\Gamma{(x + \alpha)}}{\Gamma{(x + 1)}\Gamma{(\alpha)}}

ここで、$n$個から$k$個を選ぶ組み合わせの数${}_n C_k$について式をたてます。

{}_n C_k = \frac{n!}{k!(n - k)!} = \frac{\Gamma{(n + 1)}}{\Gamma{(k + 1)\Gamma{(n - k + 1)}}}

上式において$n = x + \alpha - 1, k=x$を代入すると

{}_{x + \alpha - 1} C_x = \frac{(x + \alpha - 1)!}{x!(\alpha - 1)!} = \frac{\Gamma{(x + \alpha)}}{\Gamma{(x + 1)\Gamma{(\alpha)}}}

となり、ポアソン分布ｘガンマ分布で残った式のガンマ関数と階乗のみの部分と一致します。
つまり

= \frac{\Gamma{(x + \alpha)}}{x!\Gamma{(\alpha)}}\Biggl(\frac{\beta}{(1 + \beta)}\Biggr)^{\alpha} \Biggl(\frac{1}{(1 + \beta)}\Biggr)^x \\
= {}_{x + \alpha - 1} C_x \ \Biggl(\frac{\beta}{(1 + \beta)}\Biggr)^{\alpha} \Biggl(\frac{1}{(1 + \beta)}\Biggr)^x

この式が負の二項分布となります。

負の二項分布について

確率$p$で成功するベルヌーイ試行を繰り返し行う時、$m$回成功するまでに行う試行回数が従う分布を負の二項分布と言います。

ベルヌーイ試行の定義
- 各試行は独立
- 試行結果は失敗 or 成功
- 成功確率は一定

$n - 1$回までの試行回数で$m - 1$回成功し、その次に成功する確率が確率質量関数となる為

{}_{n - 1} C_{m - 1}\ p^{m - 1}(1 - p)^{n - m} \cdot p \\
={}_{n - 1} C_{m - 1}\ p^{m}(1 - p)^{n - m}

となります。

負の二項分布のもうひとつのモデルとして、$m$回成功するまでに行う失敗回数が従う分布も負の二項分布と言います。
$n$を失敗回数、$m$を成功回数とした時、$n + m - 1$回の試行回数までに$n$回失敗し、次に成功する確率が確率質量関数となる為

{}_{n + m - 1} C_{n} \ p^{m - 1}(1 - p)^{n} \cdot p \\
{}_{n + m - 1} C_{n} \ p^{m}(1 - p)^{n}

失敗回数が従う分布の方の負の二項分布に着目します。
$n = x, m = \alpha, p = \frac{\beta}{1 + \beta}$を代入すると

{}_{x+\alpha-1} C_{x} \ \Biggl(\frac{\beta}{(1 + \beta)}\Biggr)^{\alpha} \Biggl(\frac{1}{(1 + \beta)}\Biggr)^x

となり、ポアソン分布 x ガンマ分布の計算式の結果と一致します。

人々の単位期間あたり購入回数の期待値$\lambda$がガンマ分布$Ga(\alpha, \frac{1}{\beta})$に従う時、ある人が単位期間中に$x$回購入する確率は、
確率$\frac{\beta}{1 + \beta}$で成功するベルヌーイ試行を繰り返し行う時、$\alpha$回成功するまでに行う失敗回数の分布に等しい

という事になります。

ポアソン分布 x ガンマ分布 = 負の二項分布

先に負の二項分布を描画してみます。

python

# 負の二項分布をプロットする関数
def nbinom_exsample(n, p, random_sample):
    p = round(p, 2)
    x = np.arange(nbinom.ppf(0.01, n, p), nbinom.ppf(0.99, n, p))
    plt.figure(figsize=(14, 7))
    plt.subplot(121)
    plt.plot(x, nbinom.pmf(x, n, p), 'bo', ms=8, label=f'n={n} p={p}の負の二項分布')
    plt.vlines(x, 0, nbinom.pmf(x, n, p), colors='b', lw=5, alpha=0.5)
    plt.title(f'n={n} p={p}の負の二項分布 確率質量関数')
    plt.subplot(122)
    nbinom_random = nbinom.rvs(n, p, size=random_sample)
    plt.hist(nbinom_random, bins=50)
    plt.title(f'サンプル数={random_sample} n={n} p={p}の負の二項分布 ランダム出力')
    plt.tight_layout()
    plt.savefig(f'data_folder/サンプル数={random_sample} n={n} p={p}の負の二項分布 ランダム出力.png', bbox_inches='tight', pad_inches=0.3)
    plt.show()

n = 1, p = 0.33/(1 + 0.33) (Ga(1, 1/0.33)) の場合

n = 2, p = 0.5/1 + 0.5 (Ga(2, 1/0.5)) の場合

シミュレート

「ある人が$\lambda$を持ち、その$\lambda$を以って購入回数を$x$回起こす確率」を流れに沿って描画してみます。
設定は $\alpha = 1 ,\beta = 0.33, n = 500$です

python

# 負の二項分布が生成されていく過程をgifファイルで作成する
%matplotlib notebook
fig = plt.figure(figsize = (15, 10))
a = 1
b = 0.33
n = 500
gamma_x = np.arange(gamma.ppf(0.01, a, scale=1 / b), gamma.ppf(0.99, a, scale=1 / b))
xs = []
lambdas = []

def update(i):
    plt.clf()
    lambda_ = gamma.rvs(a, scale = 1 / b, size=1)[0]
    poisson_x = np.arange(poisson.ppf(0.01, lambda_), poisson.ppf(0.99, lambda_))
    x = poisson.rvs(lambda_, size=1)[0]
    lambdas.append(lambda_)
    xs.append(x)
    plt.subplot(311)
    plt.plot(gamma_x, gamma.pdf(gamma_x, a, scale=1 / b), color='blue')
    plt.vlines(lambda_, 0, 0.3, color='darkorange', lw=10, alpha=0.8, label=f'採択されたλ={round(lambda_, 3)}')
    plt.title(fr'$\alpha$={a} scale=1/{b}のガンマ分布 確率密度関数')
    plt.legend(loc='upper right')
    plt.subplot(312)
    plt.plot(poisson_x, poisson.pmf(poisson_x, lambda_), 'bo', ms=8)
    plt.vlines(poisson_x, 0, poisson.pmf(poisson_x, lambda_), colors='b', lw=5, alpha=0.5)
    plt.vlines(x, 0, poisson.pmf(x, lambda_), color='darkorange', lw=10, alpha=0.8, label=f'採択されたx={x}')
    plt.title(f'λ={round(lambda_, 3)}のポアソン分布 確率質量関数')
    plt.legend(loc='upper right')
    plt.subplot(313)
    sns.countplot(xs, order = np.arange(0, 31))
    plt.title(f'{i + 1}人分のサンプル')
    plt.ylim(0, 30)
    plt.xlim(-0.5, 30)
    
    plt.tight_layout()

ani = animation.FuncAnimation(fig, update, interval = 500, frames = 30)
ani.save('data_folder/nbinom_simulation_fast.gif', writer='pillow')

最終的にはこのような分布になりました。
n = 1, p = 0.33/(1 + 0.33)の負の二項分布からランダム出力した分布と似ています。

データと負の二項分布を比較する

Tafeng.csvの購買回数の分布の形に近い負の二項分布のパラメータを探してみます。

scipyには負の二項分布のfit関数が用意されていないのでfit_nbinom.pyを用いました。
事前にscipyから生成した負の二項分布のパラメータに近似するかチェックしました。

python

# fit_nbinomの精度を確かめてから元データを当て嵌める
nbinom_sample = nbinom.rvs(6, 0.2, size=500)
print(fit_nbinom(nbinom_sample))
# {'size': 6.3428032398476475, 'prob': 0.2083651677991845}
nbinom_sample = nbinom.rvs(1, 0.5, size=500)
print(fit_nbinom(nbinom_sample))
# {'size': 0.926537839394276, 'prob': 0.46038281287465993}

Tafeng.csvから出した2000年11月の購買回数データを当て嵌めてみます

python

sample_fit = fit_nbinom(dfm['size'])
print(sample_fit)
# {'size': 0.9991337981786557, 'prob': 0.5029515647557096}

n = 0.999...,p = 0.502...と出ました。
このパラメータを使って同じサンプル数だけサンプルを出してみます。

python

# fit_nbinomで求めたパラメータを使って負の二項分布から元データのサンプル分ランダムサンプルする
sample_p = sample_fit['prob']
sample_n = sample_fit['size']
rsample = nbinom.rvs(sample_n, sample_p, size=len(dfm))
# sns.countplot用にDataFrameを作成
comparison_df = pd.concat([\
    pd.DataFrame({
        'size': rsample,
        'group': 'フィットデータ'
    }),
    pd.DataFrame({
        'size': list(dfm['size']),
        'group': 'オリジナルデータ'
    })
    ], axis=0
)

描画して比較してみます。

plt.figure(figsize=(13, 8))
sns.countplot(data=comparison_df, x='size', hue='group')
plt.xticks(rotation=270)
plt.legend(loc='upper right')
plt.title('フィッティングサンプルと比較')
plt.xlabel('購買回数')
plt.ylabel('人数')
plt.savefig('data_folder/フィッティングサンプルと比較.png', bbox_inches='tight', pad_inches=0.3)
plt.show()

fit_nbinom.pyで探索したパラメータの負の二項分布と形が似ています。
購買回数の分布は負の二項分布で表すことができますね。

参考にさせて頂いたページ

Pythonコード