coupon_collectors

$N$ 種類の食玩があと $k$ 種類で揃うというときの回数の期待値 $E_k$ の漸化式を求め、$N$ 種類すべて揃うときの回数の期待値 $E_0$ を求める。


期待値 $E_k$ は最初の１種類目の回数の期待値 $1$ および、ひとつ購入して、ひとつ揃う確率の期待値 $E_{k-1}$ と揃わない確率の期待値 $E_k$ を足したものなので、以下のようになる。
\[
\begin{aligned}
E_k &= 1 + \frac{k}{N}E_{k-1} + \frac{N-k}{N}E_k\\
\end{aligned}
\]
これを変形して一般式 $E_k$ および $E_0$ を求めると以下のようになる。
\[
\begin{aligned}
E_k - \frac{N-k}{N}E_k &= \frac{k}{N}E_{k-1} + E_N\\
E_k\left(\frac{k}{N}\right) &= \frac{k}{N}E_{k-1} + 1\\
E_k &= E_{k-1} + \frac{N}{k}\\
E_k &= N\sum_{i=k}^{N}i^{-1}\\
E_0 &= N\sum_{i=0}^{N}i^{-1}\\
&= N\left(1+\frac12+\cdots+\frac1{N-1}+\frac1N\right)\\
\end{aligned}
\]


例えば、$N=8$ 種類だとおよそ $E_8=22$ なので、（あくまで）平均して $22$ 個買うと揃うことになる。


以下は、参考文献の「食玩問題」の<a href="http://aquarius10.cse.kyutech.ac.jp/~otabe/shokugan/p.pdf">導出過程</a>に示されているものであるが、最後の式と最初からの式の添字が異なっていたりするので、ここでは、多少整理した形で記録しておく。


$N$ 種類の食玩を $n$ 個買ったとして、$k$ 種類揃った確率を $p_k(n)$ と表すと、
<ol>
<li>持っている $1$ 個と同じものに当たってしまう確率 $=p_1(n+1)$</li>
<li>持っている $2$ 個と同じものに当たってしまう確率 $+$ 持っていない他の $N-1$ 個のものに当たる確率 $=p_2(n+1)$</li>
<li>持っている $3$ 個と同じものに当たってしまう確率 $+$ 持っていない他の $N-2$ 個のものに当たる確率 $=p_3(n+1)$</li>
<li>持っている $k$ 個と同じものに当たってしまう確率 $+$ 持っていない他の $N-(k-1)$ 個のものに当たる確率 $=p_k(n+1)$</li>
<li>持っている $N$ 個と同じものに当たってしまう確率 $+$ 持っていない他の $N-(N-1)$ 個のものに当たる確率 $=p_N(n+1)$</li>
</ol>
これを漸化式で表すと以下のようになり、
\[
\begin{aligned}
p_1(n+1) &= \frac{1}{N}p_1(n)\\
p_2(n+1) &= \frac{2}{N}p_2(n) + \frac{N-1}{N}p_1(n)\\
p_3(n+1) &= \frac{3}{N}p_3(n) + \frac{N-2}{N}p_2(n)\\
&\vdots\\
p_{k}(n+1) &= \frac{k}{N}p_{k}(n) + \frac{N-(k-1)}{N}p_{k-1}(n)\\
&\vdots\\
p_N(n+1) &= \frac{N}{N}p_N(n) + \frac{1}{N}p_{N-1}(n)\\
\end{aligned}
\]
これを参考文献「食玩問題」のように行列で表して対角化を施したりすると、以下のような一般式が得られることが示されている。
\[
p_m(n) = \sum_{i=1}^{m}(-1)^{m-i}\dbinom{N-i}{m-i}\dbinom{N-1}{i-1}\left(\frac{i}{N}\right)^{n-1}
\]

ここで、
\[
\begin{aligned}
\dbinom{N-i}{m-i}\dbinom{N-1}{i-1}
&= \frac{(N-i)!}{(m-i)!(N-m)!}\frac{(N-1)!}{(i-1)!(N-i)!}\\
&= \frac{(N-1)!}{(m-1)!(N-m)!}\frac{(m-1)!}{(m-i)!(i-1)!}
= \dbinom{N-1}{m-1}\dbinom{m-1}{m-i}
\end{aligned}
\]
ゆえに、

