問題集から一問無作為に選び出題するのを繰り返したときの，既出の問題数：あるいは，多面サイコロを繰り返し振ったときの，出た目の種類数

問題集から一問無作為に選んで出題するのを繰り返したときに，既出の問題数にまつわる確率と期待値を，とある経緯で知りたくなりました．これは，面，あるいは目の数がかなり多いサイコロを繰り返し振ったときに，これまでに出た目の種類を考えるのと等価になります．助けを借りつつ，この問題に取り組んだ結果を，ここに備忘録を兼ねて書き残しておきます．

きっかけ
問題の定式化
解答
別解：動的計画法により数値を得る
グラフ
おわりに

きっかけ

私は，あるウェブサービスを毎日のように利用しています．そのウェブサービスは，8000 問弱*1の問題集を使用し，サービスの利用者が問題を要求するたびに，その中から無作為に問題が一つ選ばれて出題されます*2．すでに問題集の収録問題数を超える，8000 あまりの要求がなされていますが，いまだに初出，すなわち以前に出されたことのない問題に遭遇することがあります．そこで，既出の問題数にまつわる確率と期待値を考えることにしました．

問題の定式化

上記の状況と，知りたい値を整理すると，次のようになります．

$k$ 個の問題を含む問題集がある．この問題集から無作為に 1 問選び出題することを繰り返す．出題を $m$ 回行ったときに，重複を含めずにちょうど $n$ 問が出題済みであるとする．

$k, m$ を定数として， $n = n_0$ となる確率 $P(n = n_0)$ を求めよ．また， $k = 8000, m = 8000$ として，全問が出題済みである確率 $P(n = k)$ を求めよ．
$k, m$ を定数として，出題済みの問題数 $n$ の期待値 $E(n)$ を求めよ．また， $k = 8000, m = 8000$ として，この具体的な値を求めよ．
$k$ を定数として， $n = k$ （全問が出題済み）となったときの， $m$ の期待値 $E'(m)$ を求めよ．また， $k = 8000$ として，この具体的な値を求めよ．

後の議論でも使いますが，この問題は以下と等価です．

$k$ 個の面（目）をもつ，偏りのないサイコロがある．このサイコロを $m$ 回振ったときに，重複を含めずにちょうど $n$ 種類の目が出たとする．

$k, m$ を定数として， $n = n_0$ となる確率 $P(n = n_0)$ を求めよ．また， $k = 8000, m = 8000$ として，すべての目がすでに出た確率 $P(n = k)$ を求めよ．
$k, m$ を定数として，出た目の種類数 $n$ の期待値 $E(n)$ を求めよ．また， $k = 8000, m = 8000$ として，この具体的な値を求めよ．
$k$ を定数として， $n = k$ （すべての目が出そろった）となったときの， $m$ の期待値 $E'(m)$ を求めよ．また， $k = 8000$ として，この具体的な値を求めよ．

解答

第 1 問： $n_0$ 問が既出である確率 $P(n = n_0)$ と，全問が既出である確率

[mixi]ガチャポン全種類ゲットするには・・・ - 数学 | mixiコミュニティにあるように，この確率は以下で与えられます．ただし， $S(\circ, \circ)$ は第 2 種 Stirling 数を意味します．

$\displaystyle \begin{equation} P(n = n_0) = \frac{\dbinom{k}{n_0} {n_0}! S(m, n_0)}{k^m} \tag{1} \label{1} \end{equation}$

また， $k = 8000, m = 8000$ のとき，全問がすでに出そろっている確率は， $P(n = k) \approx 9.88 \times 10 ^ {-3473}$ となります．この値は，Windows 10 上の数式処理システム Maxima で，2 分程度かけて得たものです．出題した問題数（のべ）が，問題集の収録問題数に達したからといって，全問が既出になることはまずないということがわかります．

第 2 問：既出の問題数の期待値 $E(n)$

それでは，のべ $m$ 問出題したとき，平均して何問が既出なのでしょうか．すなわち，期待値 $E(n)$ はいくらなのでしょうか．すぐ思いつく求め方としては，上式 $\eqref{1}$ を直接用いて
$\displaystyle \begin{equation} E(n) = \sum _{n_0 = 1} ^{k} n_0 P(n = n_0) \tag{2} \label{2} \end{equation}$
とすることが挙げられます．しかし，今回のように $k, m$ が数千と大きい場合には，数式処理システムに処理させてもなかなか答えが出ないおそれがあります．そこで，別のアプローチをとることにします．これは，A Collection of Dice Problems の pp.22-23 を応用したものです．ここでは，サイコロを繰り返し振るときの，出た目の種類数の問題を扱っており，上の言いかえと状況がちょうど合致します．問題 $i\:(i = 1, 2, \cdots, k)$ について，のべ $m$ 回の抽出および出題の間にそれが一度でも登場するときに $X_i = 1$ ，一度も登場しないときに $X_i = 0$ となる確率変数 $X_i$ を考えます．すると，
$\displaystyle \begin{equation} n = \sum _{i = 1} ^{k} X_i \tag{3}\label{3} \end{equation}$ となり，対称性から
$\displaystyle \begin{equation} E(n) = \sum _{i = 1} ^{k} E(X_i) = k E(X_1) \tag{4} \label{4} \end{equation}$ となります．いま， $m$ 回の間に，問題 1 が一度でも登場する場合の数は，すべての場合の数から，問題 1 以外の $k - 1$ 問だけが選ばれ続ける場合の数を除けばよいので
$\displaystyle \begin{equation} P(X_1 = 1) = 1 - \left(\frac{k - 1}{k}\right)^m \tag{5} \label{5} \end{equation}$
が成立します．よって
$\displaystyle \begin{equation} E(X_1 = 1) = 0 \cdot P(X_1 = 0) + 1 \cdot P(X_1 = 1) = 1 - \left(\frac{k - 1}{k}\right)^m \tag{6} \label{6} \end{equation}$
であるため，求める期待値は
$\displaystyle \begin{equation} E(n) = k \left[1 - \left(\frac{k - 1}{k}\right)^m \right] \tag{7} \label{7} \end{equation}$
となります．ここで $k = 8000, m = 8000$ とすると， $E(n) \approx 5057$ となり，3000 問程度が未出題という計算になります．

