食玩問題について

食玩問題とは

食品玩具、略して食玩にはおまけの玩具がついているが、いくつ買うとその玩具を全種類集められるだろうか。おまけが2種類以上であれば、いくつ買っても全部揃う確率が1になることはないが、買えば買うほど1に近づいていく。その購入数と確率の関係を考える問題である。

食玩問題を定式化をしてみると、全 $n$ 種($n \geq 1$)のおまけがついている食玩があり、そのおまけは常に等確率で出現するとする。このとき $r$ 個($r \geq n$)買ったとき、全種類揃っている確率 $p_n(r)$ を求める問題となる。

結論

早速、結論を書く (導出過程に興味がある方は、下まで読んでもらえると僕が喜びます) 。

$m = r - n$ とする。その $m$ は必要最低限の個数から余分に買った個数である。このとき、求めたい確率 $p_n(r)$ は以下のようになる。

$$p_n(r) = \frac{(n + m)!}{n^{(n + m)}} \sum^n_{k=0} {}_m\mathrm{T}_k \,\, n^{\underline{k}}$$

ここで、 $n^{\underline{k}}$ は降べき乗を表しており、 $k = 0$ のときは $1$ で、$k = 3$ のときは $n(n-1)(n-2)$ のようになる。また ${}_m\mathrm{T}_k$ は次を満たすものとする。

$${}_0\mathrm{T}_0 = 1, \,\,\,\, {}_m\mathrm{T}_k = 0 \,\,\,\, (m \lt 0 \,\,\,\, \mathrm{or} \,\,\,\, k \lt 0 \,\,\,\, \mathrm{or} \,\,\,\, m \lt k)$$

$${}_m\mathrm{T}_k = \frac{1}{m + k} ({}_{m - 1}\mathrm{T}_{k - 1} + k \,\, {}_{m - 1}\mathrm{T}_k) \,\,\,\, (m, k \geq 1)$$

プログラミングしてみた

プログラム

上記の式をもとに計算できるプログラムを作成した。全 $n$ 種のおまけで、購入数 $n$ 個から順番に確率 $p = a / b$ を初めて超えるまで計算するプログラムとした。
コマンドライン引数で、 $n, a, b$ を引き渡すプログラムとなる。

Coupon.py

from fractions import Fraction
from math import factorial
import sys


class Coupon():
    def __init__(self, n: int):
        if n <= 0:
            raise Exception("n is wrong")

        self.n = n

    @staticmethod
    def getNextCofficients(m: int, cofficients: list[Fraction]) -> list[Fraction]:
        return list(map(
            lambda same, diff, index: (diff + Fraction(index, 1) * same) / Fraction(m + index, 1),
            cofficients,
            [Fraction(0, 1)] + cofficients[:-1],
            range(len(cofficients)),
        ))

    def getProbabilityList(self, maxProb: Fraction) -> None:
        if maxProb >= Fraction(1, 1):
            raise Exception("p is wrong")

        # 最初の係数 [1, 0, 0, ..., 0] を作成
        cofficients = [Fraction(0, 1) for _ in range(self.n + 1)]
        cofficients[0] = Fraction(1, 1)

        # n に関する降べき乗 [1, n, n(n-1), ..., n!] を作成
        variables = []
        v = Fraction(1, 1)
        for i in range(self.n + 1):
            variables.append(v)
            v *= Fraction(self.n - i)

        # 最初の先頭の係数を計算
        a = Fraction(factorial(self.n), pow(self.n, self.n))

        i = 0
        while(True):
            # 係数と降べき乗とで内積をとる
            prob = Fraction(0, 1)
            for x, y in zip(cofficients, variables):
                prob += x * y

            # 先頭の係数をかける
            prob *= a

            # 結果を出力
            print(n + i, format(float(prob), '.3g'))

            # 指定した確率まで計算できたらループを抜ける
            if prob > maxProb:
                break

            # 次のループの準備
            i += 1
            cofficients = self.getNextCofficients(i, cofficients)
            a *= Fraction(self.n + i, self.n)


n = int(sys.argv[1])
p = Fraction(int(sys.argv[2]), int(sys.argv[3]))

coupon = Coupon(n)
result = coupon.getProbabilityList(p)

実行結果

下記の例は $n = 6, p = 0.99$ である。
全6種のおまけの食玩を99%の確率で全部揃えたい場合は、最低でも36個買う必要があるということがわかる。