\[
p_m(n) = \dbinom{N-1}{m-1}\sum_{i=1}^{m}(-1)^{m-i}\dbinom{m-1}{m-i}\left(\frac{i}{N}\right)^{n-1}
\]
これは参考文献の最後に示されている式と添字が異なるので、合うように逆順にする。

$i\to N-i$ とすると以下のようになる。
\[
\begin{aligned}
p_m(n)
&= \sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{i}{m-N+i}\dbinom{N-1}{N-i-1}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{i}{N-m}\dbinom{N-1}{i}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{i}{N-m}\dbinom{N-1}{N-1-i}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{i}{N-m}\dbinom{N}{N-i}\left(\frac{N-i}{N}\right)^n\\
&= \sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{i}{N-m}\dbinom{N}{i}\left(\frac{N-i}{N}\right)^n\\
\end{aligned}
\]

ここで、
\[
\begin{aligned}
\dbinom{i}{N-m}\dbinom{N-1}{i}
&= \frac{i!}{(N-m)!(i-N+m)!}\frac{(N-1)!}{i!(N-1-i)!}\\
&= \frac{(N-1)!}{(m-1)!(N-m)!}\frac{(m-1)!}{(m-N+i)!(N-1-i)!}
= \dbinom{N-1}{m-1}\dbinom{m-1}{m-N+i}
\end{aligned}
\]
ゆえに、

\[
\begin{aligned}
p_m(n)
&= \dbinom{N-1}{m-1}\sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{m-1}{m-N+i}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \dbinom{N-1}{m-1}\sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{m-1}{N-i-1}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \frac{m}{N}\dbinom{N}{m}\sum_{i=N-m}^{N-1}(-1)^{m-N+i}\frac{N-i}{m}\dbinom{m}{N-i}\left(\frac{N-i}{N}\right)^{n-1}\\
&= \dbinom{N}{m}\sum_{i=N-m}^{N-1}(-1)^{m-N+i}\dbinom{m}{N-i}\left(\frac{N-i}{N}\right)^n\\
\end{aligned}
\]


これらは、参考文献の最後の式にて $a=0, b=N$ としたもの等価である。
\[
p_{a=0,b,m}(n) = \sum_{i=b-m}^{b}(-1)^{b-m+i}\dbinom{i}{b-m}\dbinom{b}{i}\left(\frac{N-i}{N}\right)^n
\]

ここで、二項係数の２つ積の１つ目の添字の差が一方の２つ目の添字に等しい、もしくは、二項係数の２つ積の互い違いの添字が等しいとき、以下が成り立つので、
\[
\begin{aligned}
\dbinom{n}{h}\dbinom{n-h}{k} = \dbinom{n}{k}\dbinom{n-k}{h}
\Leftrightarrow&
\dbinom{n}{n-h}\dbinom{n-h}{k} = \dbinom{n}{n-k}\dbinom{n-k}{h}
\\
\end{aligned}
\]
\[
\begin{aligned}
\dbinom{n}{h}\dbinom{n-h}{k}
&= \frac{n!}{h!(n-h)!}\frac{(n-h)!}{k!(n-h-k)!}\\
&= \frac{n!}{k!(n-k)!}\frac{(n-k)!}{h!(n-k-h)!}
= \dbinom{n}{k}\dbinom{n-k}{h}
\end{aligned}
\]
ゆえに、

