幾何分布と負の二項分布 - 勉強したことメモ

幾何分布 (geometric distribution)

あっかりーんという存在感が薄い人がいたとします。あっかりーんが他人とすれ違うときに認識される確率が毎回30%だったとして、 3人目(3回目)に初めて認識される確率はいくらでしょうか。

この場合、2人に認識されない確率 $(1-0.3)^{2}$ と3人目に認識される確率 $0.3$ をかけて以下のように

$0.3\cdot(1-0.3)^{2}$

と計算すると思います。

これを抽象的に言うと、ベルヌーイ試行（※二項分布の記事参照）において初めてある事象（当たりなど）が起こるまでに要した回数をXとして、 $x$ 回目が当たり(確率 $p$ )で $x-1$ 回目までははずれ(確率 $1-p$ )であるので確率関数は

$f(x) = p(1-p)^{x-1} \ (x=1,2,3,···)$

となります。以下、python(+ jupyter notebook)であっかりーんが認識される確率 $p$ を0.1~0.7まで振って10回目までのレンジで分布をみてみました。

# geometric distribution
def geom_(x,p):
    return (1-p)**(x-1)*p

range_x = range(1, 11)

dict_geom= {}
for p in [n*0.1 for n in range(1,8)]:
    dict_geom[p]= [geom_(x, p) for x in range_x]

import pandas as pd

df_geom=pd.DataFrame(dict_geom)
df_geom

	0.1	0.2	0.3	0.4	0.5	0.6	0.7
0	0.100000	0.200000	0.300000	0.400000	0.500000	0.600000	0.700000
1	0.090000	0.160000	0.210000	0.240000	0.250000	0.240000	0.210000
2	0.081000	0.128000	0.147000	0.144000	0.125000	0.096000	0.063000
3	0.072900	0.102400	0.102900	0.086400	0.062500	0.038400	0.018900
4	0.065610	0.081920	0.072030	0.051840	0.031250	0.015360	0.005670
5	0.059049	0.065536	0.050421	0.031104	0.015625	0.006144	0.001701
6	0.053144	0.052429	0.035295	0.018662	0.007812	0.002458	0.000510
7	0.047830	0.041943	0.024706	0.011197	0.003906	0.000983	0.000153
8	0.043047	0.033554	0.017294	0.006718	0.001953	0.000393	0.000046
9	0.038742	0.026844	0.012106	0.004031	0.000977	0.000157	0.000014

matplotlibで可視化してみましょう。

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(range_x, df_geom, 'o')
plt.legend(df_geom.columns, title='p =')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

scipyパッケージにもあるので同じく確認してみましょう。

from scipy.stats import geom

list_geom= [geom(p).pmf(range(1,11)) for p in [n*0.1 for n in range(1,8)]]
df_geom_=pd.DataFrame(list_geom)
plt.plot(range_x, df_geom_.T, 'o')

負の二項分布 (negative binomial distribution)

幾何分布を拡張して、 $r$ 回当たりがでるまで試行を行なったとします。当たりが $r$ 回でるまでに要した回数を $X$ とすると確率関数は

$f(x) = {}_{x-1} C _{r-1} p^r(1-p)^{x-r}$

となり、これは負の二項分布(negative binomial distribution)またはPascal分布と呼ばれています。 $r=1$ のときは幾何分布になります。

また、 $r$ 回当たりがでるまでのはずれの回数を $X$ とした確率関数は

$f(x) = {}_{r+x-1} C _x p^r(1-p)^x$

となります。式の通り負の二項分布のパラメータは、当たりの確率 $p$ と当たり回数 $x$ の2つです。
(例)
あっかりーんが認識される確率 $p$ が0.3のとき、5人に認識される( $=r$ )まで出歩いたときの確率分布はどのような形になるでしょうか。以下でみてみましょう。

# r回当たるまで要した回数をX

from scipy.special import comb
def nebinom_s(x, r, p):
    return comb(x-1,r-1)*p**r*(1-p)**(x-r)

range_x2 = range(1,40)
r = 5     # 当てる回数
p = 0.3 # 当たる確率

nbi_s= [nebinom_s(x, r, p) for x in range_x2]
plt.plot(range_x2, nbi_s, 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

# r回当たるまでのはずれの回数をX

def nebinom_f(x, r, p):
    return comb(r+x-1,x)*p**r*(1-p)**x

nbi_f= [nebinom_f(x, r, p) for x in range_x2]
plt.plot(range_x2, nbi_f, 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

scipyパッケージにもあるのでそれでも試してみましょう。scipyのはデフォではずれの回数を $X$ とした確率関数の式を用いています。

from scipy.stats import nbinom

plt.plot(range_x2, nbinom.pmf(range_x2, r,p), 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

確率 $p$ が0.3のとき、当たり回数(認識した人の人数) $r$ を1~10まで振ってみましょう。

for temp_r in range(1,11):
    plt.plot(range_x2, nbinom.pmf(range_x2, temp_r, p), 'o', label=temp_r)
plt.legend(title='r =')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

あっかりーんを認識する人が増えるとそれまでの認識しなかった人が増えて行くことが分かりました(直感的にそうですね)。それと同時に分布の裾野が広がって曖昧さも増加しているようです。