$ python3 Coupon.py 6 99 100
6 0.0154
7 0.054
8 0.114
9 0.189
10 0.272
11 0.356
12 0.438
13 0.514
14 0.583
15 0.644
16 0.698
17 0.745
18 0.785
19 0.819
20 0.848
21 0.873
22 0.893
23 0.911
24 0.925
25 0.938
26 0.948
27 0.957
28 0.964
29 0.97
30 0.975
31 0.979
32 0.982
33 0.985
34 0.988
35 0.99
36 0.992

式の導出

行列

全 $n$ 種のおまけがあるとする。このとき、次のような $n + 1$ 個の要素からなる縦ベクトル $a_r$ を考える。

$$a_r = \begin{pmatrix*}[c] a_{r, 0} \\ a_{r, 1} \\ \vdots \\ a_{r, n} \end{pmatrix*}$$

各 $a_{r, k}$ は $r$ 個買ったとき $k$ 個揃っている確率を表す。このとき、 $k \geq 1$ に対して次の式が成り立つ。

$$a_{r + 1, k} = \frac{n - k + 1}{n} a_{r, k - 1} + \frac{k}{n} a_{r, k}$$

式の意味としては、 $r + 1$ 個買って $k$ 種類揃っている確率は、 $r$ 個買ったとき、

• $k - 1$種類揃っているところから新たに $1$ 個買って $k$ 種になる場合
• $k$ 種類揃っているところから新たに $1$ 個買ってもそのまま $k$ 種類の場合

の確率を足したものになる。 $k = 0$ の時は、$1$ 個でも買うと、確率が $0$ となるのは明らかだろう。

上記の式を行列を使って表すと以下のようになる。

$$a_{r + 1} = \begin{pmatrix*}[c] \frac{0}{n} \\ \frac{n}{n} & \frac{1}{n} & & & \text{\huge{0}} \\ & \frac{n-1}{n} & \frac{2}{n} & \\ & & \frac{n-3}{n} & \frac{3}{n} & \\ & \text{\huge{0}} & & \ddots & \ddots \\ & & & & \frac{1}{n} & \frac{n}{n} \end{pmatrix*} a_r$$

その行列を $n$ 倍したものを $A$ とおく。また、 $n+1$ 次の正方行列 $Q$ を次のようにおく。

$$Q = \begin{pmatrix*}[c] {}_n\mathrm{C}_n \\ \vdots & \ddots & & \text{\huge{0}} \\ \vdots & & (-1)^{i+j} \,\, {}_{n-j}\mathrm{C}_{n-i} & & \\ \vdots & & & \ddots \\ (-1)^n \,\, {}_n\mathrm{C}_0 & \dots & \dots & \dots & {}_0\mathrm{C}_0 \end{pmatrix*}$$

この時、次の2つが成り立つ。

$$i) \,\,\,\, Q^{-1} = \begin{pmatrix*}[c] {}_n\mathrm{C}_n \\ \vdots & \ddots & & \text{\huge{0}} \\ \vdots & & {}_{n-j}\mathrm{C}_{n-i} & & \\ \vdots & & & \ddots \\ {}_n\mathrm{C}_0 & \dots & \dots & \dots & {}_0\mathrm{C}_0 \end{pmatrix*}$$

$$ii) \,\,\,\, AQ = Q \begin{pmatrix*}[c] 0 \\ & 1 & & & \text{\huge{0}} \\ & & 2 & \\ & & & 3 & \\ & \text{\huge{0}} & & & \ddots \\ & & & & & n \end{pmatrix*}$$

証明は、ただただ計算するだけなので省く。

一般項

前節の結果をまとめると、下記の式が得られる。

$$a_r = (\frac{1}{n}A)^r a_0 = \frac{1}{n^r} Q (Q^{-1}AQ)^r Q^{-1} a_0 = \frac{1}{n^r} Q \begin{pmatrix*}[c] 0^r \\ & 1^r & & & \text{\huge{0}} \\ & & 2^r & \\ & & & 3^r & \\ & \text{\huge{0}} & & & \ddots \\ & & & & & n^r \end{pmatrix*} Q^{-1} a_0$$

最右辺を求めることによって、最終的に以下の式が得られる。

$$a_r = \frac{1}{n^r} \begin{pmatrix*}[c] * \\ \vdots & \ddots & & \text{\huge{0}} \\ \vdots & & {}_{n-j}\mathrm{C}_{n-i} \,\, \sum_{k=j}^i (-1)^{i+k} \,\, {}_{i-j}\mathrm{C}_{i-k} \,\, k^r & & \\ \vdots & & & \ddots \\ * & \dots & \dots & \dots & * \end{pmatrix*} a_0$$

