[Python]仕事で使えそうな確率まとめ – かものはしの分析ブログ

本来、正月の企画物としてアイデアを出していたものの、作成する時間がとれなかったのですが、日常的にもっと確率に従って生活できたら面白いなと思って作ってみることにしました。

二項分布

やりたいこと
オフィスにおいて、チームごとにデスク周辺の掃除当番を毎日乱数で決めているとして、1ヶ月で自分がどれくらい掃除をやらなくてはいけないのかを知りたい。あるいは、1回の掃除で済む確率とか、10回掃除する確率とかを知りたい。

これは、成功確率p（掃除当番になる確率 = 1÷チームの人数）でのN回の独立なベルヌーイ試行となり、自分が掃除当番になった場合X=1で、そうでない場合はX=0となる。
N回ベルヌーイ試行した際のXの合計は二項分布従うので、N回試行して、Xがkになる確率は以下のように示される。（チームの人数はHとする。）

$$ P(Y_N = k ) = {}_N C _k \left (\frac{1}{H} \right )^k \left ( 1- \frac{1}{H} \right )^{N-k} \\
(0 \leq k \leq N )$$

コード

from scipy.stats import binom
import matplotlib.pyplot as plt

H = 10 # チーム人数
D = 20 # 営業日

n, p = D,  1/H

# 平均、分散、歪度、尖度
mean, var, skew, kurt = binom.stats(n, p, moments='mvsk')

print('平均：',mean)

fig, ax = plt.subplots(1, 1)

x = range(D + 1)
ax.plot(x, binom.pmf(x, n, p, loc=0), 'bo', ms=8, label='binom pmf')
ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5);
ax.set_xlabel("count"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability mass function");

from scipy.stats import binom

import matplotlib.pyplot as plt

H = 10 # チーム人数

D = 20 # 営業日

n, p = D, 1/H

# 平均、分散、歪度、尖度

mean, var, skew, kurt = binom.stats(n, p, moments='mvsk')

print('平均：',mean)

fig, ax = plt.subplots(1, 1)

x = range(D + 1)

ax.plot(x, binom.pmf(x, n, p, loc=0), 'bo', ms=8, label='binom pmf')

ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5);

