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

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

トップ > 確率・統計 > T種類中のものをn回引いてx種類当てる確率

T種類中のものをn回引いてx種類当てる確率

確率・統計 Python

T種類中のものをn回引いてx種類当てる確率

今回は，T種類中のものを

  • n回引いてx種類当てる確率

を求める． (試行回数nは固定値，x種類が確率変数になる)

例えば，3種類の玩具付きお菓子あったときに，

  • 4回引いて1,2,3種類当てる確率

に該当する．

イメージ

幾何分布（おさらい）

幾何分布は，その事象が成功する確率$p$として， 確率変数$x$と置いたとき，$x$回試行してちょうど$x$回目に成功する確率を表現する確率分布になる．

ガチャの当たる確率などで試行回数がどのくらい必要になるかなどが計算できる．

$$ x \sim \mathrm{Geom}_{x}\left[ p \right] $$

関連：

T種類中のものをn回引いてx種類当てる確率

例：3種類中からn回引いて1,2,3種類当てる確率を考える

例えば3種類の場合で， 持っていないものを引く確率$p$から n回引いて1, 2, 3種類引く確率を考える．

例：n=1回

1回：最初はどれも持っていないので$p_{1}=1$で必ず未所持のものを入手できる． 幾何分布に応用して考えると，試行回数1，$p = p_{1}=1$の幾何分布と同じになる．

\mathrm{Geom}_{x=1}\left[ p=1 \right] = 1

$n=1$回引いて，1種類当てる確率$P(x=1|n=1)$で表すと，

P(x=1 | n=1) = \mathrm{Geom}_{x=1}\left[ p=1 \right] = 1

最初の1回は必ず100％で1種類目が手に入る

例：n=2回

2回：ここから被りが発生する． 持っていないものは残り2種類より,確率は$p_{2} = \frac{3-1}{3}=\frac{2}{3}$． 被る確率は$1-\frac{2}{3}=\frac{1}{3}$， 幾何分布で考えると，

\mathrm{Geom}_{x=1}\left[ p=2/3 \right] = \frac{2}{3}

2種類目を引く場合

2回引いて2種類引く確率は，

\begin{eqnarray} P(x=2|n=2) &=& P(x=1|n=1) \cdot \mathrm{Geom}_{x=1}\left[ p=2/3 \right] \\\\ &=& 1 \cdot \frac{2}{3} = \frac{2}{3} \\\\ \end{eqnarray}

2回引いて2種類引く確率

2回目までに1種類引く確率は，

\begin{eqnarray} P(x=1|n=2) &=& P(x=1|n=1) \cdot ( 1 - \mathrm{Geom}_{x=1}\left[ p=2/3 \right])\\\\ &=& 1 \cdot (1 - \frac{2}{3}) = \frac{1}{3} \end{eqnarray}

2回引いて1種類引く確率

例：n=3回

ここから3種類全部引ける可能性が出てくる． 1種類だけ引いているか，すでに2種類引いているかで話が変わる． すでに2種類引いて3種類目を引く確率は$p_{3} = \frac{3-2}{3} = \frac{1}{3}$. 幾何分布で考えると，

\mathrm{Geom}_{x=1}\left[ p=1/3 \right] = \frac{1}{3}

３種類目を引く場合

3回引いて全3種類引く確率は

\begin{eqnarray} P(x=3|n=3) &=& P(x=2|n=2) \cdot \mathrm{Geom}_{x=1}\left[ p=1/3 \right] \\\\ &=& 1 \cdot \frac{2}{3} \cdot \frac{1}{3} \\\\ &=& \frac{2}{9} \\\\ &=& 0.222 \cdots \end{eqnarray}

3回引いて全3種類引く確率

3回引いて2種類引く確率は，2回目に2種類目を引く場合と3回目に2種類目を引く場合の2通りある．よって，

\begin{eqnarray} P(x=2|n=3) &=& P(x=2|n=2) \cdot (1 - \mathrm{Geom}_{x=1}\left[ p=1/3 \right]) + P(x=1|n=2) \cdot \mathrm{Geom}_{x=1}\left[ p=2/3 \right] \\\\ &=& \frac{2}{3} \cdot (1 - \frac{1}{3}) + \frac{1}{3} \cdot \frac{2}{3} \\\\ &=& \frac{4}{9} + \frac{2}{9} = \frac{6}{9} \end{eqnarray}

3回引いて2種類引く確率