ここで、初項 $a_0$ を $a_{0, 0}$ だけ $1$ で他が $0$ である縦ベクトルとし、上記式で求められた縦ベクトル $a_r$ の最終行をみると、今回求めたかった $p_n(r)$ となる。つまり、上記の行列の $(i=)n$ 行 (j=)0 列の元になるので、

$$p_n(r) = \frac{1}{n^r} \sum_{k=0}^n (-1)^{n-k} \,\, {}_{n}\mathrm{C}_{k} \,\, k^r$$

となる(表記が綺麗になるように、うまいこと式変形をしていることに注意)。

系

$a_0$ を $a_{0, j}$ だけ $1$ で他が $0$ である縦ベクトルとしたら、最初 $j$ 個持っている状態となり、求められた $a_r$ の第 $i$ 行を見ると $i$ 個揃っている確率となっている。つまり、元の行列の $i$ 行 $j$ 列の元がその確率となる。

まとめると、全 $n$ 種類のおまけの食玩があり、最初 $a$ 個持っている。追加で $r$ 個買ったときちょうど $b$ 個揃っている確率は以下になる。

$$\frac{1}{n^r} \,\, {}_{n-a}\mathrm{C}_{n-b} \,\, \sum_{k=a}^b (-1)^{b+k} \,\, {}_{b-a}\mathrm{C}_{b-k} \,\, k^r$$

考察

今までの議論によって $p_n(r)$ を求めることができたが、その式を使って実際に値を求める際に少々不都合が起きる。まず、二項係数 $\mathrm{C}$ という、簡単に大きくなる値がある点である。それに加えて、小さい個数から順次求める際に、前の個数の計算式を使って次の計算ができないことである。

上記で求めた式をさらに変形して、計算がしやすい形にできないかどうかを見ていく。

式変形をしていく

式変形

二項定理から以下の式が成り立つ。

$$(e^x - 1)^n = \sum_{k=0}^n (-1)^{n-k} \,\, {}_n\mathrm{C}_k \,\, (e^x)^k$$

その両辺を $x$ で $r$ 回微分し、 $x \to 0$ に極限をとる。

$$\lim_{x \to 0} \frac{d^r}{dx^r} (e^x - 1)^n = \lim_{x \to 0} \sum_{k=0}^n (-1)^{n-k} \,\, {}_n\mathrm{C}_k \,\, k^r e^{kx} = \sum_{k=0}^n (-1)^{n-k} \,\, {}_n\mathrm{C}_k \,\, k^r = n^r p_n(r)$$

ここで、 $m = r - n$ とおくと、両辺の無限級数の係数を比較することで、

$$\frac{1}{r!} \lim_{x \to 0} \frac{d^r}{dx^r} (e^x - 1)^n = \frac{1}{m!} \lim_{x \to 0} \frac{d^m}{dx^m} \Big( \frac{e^x - 1}{x} \Big)^n$$

が成り立つことがわかる。ゆえに、

$$t_m(n) = \frac{1}{m!} \lim_{x \to 0} \frac{d^m}{dx^m} \Big( \frac{e^x - 1}{x} \Big)^n$$

と置いてやると、

$$p_n(r) = \frac{(n + m)!}{n^{n + m}} \frac{1}{m!} \lim_{x \to 0} \frac{d^m}{dx^m} \Big( \frac{e^x - 1}{x} \Big)^n = \frac{(n + m)!}{n^{n + m}} t_m(n)$$

と変形できる。この $t_m(n)$ を求めることに帰着できた。

$t_m(n)$ の意味としては、$\Big( \frac{e^x - 1}{x} \Big)^n$ をマクローリン展開して、$m$ 次の係数を導出しているに過ぎない。それがわかると複素積分でも求める方法もある。複素関数として、

$$f(z) = \Big( \frac{e^z - 1}{z} \Big) ^ n$$

と $f$ を定義してやると、

$$t_m(n) = \frac{1}{2 \pi i} \int_{C} \frac{f(z)}{z^{m + 1}} dz = \frac{1}{2 \pi i} \int_{C} \frac{(e^z - 1)^n}{z^{n + m + 1}} dz$$

と書くことができる。ただし、積分路 $C$ は原点を含む単純閉曲線とする。

実際に求めるにあたり、上の式は実質的に留数定理から導出したものになるので、留数定理を使わずに求める必要がある。複素積分の定義にしたがって求めようとしたが、途中で式が大変なことになり($\arctan$ が出てきたりする)、収拾がつかなくなった。

そのアプローチで成功すると、離散的な問題が複素積分で解けたということになり面白かったが現実はそう甘くなかった。

$t_m(n)$ を求めていく

最初の方で出てきた漸化式に $k=n$ を代入してみると、