ax.set_xlabel("count"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability mass function");

平均では20営業日で2回は掃除することになりますが、4回以上掃除しないといけない確率が10%以上はあることがわかります。掃除が嫌いな人は、苦しみの上限をイメージできると楽なのではないかと思うので、確率に従い掃除当番を甘んじて受け入れましょう。

参考情報
scipy.stats.binom

matplotlib.pyplot.grid

多項分布

やりたいこと
先ほどのように乱数で掃除当番を割り当て、それをランキングにし、一位を掃除大臣として称号を付与するものとする。10名いるもののうち、上位のメンバーA・メンバーBとで累積の掃除回数がそれぞれ10回、6回としたもとで次の月に掃除大臣が入れ替わる確率を計算したい。

コード

import math

H = 10 # チーム人数
D = 20 # 営業日

# あるメンバーがa回当たり、あるメンバーがb回当たる確率
def calc_prob_multinomial(a, b, H=H, D=D):
    return math.factorial(D) / (math.factorial(a) * math.factorial(b) * math.factorial(D - a - b)) * (1 / H)**(a) * (1 / H)**(b) * (1 - 2 / H )**(D - a - b)

base_a = 10 # メンバーaの現在の掃除回数
base_b = 6 # メンバーbの現在の掃除回数
gap = base_a - base_b # 差分

prob_sum = 0.0
for i in range(D + 1):
    for j in range(1, D + 1):
        b = gap + i + j  # 差分に1を足すことで負けるケースを計算（b = jにしたら確率の合計は1になった。）
        if b + i  < D + 1: # 合計が月数を超えないようにする。
            prob_sum = prob_sum + calc_prob_multinomial(i, b, H=H, D=D) # 確率を足す。
        
print('sum:',prob_sum)

import math

H = 10 # チーム人数

D = 20 # 営業日

# あるメンバーがa回当たり、あるメンバーがb回当たる確率

def calc_prob_multinomial(a, b, H=H, D=D):

return math.factorial(D) / (math.factorial(a) * math.factorial(b) * math.factorial(D - a - b)) * (1 / H)**(a) * (1 / H)**(b) * (1 - 2 / H )**(D - a - b)

base_a = 10 # メンバーaの現在の掃除回数

base_b = 6 # メンバーbの現在の掃除回数

gap = base_a - base_b # 差分

prob_sum = 0.0

for i in range(D + 1):

for j in range(1, D + 1):

b = gap + i + j # 差分に1を足すことで負けるケースを計算（b = jにしたら確率の合計は1になった。）

if b + i < D + 1: # 合計が月数を超えないようにする。

prob_sum = prob_sum + calc_prob_multinomial(i, b, H=H, D=D) # 確率を足す。

print('sum:',prob_sum)

計算したところ、1.2%くらいになった。このケースにおいてAさんがBさんに次の月に掃除大臣の座を奪われる可能性は低いらしい。

参考情報
scipy.special.comb

誕生日の問題

やりたいこと
あまりにも有名な確率の話で、誕生日のパラドックスとも呼ばれています。誕生日が同じ人が部署内にいる確率が0.5を超える人数は何人か？

コード

import matplotlib.pyplot as plt

N = 50 # メンバーの数の上限
D = 365 # 誕生日の数

def calc_prob_birth(N=N, D=D):
    prob = 1.0
    for i in range(1, N):
        prob = (1 - i /D) * prob

    return 1 - prob

fig, ax = plt.subplots(1, 1)

x = range(2 ,N + 1)
probability = list(map(lambda x:calc_prob_birth(N=x) , x))
ax.plot(x, probability, 'bo', ms=8)
ax.hlines(y=0.5, xmin= 0, xmax=N, color='r');
ax.set_xlabel("the number of classmate"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability");

import matplotlib.pyplot as plt

N = 50 # メンバーの数の上限

D = 365 # 誕生日の数

def calc_prob_birth(N=N, D=D):

prob = 1.0

for i in range(1, N):

prob = (1 - i /D) * prob

return 1 - prob

fig, ax = plt.subplots(1, 1)

x = range(2 ,N + 1)

probability = list(map(lambda x:calc_prob_birth(N=x) , x))

ax.plot(x, probability, 'bo', ms=8)

ax.hlines(y=0.5, xmin= 0, xmax=N, color='r');

ax.set_xlabel("the number of classmate"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability");

参考情報
誕生日のパラドックス
 matplotlib.pyplot.hlines

ファーストサクセス分布

やりたいこと
アルバイトの面接を受けるとして、その合格率が1/30だとして、何回くらい面接しても合格できないものだろうか。苦しみの上限を知りたい。20回、30回面接を受けた時に合格できてない確率はどれくらいだろうか。何回くらい面接を受けたら合格できていない確率を1%未満に下げれるだろうか。
成功確率がpのベルヌーイ試行を独立に何回も行うとき、初めて成功数までに失敗した回数をXとして、初めて成功した際の試行回数をYとした場合、その確率は$ k \geq 1 $として、

$$ P(Y = k) = P(X=k-1)=p(1-p)^{k-1} $$

と表すことができる。これをパラメータpのファーストサクセス分布と呼ぶ。（幾何分布とも呼ぶ。）

コード

J = 1/ 30 # 面接合格率

def calc_prob_first_success(J=J, k=k):
    return (J) * (1- J )**(k-1)

fig, ax = plt.subplots(1, 1)

x = range(1 , 10000) # 無視できるくらい伸ばしておく。
probability = list(map(lambda x:calc_prob_first_success(k=x,) , x))
probability_sum = list(map(lambda x:1 - sum(probability[:x]) , x))
length = 90 # 90回の失敗までを見る。

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)
ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("the probability we can not success on the number of failure");

J = 1/ 30 # 面接合格率

def calc_prob_first_success(J=J, k=k):

return (J) * (1- J )**(k-1)

fig, ax = plt.subplots(1, 1)

x = range(1 , 10000) # 無視できるくらい伸ばしておく。

probability = list(map(lambda x:calc_prob_first_success(k=x,) , x))

probability_sum = list(map(lambda x:1 - sum(probability[:x]) , x))

length = 90 # 90回の失敗までを見る。

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)

ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("the probability we can not success on the number of failure");

20回面接をしても合格できない確率はおよそ50%。30回面接しても合格できない確率はおよそ36%。合格できない確率が1%を切るのは面接回数が135回になった時。このレベルを追い求めると心が折れる面接回数なのかもしれない。

参考情報
弱点克服大学生の確率・統計

負の二項分布

やりたいこと
ある営業がリンゴを売る仕事をしておりタウンページで飲食店に電話をかけているとする、商談の予約が取れる確率(p)が0.05として、商談予約が5(n)回起きるまでに失敗がk回生じる確率について知りたい。いったい、この営業は飲食店にどれくらい商談をお断りをされる覚悟で臨めばいいのか。

このような事象は負の二項分布に従うので、
$$ P(X = k) = {}_{n – 1 + k} C _k p^{n-1}(1-p)^k $$
という確率に従う。
なお、この分布の期待値は
$$ E(X) = \frac{n(1-p)}{p} $$
に従う。

コード

import matplotlib.pyplot as plt
import math

prob = 1/ 20 # リンゴ営業の商談予約できる確率
n = 5 # 営業の目標である商談予約の回数

# 負の二項分布の確率密度関数
def calc_prob_nb(prob=prob, K=K, n=n):
    return math.factorial(K+n-1) / (  math.factorial(n-1) * math.factorial(  K+n-1 - (n-1) )) *(prob)**n * ( 1- prob )**K

fig, ax = plt.subplots(1, 1)

x = range(1 , 200)
probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x))
length =200