\[
p_{a=0,b,m}(n) = \dbinom{b}{m}\sum_{i=b-m}^{b}(-1)^{b-m+i}\dbinom{m}{m-b+i}\left(\frac{N-i}{N}\right)^n
\]
特に $m=b$ であれば：
\[
p_{a=0,b,m=b}(n) = \sum_{i=0}^{b}(-1)^{i}\dbinom{b}{i}\left(\frac{N-i}{N}\right)^n
\]
念の為、参考文献の最後の式をここでの表し方で記しておく。すでに $a$ 種類揃っている $N$ 種類の食玩のうち $b$ 種類揃える。さらに $n$ 個買ったとして、$m$ 種類揃う確率は以下の通りである。
\[
p_{a,b,m}(n) = \sum_{i=b-m}^{b-a}(-1)^{b-m+i}\dbinom{i}{b-m}\dbinom{b-a}{i}\left(\frac{N-i}{N}\right)^n
\]
\[
p_{a,b,m}(n) = \dbinom{b-a}{m-a}\sum_{i=b-m}^{b-a}(-1)^{b-m+i}\dbinom{m-a}{m-b+i}\left(\frac{N-i}{N}\right)^n
\]
特に $m=b$ であれば：
\[
p_{a,b,m=b}(n) = \sum_{i=0}^{b-a}(-1)^i\dbinom{b-a}{i}\left(\frac{N-i}{N}\right)^n
\]