$$a_{r + 1, n} = \frac{1}{n} a_{r, n - 1} + a_{r, n}$$

となる。ここで $a_{r, n - 1}$ というのは、 $r$ 個買った時にちょうど $n-1$ 個揃っている確率になるので、系の節で触れた式を使うと、

$$\begin{aligned} a_{r, n - 1} &= \frac{1}{n^{n + m}} \,\, {}_{n}\mathrm{C}_{1} \,\, \sum_{k=0}^{n-1} (-1)^{n-1+k} \,\, {}_{n - 1}\mathrm{C}_{n - 1 -k} \,\, k^r \\ &= \frac{1}{n^{n + m}} \,\, n \,\, \sum_{k=0}^{n-1} (-1)^{n-1+k} \,\, {}_{n - 1}\mathrm{C}_{k} \,\, k^r \\ &= \frac{1}{n^{n + m - 1}} (n - 1)^r p_{n-1}(r) \\ &= \frac{(n-1)^r}{n^{n + m - 1}} \frac{(n - 1 + m + 1)!}{(n - 1)^{n - 1 + m + 1}} t_{m + 1}(n - 1) \\ &= \frac{(n+m)!}{n^{n + m - 1}} t_{m + 1}(n - 1) \end{aligned}$$

となるので、元の漸化式に代入すると、

$$\frac{(n + m + 1)!}{n^{n + m + 1}} t_{m + 1}(n) = \frac{(n + m)!}{n^{n + m}} t_{m + 1}(n - 1) + \frac{(n + m)!}{n^{n + m}} t_m(n)$$

となることがわかる。両辺 $n^{n+m}/(n-1)!$ 倍してやって移行し、

$$\frac{(n + m + 1)!}{n!} t_{m + 1}(n) - \frac{(n + m)!}{(n - 1)!} t_{m + 1}(n - 1) = \frac{(n + m)!}{(n - 1)!} t_m(n)$$

$n$ を $1$ から $n$ まで動かしたときの和をとってみると、

$$\sum_{k=1}^n \Big(\frac{(k + m + 1)!}{k!} t_{m + 1}(k) - \frac{(k + m)!}{(k - 1)!} t_{m + 1}(k - 1) \Big) = \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} t_m(k)$$

左辺がうまいこと打ち消しあって、以下の式が得られる。

$$\frac{(n + m + 1)!}{n!} t_{m + 1}(n) - (1 + m)! \,\, t_{m + 1}(0) = \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} t_m(k)$$

ここで、

$$t_m(0) = \frac{1}{m!} \lim_{x \to 0} \frac{d^m}{dx^m} \Big( \frac{e^x - 1}{x} \Big)^0 = \frac{1}{m!} \lim_{x \to 0} \frac{d^m}{dx^m} 1 = 0$$

$$t_0(n) = \frac{1}{0!} \lim_{x \to 0} \frac{d^0}{dx^0} \Big( \frac{e^x - 1}{x} \Big)^n = \lim_{x \to 0} \Big( \frac{e^x - 1}{x} \Big)^n = 1^n = 1$$

となっているので、

$$t_{m + 1}(n) = \frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} t_m(k), \,\,\,\, t_0(n) = 1$$

という関係式が得られる。

和を求めていくための準備

和を求めやすくするために、降べき乗 $n^{\underline{k}} \,\, (k \geq 0)$ を導入する。これは次で定義されるものである。

$$n^{\underline{k}} = \begin{cases} 1 & k = 0 \\ n * (n - 1)^{\underline{k-1}} & \mathrm{o.w.} \end{cases}$$

例として $n^{\underline{3}} = n(n-1)(n-2)$ となる。また、その例のときに $n = 2$ とすると $0$ になるように、 $k \gt n$ であれば $0$ となる。

まず手始めに、 $t_m(n) = n^{\underline{i}} \,\, (i \geq 0)$ として求めていく。 $i \geq 1$ であれば、

$$\begin{aligned} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} k^{\underline{i}} &= \sum_{k=1}^n k(k+1) \cdots (k+m) * k (k - 1) \cdots (k - i + 1) \\ &= \sum_{k=1}^n k * (k - i + 1) \cdots (k-1)k(k+1) \cdots (k+m) \\ &= \sum_{k=1}^n (k - i) * (k - i + 1) \cdots (k-1)k(k+1)\cdots (k+m) \\ & \,\,\,\,\,\,\,\, + i * \sum_{k=1}^n (k - i + 1) \cdots (k-1)k(k+1) \cdots (k+m) \end{aligned}$$

となる。ここで $c, d \geq 0$ としたとき、

