はしくれエンジニアもどきのメモ

情報・Web系技術・Englishの勉強メモ・備忘録です。

T種類中のものをx種類引くのに試行回数n回かかる確率

T種類中のものをx種類引くのに試行回数n回かかる確率

前回は,T種類中のものをn回引いてx種類当てる確率(n:固定値,x:確率変数)を求めたが,

cartman0.hatenablog.com

これを利用してT種類中のものをx種類引くのに試行回数n回かかる確率$P(n|x)$(つまり,前回の条件付確率の条件を反転させた確率)を求めてみる.

この確率分布は,試行回数の分布なので幾何分布や負の二項分布に近い分布になる(試行回数が増えれば増えるほどその間失敗が増えるので確率値が減衰していくような分布).より,幾何・負の2項分布のような説明に置き換えると 「 T種類中のものをx種類引くのに試行回数n回かかる確率」 = 「T種類中のものから試行回数n回目にx種類目を引く確率」になる,

今回は,

  • x種類が固定値
  • 試行回数nが確率変数になる.

例えば,T=3種類の玩具付きお菓子があったときに,

  • 全3種類中の2種類を当てるのに2,3,4,...,∞回かかる確率

に該当する.

イメージ

T種類中のものを全T種類引くのにかかる試行回数の期待値

幾何分布の期待値を応用することで,コンプリートにかかる期待値は比較的簡単に求めることができる.

これはクーポン・コレクター問題とも言うらしい?(要調査)

幾何分布の期待値は,例えば,

  • コインで表が出る場合、$p=\frac{1}{2}$より,期待値はその逆数の2回
  • サイコロで1が出る場合、$p=\frac{1}{6}$より,期待値はその逆数は6回

例として,T=3種類の場合に全種類揃えるのにかかる試行回数の期待値を考える.

1種類目:100%新しいものが引けるのでp=1の幾何分布$\mathrm{Geom}_{x}\left[p=1\right]$になる. 期待値は1回となる.

2種類目:新しいものが引ける場合はp=2/3の幾何分布$\mathrm{Geom}_{x}\left[p=2/3\right]$になる. 期待値は1/p=3/2=1.5回となる. 累計すると1種類目も含めて,$1+3/2=2.5$回となる.

3種類目:最後の1種類が引ける場合はp=1/3の幾何分布$\mathrm{Geom}_{x}\left[p=1/3\right]$になる. 期待値は3回となる. 累計すると1,2回目も含めて,$1+3/2+3/1=5.5$回となる.

つまり,全3種類の場合は5-6回引けば平均的にはコンプすることになる.

試行回数の期待値を一般化して考えると,

  • 1種類目を引くのにかかる試行回数の期待値=T/T=1
  • 2種類目を引くのにかかる期待値=T/(T-1)
  • x種類目を引くのにかかる期待値=T/(T - (x-1) )
  • T種類目を引くのにかかる期待値=T/1

これらの総和より,


\begin{eqnarray}
E[ x ]
&=& \frac{T}{T} + \frac{T}{T-1} + \cdots + \frac{T}{T - (x-1) }  + \cdots + \frac{T}{1 } \\\\
&=& \sum_{i=T}^{1} \frac{T}{i}  \\\\
&=& T \sum_{i=1}^{T} \frac{1}{i}  \\\\
\end{eqnarray}

$\sum_{i=1}^{T} \frac{1}{i$}を簡単化する公式はない?ようなので総和計算が必要になる.

pythonコードにすると以下で計算できる.

import scipy as sp

T=3
T * sp.sum(1 / sp.arange(1, T+1))
> 5.5

参考:

http://www.math.kobe-u.ac.jp/HOME/higuchi/h25kogi/13prob1ex-4-ans.pdf

T=1...10種類コンプにかかる期待値

T=1...10種類コンプにかかる期待値を表にまとめる.

種類数 試行回数の期待値
1 1.0
2 3.0
3 5.5
4 8.33
5 11.41
6 14.7
7 18.15
8 21.74
9 25.46
10 29.28

10種類は約30回試行で平均的にコンプできることになる.

T種類中からn回引いてx種類当てる確率(おさらい)

T種類中からn回引いてx種類当てる確率を一般化すると次の2通りの確率を考慮して計算できる.

  • n-1回引いてx-1種類当てている状態から次のn回目でx種類目を当てる確率
  • n-1回引いてx種類当てている状態から次のn回目で外す確率