第 3 問：全問が出題済みになるとき $(k = n)$ ，のべ出題回数 $m$ の期待値 $E'(m)$

これはwikipedia:クーポンコレクター問題の一種です．求める期待値は，全 $k$ 種類の異なるクーポンがあるとき，各種類のクーポンを各 1 回以上引くまでに，クーポンを引かねばならない回数 $m$ の期待値と言い換えられます．よって
$\displaystyle \begin{equation} E'(m) = k\sum_{i = 1} ^ k \frac{1}{i} \tag{8} \label{8} \end{equation}$
となります．ここで $k = 8000$ とすると， $E'(m) \approx 76516$ となります．

別の方法として，第 2 問の結果を利用するというものがあります．式 $\eqref{7}$ を用い， $E(n) = k$ を $m$ について解くと，それが求める期待値になります*3．

別解：動的計画法により数値を得る

第 1 問では， $n_0$ 問が既出である確率 $P(n = n_0)$ を，第 2 問では既出の問題数の期待値 $E(n)$ を，それぞれ数式の形で求めました．その際に， $k, m$ が数千に達するために，得られた式 $\eqref{1}$ を用いて，それぞれの $n_0$ に対して具体的に $P(n = n_0)$ の値を求め，そこから $E(n)$ を得るのは困難であると述べました．しかし，この点について，ininsanus さんより Twitter 上で，動的計画法を用いることで，それぞれの $n_0$ に対する具体的な計算が可能になると，C++ による実装例とともに教えていただきました．これを私が Python 3 で書き換えたものを以下に示します．

#!/usr/bin/env python3


EPS = 1e-15


def used_quizs_prob(num_quizs, cur_requests, eps=EPS):
    k = num_quizs
    m = cur_requests
    dp = [0 for _ in range(k + 1)]
    dp[0] = 1.0
    for i in range(m):
        ndp = [0 for _ in range(k + 1)]
        for j in range(0, k + 1):
            if dp[j] < eps:
                continue
            ratio = j / k
            ndp[j] += ratio * dp[j]
            if j + 1 <= k:
                ndp[j + 1] += (1 - ratio) * dp[j]
        dp = ndp[:]
    exp = sum(j * dp[j] for j in range(k + 1))
    return dp, exp


def main():
    num_quizs = int(input())
    cur_requests = int(input())
    _, exp = used_quizs_prob(num_quizs, cur_requests)
    print("{:.6f}".format(exp))


if __name__ == "__main__":
    main()

ここでは，dp[k 問のうち，出題済みの問題数] := その確率 となっています．いま， $j$ 問が出題済みのとき，出題済みの問題が出る確率は $j/k$ ，新しい問題が出る確率は $1- j/k$ となります．この性質を利用して，確率の値を順次確定させることができます．

グラフ

上の動的計画法による方法を利用し， $k = m = 8000$ として $P(n = n_0)$ を図示したものが次の図です．期待値 $E(n) \approx 5057$ を中心として，急峻なピークがあることが分かります．

f:id:hamukichi_nbr:20190824194059p:plain — 図 2：既出の問題数 $n$ が $n_0$ である確率 $P(n = n_0)$ の $n_0$ 依存性．ただし $k = m = 8000$ とした．

図 2：既出の問題数 $n$ が $n_0$ である確率 $P(n = n_0)$ の $n_0$ 依存性．ただし $k = m = 8000$ とした．

また，上式 $\eqref{7}$ を用いて，期待値 $E(n)$ の $m$ 依存性を図示したものが次の図です．

f:id:hamukichi_nbr:20190824195413p:plain — 図 2：既出の問題数 $n$ の期待値 $E(n)$ の $m$ 依存性．ただし $k = 8000$ とした．

図 2：既出の問題数 $n$ の期待値 $E(n)$ の $m$ 依存性．ただし $k = 8000$ とした．

おわりに

ほんの思いつきで考えた問題ですが，ご教示をいただいたおかげで，我ながら有意義なものになったと思います．「クーポンコレクター問題」などの名前が出てくるように，あまり目新しいものではないようですが，私としては確率や期待値への理解を深めるよい機会になりました．動的計画法についてはあまり得意ではないので，上記で行っていることをより明確に言語化する，あるいは上記の類題を解くなどし，自分でも一から考えられるようにしたいものです．

*1:この問題を考えた当時の値です．現在はさらに増加しています．

*2:実際には状況はより複雑ですが，簡単のためにこうしておきます．

*3:数学的な妥当性は未証明です．

はむ吉（のんびり）の練習ノート

プログラミングとことばに関する話題を中心に，思いついたこと，試してみたこと，学んだことを，覚え書きを兼ねてまとめます．その際に役立った，技術書，参考書，辞書，機器などの紹介も行います．

問題集から一問無作為に選び出題するのを繰り返したときの，既出の問題数：あるいは，多面サイコロを繰り返し振ったときの，出た目の種類数