$$\begin{aligned} &\sum_{k=1}^n (k-c) \cdots k \cdots (k+d) \\ &= \frac{1}{c+d+2} \\ & \,\,\,\,\,\,\,\, * \sum_{k=1}^n \Big( (k-c) \cdots k \cdots (k+d)(k+d+1) - (k-c-1)(k-c) \cdots k \cdots (k+d) \Big) \\ &= \frac{1}{c+d+2} (n-c) \cdots (n+d)(n+d+1) \end{aligned}$$

であるので、

$$\begin{aligned} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} k^{\underline{i}} &= \frac{1}{m + i + 2} (n - i) \cdots (n + m + 1) \\ & \,\,\,\,\,\,\,\, + i * \frac{1}{m + i + 1} (n - i + 1) \cdots (n + m + 1) \\ &= \Big( \frac{n - i}{m + i + 2} + \frac{i}{m + i + 1} \Big) (n - i + 1) \cdots (n + m + 1) \\ \end{aligned}$$

が求められる。両辺 $n!/(n + m + 1)!$ をかけてやると、

$$\begin{aligned} &\frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} k^{\underline{i}} \\ &= \Big( \frac{n - i}{m + i + 2} + \frac{i}{m + i + 1} \Big) \frac{n! (n - i + 1) \cdots (n + m + 1)}{(n + m + 1)!} \\ &= \Big( \frac{n - i}{m + i + 2} + \frac{i}{m + i + 1} \Big) \frac{n!}{(n - i)!} \\ &= \frac{1}{m + i + 2} \frac{n!}{(n - i - 1)!} + \frac{i}{m + i + 1} \frac{n!}{(n - i)!} \\ &= \frac{1}{m + i + 2} n^{\underline{i + 1}} + \frac{i}{m + i + 1} n^{\underline{i}} \end{aligned}$$

となることがわかる。一方 $i = 0$ の時は、

$$\begin{aligned} &\frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} \\ &= \frac{n!}{(n + m + 1)!} \sum_{k=1}^n k(k+1)\cdots (k+m) \\ &= \frac{n!}{(n + m + 1)!} \frac{1}{m + 2} n(n + 1) \cdots (n + m + 1) \\ &= \frac{1}{m + 2} n^{\underline{1}} \end{aligned}$$

となり、 $i \geq 1$ として求めた式に $i = 0$ を代入したものになるので、$i \geq 1$ という条件は $i \geq 0$ として問題ない。

和を求める

前節の結果をみると、降べき乗 $n^{\underline{i}}$ に対して和を求めると $n^{\underline{i + 1}}$ と $n^{\underline{i}}$ の項が出てきた。つまり、次数が1つ上がったことになる。 $t_0(n) = 1$ であることから、 $t_m(n)$ は最大次数が $m$ 次の降べき乗の多項式となる。そこで $c_{m, i} \,\, (0 \leq i \leq m)$ を用意すると、

$$t_m(n) = \sum_{i=0}^m c_{m, i} \,\, n^{\underline{i}}$$

と書くことができる(このとき $c_{0, 0} = 1$ である)ので、

$$\begin{aligned} t_{m + 1}(n) &= \frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} t_m(k) \\ &= \frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} \sum_{i=0}^m c_{m, i} \,\, k^{\underline{i}} \\ &= \sum_{i=0}^m c_{m, i} \,\, \frac{n!}{(n + m + 1)!} \sum_{k=1}^n \frac{(k + m)!}{(k - 1)!} k^{\underline{i}} \\ &= \sum_{i=0}^m c_{m, i} \,\, \Big( \frac{1}{m + i + 2} n^{\underline{i + 1}} + \frac{i}{m + i + 1} n^{\underline{i}} \Big) \\ &= c_{m, m} \,\, \frac{1}{m + m + 2} n^{\underline{m + 1}} \\ & \,\,\,\,\,\,\,\, + \sum_{i=1}^{m} \Big( \frac{c_{m, i - 1}}{m + i + 1} n^{\underline{i}} + \frac{i * c_{m, i}}{m + i + 1} n^{\underline{i}} \Big) \\ & \,\,\,\,\,\,\,\, + c_{m, 0} \,\, \frac{0}{m + 1} n^{\underline{0}} \\ &= c_{m, m} \,\, \frac{1}{2(m + 1)} n^{\underline{m + 1}} \\ & \,\,\,\,\,\,\,\, + \sum_{i=1}^{m} \frac{c_{m, i - 1} + i * c_{m, i}}{m + i + 1} n^{\underline{i}} \\ & \,\,\,\,\,\,\,\, + 0 \end{aligned}$$

となることがわかる。 $t_{m + 1}(n)$ と係数比較することによって、