\begin{eqnarray}
\left\{
\begin{array}{l}
P(x|n) &=& P(x-1 | n-1) \cdot \frac{T-x}{T} + P(x|n-1) \cdot \frac{x}{T} \\\\
P(x=1 | n=1) &=& 1\\\\
P(x | n) &=& 0, ( \text{for } x > n )
\end{array}
\right.
\end{eqnarray}

関連:

T種類中のものをx種類引くのに試行回数n回かかる確率

例:3種類中から2種類を引くのに1,2,...,n回かかる確率を考える

例えば3種類の場合で, 持っていないものを引く確率$p$から 2種類引くのにn回かかる確率を考える. 幾何分布的に言えば,n回目にちょうど2種類目を引く確率になる.

例:2種類目を引くのにn=2回かかる確率

2種類引くので、最低でも2回かかる.

最初の1回に1種類を引いてる状態P(x=1|n=1)から, 次の2回目に2種類目を引くと考えればいい.


\begin{eqnarray}
P(n=2|x=2)
&=& P(x=1|n=1)  \cdot \frac{2}{3} \\\\
&=& 1\cdot \frac{2}{3}
\end{eqnarray}

2種類目を引くのにn=2回かかる場合

例:2種類目を引くのにn=3回かかる場合

2回引いて1種類引いてる状態P(x=1|n=2)から, 3回目に2種類目を引いてる場合を考えればいい.


P(n=3|x=2)
&=& P(x=1|n=2)  \cdot \frac{2}{3} \\\\
&=& \frac{1}{3} \cdot \frac{2}{3}

2種類目を引くのにn=3回かかる場合

例:2種類目を引くのにn=4回かかる場合

3回引いて1種類引いてる状態P(x=1|n=3)から, 4回目に2種類目を引いてる場合を考えればいい.


P(n=4|x=2)
&=& P(x=1|n=3)  \cdot \frac{2}{3} \\\\
&=& \frac{1}{9} \cdot \frac{2}{3}

2種類目を引くのにn=4回かかる場合

5回,..., x回も同様に考えることができる.

例:3種類中から3種類を引くのに1,2,...,n回かかる確率を考える

同様に,3種類(全種類)引くのにn回かかる確率を考える.

例:3種類目を引くのにn=3回かかる確率

3種類引くので、最低でも3回かかる.被りが発生しない場合になる.

2回で2種類を引いてる状態P(x=2|n=2)から, 次の3回目に3種類目を引くと考えればいい.


\begin{eqnarray}
P(n=3|x=3)
&=& P(x=2|n=2)  \cdot \frac{1}{3} \\\\
&=& \frac{2}{3} \cdot \frac{1}{3}
\end{eqnarray}

3種類目を引くのにn=3回かかる場合

例:3種類目を引くのにn=4回かかる確率

3回で2種類を引いてる状態P(x=2|n=3)から, 次の4回目に3種類目を引くと考えればいい.


\begin{eqnarray}
P(n=4|x=3)
&=& P(x=2|n=3)  \cdot \frac{1}{3} \\\\
&=& \frac{6}{9} \cdot \frac{1}{3}
\end{eqnarray}

3種類目を引くのにn=4回かかる場合

T種類中のものをx種類当てるのにn回かかる確率の一般化

上記の例のように, n-1回引いてx-1種類引いてる状態P(x-1|n-1)から, n回目にx種類目を引く確率$(T - (x-1) ) / T$で計算できる.

これを式で表すと,


\begin{eqnarray}
P(x|n) &=& P(x-1 | n-1) \cdot \frac{T - ( x - 1 ) }{T} \\\\
\end{eqnarray}

一般化

実装コード

import scipy as sp
import scipy.stats as stats

prob_dict = {}

def prob_type_cond_try(m, n, T):
    if n == 1 and m == 1:
        return 1
    if n <= 0 or m <= 0:
        return 0
    if n < m:
        return 0

    if (m, n) in prob_dict.keys():
        return prob_dict[(m, n)]

    res = prob_type_cond_try(m, n-1, T) * (m / T) + prob_type_cond_try(m-1, n-1, T) * (1 - (m-1) / T)
    prob_dict[(m, n)] = res
    return res

def prob_try_cond_type(n, m, T):
    prob_dict = {}
    return prob_type_cond_try(m-1, n-1, T) * (1 - (m-1) / T)
prob_try_cond_type(3, m=2, T=3)
> 0.22222222222222224

3種類中の2種類引くのにx回かかる分布(ヒストグラム)を描いてみる.

import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inline

m, T = 2, 3
N = 10
n_ = sp.arange(1, N+1)
p = sp.zeros_like(n_, dtype=float)
for i in n_:
    p[i-1] = prob_try_cond_type(i, m, T)

# print(p)
plt.bar(n_, p)
plt.grid(True)
plt.xlabel("x")
plt.ylabel("p")
plt.title(f"dist $P(n | m={m}, T={T})$")
plt.show()

3種類中の2種類引くのにx回かかる分布(ヒストグラム)