以下にいくつかの数値解の計算結果を示す。
<form name="ccp" style="white-space: nowrap" onInput='
if (parseInt(parameter_N.value) > 0 &&
 parseInt(parameter_n.value) >= 0 &&
 parseInt(parameter_N.value) >= parseInt(parameter_b.value) &&
 parseInt(parameter_b.value) >= parseInt(parameter_m.value) &&
 parseInt(parameter_m.value) >= parseInt(parameter_a.value) &&
 parseInt(parameter_a.value) >= 0) {
 parameters_condition.removeAttribute("style");
}
else {
 parameters_condition.style = "color: red;";
};
'>
$N$: <input type="text" name="parameter_N" value="8"/> 
$a$: <input type="text" name="parameter_a" value="0"/>
$b$: <input type="text" name="parameter_b" value="8"/> 
$m$: <input type="text" name="parameter_m" value="8"/>
$n$: <input type="text" name="parameter_n" value="22"/> 
<output class="small" name="parameters_condition">$N\gt 0, n\ge 0, N\ge b\ge m\ge a\ge 0$</output> 
<input type="button" value="計算" class="mt-2 mb-2" onClick='
function binomial(n, k) {
 let r = 1;
 for (let i=1; i<=k; ++i)
 r *= (n+1-i)/i;
 return r;
}
/*
for (let n=0; n<=10; ++n) {
 let s = ""
 for (let k=0; k<=n; ++k) {
 s += `${binomial(n, k)} `
 }
 console.log(s)
}
*/
function ccp0d(N, m, n) {
 let p = Array(N+1);
 for (let k=0; k<=N; ++k) {
 p[k] = 0;
 }
 p[1] = 1;
 for (let j=1; j<n; ++j) {
 for (let k=N; k>=1; --k) {
 p[k] = k/N*p[k] + (N-(k-1))/N*p[k-1];
 }
 }
 return p[m];
}
function ccp0(N, m, n) {
 let r = 0;
 if (n < 1) return r;
 for (let i=1; i<=m; ++i)
 //r += Math.pow(-1, m-i)*binomial(N-i, m-i)*binomial(N-1, i-1)*Math.pow(i/N, n-1);
 r += binomial(N-1, m-1)*Math.pow(-1, m-i)*binomial(m-1, m-i)*Math.pow(i/N, n-1);
 return r;
}
function ccp1(N, m, n) {
 let r = 0;
 if (n < 1) return r;
 for (let i=N-m; i<=N-1; ++i)
 //r += Math.pow(-1, m-N+i)*binomial(i, N-m)*binomial(N-1, i)*Math.pow((N-i)/N, n-1);
 r += binomial(N-1, m-1)*Math.pow(-1, m-N+i)*binomial(m-1, N-1-i)*Math.pow((N-i)/N, n-1);
 return r;
}
function _ccp1(N, m, n) {
 let r = 0;
 if (n < 1) return r;
 for (let i=N-m+1; i<=N; ++i)
 //r += Math.pow(-1, m-N+i-1)*binomial(i-1, N-m)*binomial(N-1, i-1)*Math.pow((N-i+1)/N, n-1);
 r += binomial(N-1, m-1)*Math.pow(-1, m-N+i-1)*binomial(m-1, N-i)*Math.pow((N-i+1)/N, n-1);
 return r;
}
function ccpd(N, a, b, m, n) {
 let p = Array(N+1);
 for (let k=0; k<=N; ++k) {
 p[k] = 0;
 }
 p[a] = 1;
 for (let j=0; j<n; ++j) {
 for (let k=b; k>=a; --k) {
 p[k] = (N-b+k)/N*p[k] + (a<k ? (b-k+1)/N*p[k-1] : 0);
 }
 }
 return p[m];
}
function ccp(N, a, b, m, n) {
 let r = 0;
 for (let i=b-m; i<=b-a; ++i)
 //r += Math.pow(-1, m-b+i)*binomial(i, b-m)*binomial(b-a, i)*Math.pow((N-i)/N, n);
 r += binomial(b-a, m-a)*Math.pow(-1, (m-b)+i)*binomial(m-a, m-b+i)*Math.pow((N-i)/N, n);
 return r;
}
output_p_m[0].value = `${ccp0(parseInt(parameter_N.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
output_p_m[1].value = `${ccp1(parseInt(parameter_N.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
output_p_m[2].value = `${ccp0d(parseInt(parameter_N.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
output_p_a_b_m[0].value = `${ccp(parseInt(parameter_N.value), 0, parseInt(parameter_N.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
output_p_a_b_m[1].value = `${ccp(parseInt(parameter_N.value), parseInt(parameter_a.value), parseInt(parameter_b.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
output_p_a_b_m[2].value = `${ccpd(parseInt(parameter_N.value), parseInt(parameter_a.value), parseInt(parameter_b.value), parseInt(parameter_m.value), parseInt(parameter_n.value))}`;
'/> 
$p_m(n)$: <output name="output_p_m" style="display: inline-block; width: 12rem;">0.6242501484436673</output>
… 参考文献の最初の一般式 
$p_m(n)$: <output name="output_p_m" style="display: inline-block; width: 12rem;">0.6242501484436674</output>
… その添字を逆順にした式 
$p_m(n)$: <output name="output_p_m" style="display: inline-block; width: 12rem;" style="display: inline-block; width: 12rem;">0.6242501484436673</output>
… 冒頭の漸化式 $p_1(0)=1$ 
$p_{a,b,m}(n)$: <output name="output_p_a_b_m" style="display: inline-block; width: 12rem;">0.6242501484436674</output>
… 次式で $a=0, b=N$ とした式 
$p_{a,b,m}(n)$: <output name="output_p_a_b_m" style="display: inline-block; width: 12rem;">0.6242501484436674</output>
… 参考文献の最後の一般式 
$p_{a,b,m}(n)$: <output name="output_p_a_b_m" style="display: inline-block; width: 12rem;">0.6242501484436674</output>
… 末尾の漸化式 $p_a(0)=1$ 
</form>


ちなみに、元々の漸化式は以下であり、
\[
\begin{aligned}
p_a(n+1) &= \frac{N-b+a+0}{N}p_a(n)\\
p_{a+1}(n+1) &= \frac{N-b+a+1}{N}p_{a+1}(n) + \frac{b-a-0}{N}p_{a}(n)\\
p_{a+2}(n+1) &= \frac{N-b+a+2}{N}p_{a+2}(n) + \frac{b-a-1}{N}p_{a+1}(n)\\
p_{a+3}(n+1) &= \frac{N-b+a+3}{N}p_{a+3}(n) + \frac{b-a-2}{N}p_{a+2}(n)\\
&\vdots\\
p_{k}(n+1) &= \frac{N-b+k}{N}p_{k}(n) + \frac{b-k+1}{N}p_{k-1}(n)\\
&\vdots\\
p_{b}(n+1) &= \frac{N}{N}p_{b}(n) + \frac{1}{N}p_{b-1}(n)\\
\end{aligned}
\]
これを行列で表して対角化を施したりすると、上記のような一般式が得られることになる。

coupon_collectors_problem