$$c_{m + 1, m + 1} = \frac{1}{2(m + 1)} c_{m, m}$$

$$c_{m + 1, i} = \frac{1}{m + i + 1} (c_{m, i - 1} + i * c_{m, i}) \,\,\,\, (1 \leq i \leq m)$$

$$c_{m + 1, 0} = 0$$

という漸化式が得られる。

$\mathrm{T}$ の定義

前節に出てきた $c_{m, i}$ を ${}_m\mathrm{T}_i$ という表記にしてみる。 $0 \leq i \leq m$ に対して ${}_m\mathrm{T}_k$ を、

$${}_0\mathrm{T}_0 = 1, \,\,\,\, {}_m\mathrm{T}_0 = 0 \,\,\,\, (m \geq 1), \,\,\,\, {}_m\mathrm{T}_m = \frac{1}{2^m m!} \,\,\,\, (m \geq 0)$$

$${}_m\mathrm{T}_k = \frac{1}{m + k} ({}_{m - 1}\mathrm{T}_{k - 1} + k \,\, {}_{m - 1}\mathrm{T}_k) \,\,\,\, (m, k \geq 1)$$

定義する。二項係数 ${}_n\mathrm{C}_k$ みたいな表記にしたのは、パスカルの三角形と同じような三角形を作ることができるからである。

$$\begin{array}{ccccccccccccc} & & & & & & 1 & & & & & & \\ & & & & & 0 & & \frac{1}{2} & & & & & \\ & & & & 0 & & \frac{1}{6} & & \frac{1}{8} & & & & \\ & & & 0 & & \frac{1}{24} & & \frac{1}{12} & & \frac{1}{48} & & & \\ & & 0 & & \frac{1}{120} & & \frac{5}{144} & & \frac{1}{48} & & \frac{1}{384} & & \\ & 0 & & \frac{1}{720} & & \frac{1}{90} & & \frac{7}{576} & & \frac{1}{288} & & \frac{1}{3840} & \\ 0 & & \frac{1}{5040} & & \frac{17}{5760} & & \frac{137}{25920} & & \frac{1}{384} & & \frac{1}{2304} & & \frac{1}{46080} \\ \end{array}$$

一例を見てみると、

$${}_4\mathrm{T}_2 = \frac{1}{4 + 2}({}_3\mathrm{T}_1 + 2 \,\, {}_3\mathrm{T}_2) = \frac{1}{6}(\frac{1}{24} + 2 \frac{1}{12}) = \frac{5}{144}$$

のように上段の左右2項を使って導出することができる。つまり、上段があれば次の段の計算ができることになる。

最終的な導出結果

今まで計算してきた結果をまとめると、下記の式として表現することができる。

$$p_n(r) = \frac{(n + m)!}{n^{(n + m)}} \sum^m_{k=0} {}_m\mathrm{T}_k \,\, n^{\underline{k}}$$

プログラムを書くにあたり変形

降べき乗の性質から $k \gt n$ のとき $n^{\underline{k}} = 0$ となるので、 $p_n(r)$ を求める際には最大 $n$ 項まで求めればいいことになる。また、 $\mathrm{T}$ が定義できない場合は $0$ としてやることで、常に $n$ 項まででよくなる。

以上をまとめると、

$$p_n(r) = \frac{(n + m)!}{n^{(n + m)}} \sum^n_{k=0} {}_m\mathrm{T}_k \,\, n^{\underline{k}}$$

という式に変形できる。これにより、冒頭で述べた式が導出できた。

正論は時に人を傷つける

ここである人は言った。

プログラムで求めるのであれば、最初の漸化式で求めればよくね？

$$a_{r + 1, k} = \frac{n - k + 1}{n} a_{r, k - 1} + \frac{k}{n} a_{r, k}$$

上記の式が成り立っているので、下記のようにコードを書いたら、

SeironIsNotGoodForPeople.py

from fractions import Fraction
import sys

class SeironIsNotGoodForPeople():
    def __init__(self, n: int):
        if n <= 0:
            raise Exception("n is wrong")

        self.n = n

    def getNextProbs(self, probs: list[Fraction]) -> list[Fraction]:
        return list(map(
            lambda same, diff, index: ((self.n - Fraction(index, 1) + Fraction(1, 1)) * diff + Fraction(index, 1) * same) / Fraction(self.n, 1),
            probs,
            [Fraction(0, 1)] + probs[:-1],
            range(len(probs)),
        ))

    def getProbabilityList(self, maxProb: Fraction) -> None:
        if maxProb >= Fraction(1, 1):
            raise Exception("p is wrong")

        # 最初の確率 [1, 0, 0, ..., 0] を作成
        probs = [Fraction(0, 1) for _ in range(self.n + 1)]
        probs[0] = Fraction(1, 1)

        i = 0
        while(True):
            # 結果を出力
            print(i, list(map(lambda x: format(float(x), '.3g'), probs)))

            # 指定した確率まで計算できたらループを抜ける
            if probs[-1] > maxProb:
                break

            # 次のループの準備
            i += 1
            probs = self.getNextProbs(probs)


n = int(sys.argv[1])
p = Fraction(int(sys.argv[2]), int(sys.argv[3]))

coupon = SeironIsNotGoodForPeople(n)
result = coupon.getProbabilityList(p)

コード量も少なく、i ループ目の probs[k] には i 個買ったときにちょうど k 個揃っている確率が入っているので、

$ python3 SeironIsNotGoodForPeople.py 6 99 100
0 ['1', '0', '0', '0', '0', '0', '0']
1 ['0', '1', '0', '0', '0', '0', '0']
2 ['0', '0.167', '0.833', '0', '0', '0', '0']
3 ['0', '0.0278', '0.417', '0.556', '0', '0', '0']
4 ['0', '0.00463', '0.162', '0.556', '0.278', '0', '0']
5 ['0', '0.000772', '0.0579', '0.386', '0.463', '0.0926', '0']
6 ['0', '0.000129', '0.0199', '0.231', '0.502', '0.231', '0.0154']
7 ['0', '2.14e-05', '0.00675', '0.129', '0.45', '0.36', '0.054']
8 ['0', '3.57e-06', '0.00227', '0.069', '0.365', '0.45', '0.114']
9 ['0', '5.95e-07', '0.000759', '0.036', '0.278', '0.497', '0.189']
10 ['0', '9.92e-08', '0.000254', '0.0185', '0.203', '0.506', '0.272']
11 ['0', '1.65e-08', '8.46e-05', '0.00943', '0.145', '0.49', '0.356']
12 ['0', '2.76e-09', '2.82e-05', '0.00477', '0.101', '0.456', '0.438']
13 ['0', '4.59e-10', '9.41e-06', '0.0024', '0.0698', '0.414', '0.514']
14 ['0', '7.66e-11', '3.14e-06', '0.00121', '0.0477', '0.368', '0.583']
15 ['0', '1.28e-11', '1.05e-06', '0.000606', '0.0324', '0.323', '0.644']
16 ['0', '2.13e-12', '3.48e-07', '0.000304', '0.0219', '0.28', '0.698']
17 ['0', '3.54e-13', '1.16e-07', '0.000152', '0.0148', '0.24', '0.745']
18 ['0', '5.91e-14', '3.87e-08', '7.61e-05', '0.00992', '0.205', '0.785']
19 ['0', '9.85e-15', '1.29e-08', '3.81e-05', '0.00665', '0.174', '0.819']
20 ['0', '1.64e-15', '4.3e-09', '1.91e-05', '0.00445', '0.148', '0.848']
21 ['0', '2.74e-16', '1.43e-09', '9.53e-06', '0.00298', '0.124', '0.873']
22 ['0', '4.56e-17', '4.78e-10', '4.77e-06', '0.00199', '0.105', '0.893']
23 ['0', '7.6e-18', '1.59e-10', '2.38e-06', '0.00133', '0.0879', '0.911']
24 ['0', '1.27e-18', '5.31e-11', '1.19e-06', '0.000887', '0.0737', '0.925']
25 ['0', '2.11e-19', '1.77e-11', '5.96e-07', '0.000592', '0.0617', '0.938']
26 ['0', '3.52e-20', '5.9e-12', '2.98e-07', '0.000395', '0.0516', '0.948']
27 ['0', '5.86e-21', '1.97e-12', '1.49e-07', '0.000264', '0.0431', '0.957']
28 ['0', '9.77e-22', '6.56e-13', '7.45e-08', '0.000176', '0.036', '0.964']
29 ['0', '1.63e-22', '2.19e-13', '3.73e-08', '0.000117', '0.0301', '0.97']
30 ['0', '2.71e-23', '7.29e-14', '1.86e-08', '7.82e-05', '0.0251', '0.975']
31 ['0', '4.52e-24', '2.43e-14', '9.31e-09', '5.21e-05', '0.021', '0.979']
32 ['0', '7.54e-25', '8.09e-15', '4.66e-09', '3.48e-05', '0.0175', '0.982']
33 ['0', '1.26e-25', '2.7e-15', '2.33e-09', '2.32e-05', '0.0146', '0.985']
34 ['0', '2.09e-26', '8.99e-16', '1.16e-09', '1.54e-05', '0.0122', '0.988']
35 ['0', '3.49e-27', '3e-16', '5.82e-10', '1.03e-05', '0.0101', '0.99']
36 ['0', '5.82e-28', '9.99e-17', '2.91e-10', '6.87e-06', '0.00845', '0.992']