ax.plot(x[:length], probability[:length], 'bo', ms=8)
ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1);
ax.set_title("Probability until the success happens on five times");

import matplotlib.pyplot as plt

import math

prob = 1/ 20 # リンゴ営業の商談予約できる確率

n = 5 # 営業の目標である商談予約の回数

# 負の二項分布の確率密度関数

def calc_prob_nb(prob=prob, K=K, n=n):

return math.factorial(K+n-1) / ( math.factorial(n-1) * math.factorial( K+n-1 - (n-1) )) *(prob)**n * ( 1- prob )**K

fig, ax = plt.subplots(1, 1)

x = range(1 , 200)

probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x))

length =200

ax.plot(x[:length], probability[:length], 'bo', ms=8)

ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1);

ax.set_title("Probability until the success happens on five times");

75回くらいの失敗で確率自体のピークがきている。期待値だと95回だが、確率のピークはもっと手前にある。

累積確率を知りたい場合はこちらのコード

fig, ax = plt.subplots(1, 1)

x = range(1 , 200) # 無視できるくらい伸ばしておく。
probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x))
probability_sum = list(map(lambda x:sum(probability[:x]) , x))
length =200

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)
ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1);
ax.set_title("Cumulative probability until the success happens on five times");

fig, ax = plt.subplots(1, 1)

x = range(1 , 200) # 無視できるくらい伸ばしておく。

probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x))

probability_sum = list(map(lambda x:sum(probability[:x]) , x))

length =200

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)

ax.set_xlabel("the number of failure"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1);

ax.set_title("Cumulative probability until the success happens on five times");

80%のところで125回なので、5人に一人は期待値の95回を30回くらいは超えてしまってもおかしくないということになる。そして、5人に一人くらいは50〜60回で達成できるラッキーな人もいる。

参考情報
弱点克服大学生の確率・統計

ポアソン分布

やりたいこと
先ほど、営業が飲食店に売ったリンゴの中に芋虫が住み付いているケースが500回に1回の頻度で発生し、営業と顧客の間でトラブルに発展するとする。年間に1200の飲食店と取引を行うが、トラブルはどれくらい発生するものなのか。
顧客とのトラブルが発生する回数がYだとした際に、ポアソン分布に従うとした場合、
$$P(Y = k) = \frac{\lambda^k}{k!} e^{-\lambda} $$
となる。ここで、
$$ \lambda = 1200 \times \frac{1}{500} $$
とおく。
なお、期待値と分散はともに
$$ E(Y) = V(Y) = \lambda $$
となる。

コード

import matplotlib.pyplot as plt
import math
import numpy as np

lamda = 1200 * (1 /500 )

def calc_prob_poisson(lamda=lamda, K=K):
    return lamda**K / math.factorial(K) * np.exp(-lamda)

fig, ax = plt.subplots(1, 1)

x = range(0 , 10)
probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x))
length =10

ax.plot(x[:length], probability[:length], 'bo', ms=8)
ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1);

import matplotlib.pyplot as plt

import math

import numpy as np

lamda = 1200 * (1 /500 )

def calc_prob_poisson(lamda=lamda, K=K):

return lamda**K / math.factorial(K) * np.exp(-lamda)

fig, ax = plt.subplots(1, 1)

x = range(0 , 10)

probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x))

length =10

ax.plot(x[:length], probability[:length], 'bo', ms=8)

ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1);

トラブルが2回の時が確率が一番高い。0回で済む確率も0.1に近いくらいはある。

fig, ax = plt.subplots(1, 1)

x = range(0 , 10)
probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x))
probability_sum = list(map(lambda x:sum(probability[:x]) , x))
length = 10

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)
ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability");
ax.grid(color='b', linestyle='-', linewidth=0.1);
ax.set_title("Cumulative probability");

fig, ax = plt.subplots(1, 1)

x = range(0 , 10)

probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x))

probability_sum = list(map(lambda x:sum(probability[:x]) , x))

length = 10

ax.plot(x[:length], probability_sum[:length], 'bo', ms=8)

ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability");