3回引いて1種類引く確率（つまり最初の1回以外全部外す確率）：

  • 2回引いて1種類引く確率$P(x=1|n=2)$ から次に外す確率
\begin{eqnarray} P(x=1|n=3) &=& P(x=1|n=2) \cdot \left( 1 - \mathrm{Geom}_{x=1}\left[ p=2/3 \right] \right) \\\\ &=& \frac{1}{3} \cdot \frac{1}{3} = \frac{1}{9} \end{eqnarray}

（図は省略）

例：n=4回

3種類引く確率は以下2通りある．

  • 3回目で2種類引いている確率$P(x=2|n=3)$から次に3種類目を引く場合
  • 3回目で3種類引いている確率$P(x=3|n=3)$から次に外す場合
\begin{eqnarray} P(x=3 | n=4) &=& P(x=2 | n=3) \cdot \mathrm{Geom}_{x=1}\left[ p=1/3 \right] + P(x=3 | n=3) \cdot \left(1 - \mathrm{Geom}_{x=1}\left[ p=0 \right] \right) \\\\ &=& \frac{6}{9} \cdot \frac{1}{3} + \frac{2}{9} \cdot 1 \\\\ &=& \frac{4}{9} = \frac{12}{27} \end{eqnarray}

4回引いて3種類引く確率

2種類引く確率は以下2通り

  • 3回目で1種類引いている確率$P(x=1|n=3)$から次に2種類目を引く場合
  • 3回目で2種類引いている確率$P(x=2|n=3)$から次に外す場合
\begin{eqnarray} P(x=2 | n=4) &=& P(x=1 | n=3) \cdot \mathrm{Geom}_{x=1}\left[ p=2/3 \right] + P(x=2 | n=3) \cdot \left(1 - \mathrm{Geom}_{x=1}\left[ p=\frac{1}{3} \right] \right) \\\\ &=& \frac{1}{9} \cdot \frac{2}{3} + \frac{6}{9} \cdot \frac{2}{3} \\\\ &=& \frac{2}{27} + \frac{12}{27} = \frac{14}{27} \end{eqnarray}

4回引いて2種類引く確率

1種類引く確率：

  • 3回目で1種類引いている確率$P(x=1|n=3)$から次に外す場合
\begin{eqnarray} P(x=1|n=4) &=& P(x=1 | n=3) \cdot \left( 1 - \mathrm{Geom}_{x=1}\left[ p=2/3 \right] \right) \\\\ &=& \frac{1}{9} \cdot \frac{1}{3} = \frac{1}{27} \end{eqnarray}

4回引いて1種類引く確率

5回,..., x回までに全3種類揃う確率も同様に考えることができるが， 1種類引いた時点，2種類引いた時点，3種類引いた時点で何回被るかで確率が変わるので漸化式の表現が必要になる．

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}

実装コード

厳密値

再帰で書ける． なお，メモ化しないとn=1まで辿って（再帰して）いくので時間がかかる．

import scipy as sp
import scipy.stats as stats

prob_dict = {}

def prob(N, m, T):
    if N == 1 and m == 1:
        return 1
    if N <= 0 or m <= 0:
        return 0
    if N < m:
        return 0

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

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

T=3種類でn=4回引いたときに1，2，3種類当たる確率は以下になる．

print(prob(4, 1, T=3))
> 0.037037037037037035

print(prob(4, 2, T=3))
> 0.5185185185185186

print(prob(4, 3, T=3))
> 0.44444444444444453

モンテカルロで解く

厳密値より精度は落ちるが， 再帰の計算量を気にしなくて済む．

また，容易に累積密度（x<=2種類以下などの確率）も計算できる

import scipy as sp
import scipy.stats as stats


def count_monte(T, x, n, Iter):
    numbers = sp.arange(1, T+1, 1)

    try_counts = sp.empty(Iter)
    for i in sp.arange(Iter):
        s = set(sp.random.choice(numbers, size=n))
        count = len(s)
        try_counts[i] = count
    return try_counts


def prob_monte(T, x, n, Iter=500000):
    try_counts = count_monte(T, x, n, Iter)
    return len(try_counts[try_counts == x]) / len(try_counts)

def cumsum_prob_monte(T, x, n, Iter=500000):
    try_counts = count_monte(n, x, T, Iter)
    return len(try_counts[try_counts <= x]) / len(try_counts)

> print(prob_monte(T=3, x=2, n=4))
0.51855

> print(cumsum_prob_monte(T=3, x=3, n=4))
1.0