最初のコードより簡潔で情報量増えていいんじゃないの？

ここで声を大にして言いたいことがあります。
.
.
.
.
.
.
.
.
.
.
.
.
「正論は時に人を傷つける」

数学的な意味

正直なところ、確率を求めるだけなら最初の漸化式で良い。ただ、食玩問題の性質を探っていくには数式を見て行かないとわからない。

どういうことかというと、最初に出した漸化式、

$$a_{r + 1, k} = \frac{n - k + 1}{n} a_{r, k - 1} + \frac{k}{n} a_{r, k}$$

は $n$ と $r$ に依存していて、項の前後の関係しかわからない。そこから、行列を用いて一般項、

$$p_n(r) = \frac{1}{n^r} \sum_{k=0}^n (-1)^{n-k} \,\, {}_{n}\mathrm{C}_{k} \,\, k^r$$

を求めたことにより、 $r$ 個買った時の関係の見通しが良くなった。買う個数が増えるにあたり、累乗部分が大きくなっていくことが読み取れる。そこからさらに式変形をして、

$$p_n(r) = \frac{(n + m)!}{n^{(n + m)}} \sum^m_{k=0} {}_m\mathrm{T}_k \,\, n^{\underline{k}}$$

を得た。その式は、本質的な部分は ${}_m\mathrm{T}_k$ に帰着させることになった。その $m$ は $r-n$ で求められるものではあるが、必要最低限から追加で買った個数という意味としてみると $n$ には依存しない。つまり ${}_m\mathrm{T}_k$ は $n$ に依存しないものとなった。

$n$ に依存しているかどうかは、上で紹介した2つのプログラムで、前者の getNextCofficients() は staticmethod になっており、後者の getNextProbs() はそうなっていないところに表れている。

${}_m\mathrm{T}_k$ は必要最低限から追加で買った個数 $m$ だけに依存し、おまけの数 $n$ には依存しない。それにより、より食玩問題の本質を表していると思っている。ただ、この部分はまだまだ研究中のところがあるので、よく分かっていない。

研究の一例として、 ${}_m\mathrm{T}_k$ の性質を $k$ を固定して $m$ の式で表していくということをしている。

$$\begin{aligned} {}_m\mathrm{T}_0 &= \begin{cases} 1 & m = 0 \\ 0 & \mathrm{o.w.} \end{cases} \\ {}_m\mathrm{T}_1 &= \frac{1}{(m + 1)!} \\ {}_m\mathrm{T}_2 &= \frac{1}{(m + 2)!}(2^{m + 1} - (m + 3)) \\ {}_m\mathrm{T}_3 &= \frac{1}{(m + 3)!}(\frac{1}{2} 3^{m + 2} - (m + 5) 2^{m + 1} + \frac{1}{2}(m^2 + 7m + 13)) \\ \vdots \\ {}_m\mathrm{T}_{m - 1} &= \frac{1}{3 * 2^{m - 1} (m - 2)!} \\ {}_m\mathrm{T}_{m} &= \frac{1}{2^m m!} \end{aligned}$$

上の式を眺めてみると、途中で出てきた式、

$$p_n(r) = \frac{1}{n^r} \sum_{k=0}^n (-1)^{n-k} \,\, {}_{n}\mathrm{C}_{k} \,\, k^r$$

に似たような部分が出てきていそうではある(その式から導出しているので、当たり前っちゃ当たり前かもだが)。

おわりに

食玩問題については15年ほど前から研究している。というのも、当時AKB48がランダムでメンバーの写真が入っているグッズを出し、全部揃えるのは不可能じゃないかとネット上で話題になったことがあった。感覚的に不可能に近いことは分かるが、どのくらい不可能なのか？どのくらい買えばいいのか？という点が気になった。それが研究を始めたきっかけである(ちなみに僕はAKB48のファンでも何でもないです)。

その問題の答えとしては(全48種のおまけでは)、全部揃っている確率を50%を越すには204個、99%を越すには403個買う必要があることが、上記プログラムを動かすことでわかる。

この食玩問題の研究をどこかでまとめたいなと思い、プログラムを書くことによって技術的な話だろ！とすることで、 Zenn に投稿にすることにした。書き終わって改めてみると「完全に数学の話じゃん！」となったことは内緒。