ax.grid(color='b', linestyle='-', linewidth=0.1);

ax.set_title("Cumulative probability");

5人に一人は4回以上のトラブルに巻き込まれる気の毒な営業がいそうですね。まぁ、リンゴの中は開けてみないとわかりませんから、顧客トラブルがあっても恨みっこなしで頑張りましょう。

参考情報
弱点克服大学生の確率・統計

フィギュア集め問題

やりたいこと
現場猫フィギュアのガチャポンに関心があるとする。N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率はどれくらいだろうか？（どの現場猫のフィギュアも出現確率は同じものとする。）
また、N+1回ガチャポンを回したとき、N種類の全部が揃う確率はどれくらいだろうか？あるいは、m回ガチャポンを回してそろったフィギュアの数の期待値は？そして、Y回目のガチャポンではじめてN種類の現場猫の全部揃う、Yの期待値は？

・N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率：
$$\frac{N!}{N^N}$$

・N+1回ガチャポンを回したとき、N種類の全部が揃う確率：
$$ N \times \frac{(N+1)!}{2} \times \frac{1}{N^{N+1}} $$

・m回ガチャポンを回してそろったフィギュアの種類$ X_m $の期待値：
m回ガチャポンを回してそろったフィギュアの種類$ X_m $を以下のように表す。
$$X_m = Z_1 + Z_2 + \dots + Z_N$$
m回ガチャポンを引いて、当たらないということは
$$ P(Z_i = 0) = \left ( \frac{N – 1}{N} \right )^m $$
ということになる。
これをもとに$ Z_i $の期待値を計算すると

$$ E(Z_i) = 1 \times \left ( 1 – \left ( \frac{N – 1}{N} \right )^m \right ) + 0 \times \left ( \frac{N – 1}{N} \right )^m \\
= 1 – \left( \frac{N-1}{N} \right )^m $$
よって、$ X_m $の期待値は
$$ E(X_m) = E( Z_1 + \dots + Z_N) \\
= N \left ( 1 – \left( \frac{N-1}{N} \right )^m \right )$$

となる。

・Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値：
i-1種類からi種類になるまでのガチャポンの回す回数を$ T_i $として、その階差がファーストサクセス分布に従うとする。
$$ T_i – T_{i-1} \sim F_s \left ( \frac{N – i + 1}{N} \right ) $$
その場合、Yの期待値は以下のようになる。

$$ E(Y) = E( T_N ) \\
= E( (T_{N} – T_{N-1}) + (T_{N-1} – T_{N-2}) + \dots + (T_{2} – T_{1}) + T_{1} ) \\
= E \left ( F_s \left ( \frac{1}{N} \right ) + F_s \left ( \frac{2}{N} \right ) \dots + F_s \left ( \frac{N-1}{N} \right ) + 1 \right ) \\
= N \left ( 1 + \frac{1}{2} + \dots + \frac{1}{N} \right )
$$

コード

import matplotlib.pyplot as plt
import math

N = 5 # 現場猫フィギュアの種類
m = 5 # 回すガチャポンの回数m

# N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率
print(math.factorial(N) / (N ** N))


# N+1回ガチャポンを回したとき、N種類の全部が揃う確率
print(N * math.factorial(N+1) / 2 * (1/N**(N+1)) )


# m回ガチャポンを回してそろったフィギュアの数のXの期待値
print(N * (1 - ((N-1)/N)**m  ))


# Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値
def harmonic_series_sum(n): 
    i = 1
    s = 0.0
    for i in range(1, n+1): 
        s = s + 1/i; 
    return s; 

print(N * (harmonic_series_sum(N)))

import matplotlib.pyplot as plt

import math

N = 5 # 現場猫フィギュアの種類

m = 5 # 回すガチャポンの回数m

# N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率

print(math.factorial(N) / (N ** N))

# N+1回ガチャポンを回したとき、N種類の全部が揃う確率

print(N * math.factorial(N+1) / 2 * (1/N**(N+1)) )

# m回ガチャポンを回してそろったフィギュアの数のXの期待値

print(N * (1 - ((N-1)/N)**m ))

# Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値

def harmonic_series_sum(n):

i = 1

s = 0.0

for i in range(1, n+1):

s = s + 1/i;

return s;

print(N * (harmonic_series_sum(N)))

N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率：
0.0384
N+1回ガチャポンを回した時、N種類の現場猫フィギュア全部が揃う確率：
0.1152
m回ガチャポンを回して揃った現場猫フィギュアの種類Xの期待値：
3.3615
Y回目のガチャポンではじめてN種類の現場猫フィギュアが揃う際のYの期待値：
11.4166

参考情報
弱点克服大学生の確率・統計

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル