はじめに

前回の記事は出先で書いていたので，まず結論を得るために速攻でシミュレーションプログラムを作成して求めましたが，今回は数学の問題としてじっくり落ち着いて考えてみることにします。

動画の中でしのけん氏自身が述べていましたが，この種のコレクターズアイテムは最後のほうにはひたすら当たりを待つ忍耐ゲームとなります。全 30 種のうち 29 種までを集めた状態まで達すると，残る 1 種類を確率 $1/30$ で待ち続けることになるからです。その手前の 28 種までを集めた状態を考えると，確率 $2/30$ で待ち続けることになります。さらにその手前の 27 種を集めた状態では確率 $3/30$ で待ち続けます。この調子で初期状態まで遡って考えると何だか解けそうです。そう，いわゆる＊状態遷移モデル＊を考えます。

状態遷移モデルを考える

ゼロ種類揃えた状態（初期状態）から一種類揃えた状態への遷移確率は（当たり前ですが）100% です。次に二種類揃えた状態へ一箱で済む確率は $29/30$ ですが，$1/30$ の確率で同じ種類を引いてしまう可能性があります。二箱要する確率は $(29/30) \times (1/30)$ であり，三箱要する確率は $(29/30) \times (1/30)^2$，以降 $m$ 箱要する確率は $(29/30) \times (1/30)^{m - 1}$ となります。

以下，同様にして $n$ 種類揃えた状態から $n + 1$ 種類揃えた状態へ一箱で遷移できる確率は $(1 - n/30)$ であり，既に持っているものと重複する確率は $n/30$ ですから，$m$ 箱要する確率は $(1 - n/30) \times (n/30)^{m - 1}$ となります。

$n$ 種類集めるために要した箱数 $m$ の確率 $P_n(m)$ とおくと，初項の $P_1(m)$ は

P_1(m) =\left\lbrace
\begin{matrix}
1 & m \le 1 \\
0 & m > 1
\end{matrix} \right. \tag{1}

となり，以降，下記の漸化式によって求めることができます。

P_{n + 1}(m) = \left \lbrace \begin{matrix}
0 & m \le n \\
\displaystyle \sum_{k = 1}^{m - 1} P_{n}(m - k) \left( 1 - \frac{n}{30} \right) \left( \frac{n}{30} \right) ^{k - 1} & m > n
\end{matrix} \right . \tag{2}

最終的に $P_{30}(m)$ を求めるのが目標です。これを手計算で求めるのは大変そうですが，計算機を使えるのなら楽勝でしょう。

軽く計算してみる

全 30 種を重複なしで 30 箱だけでコンプリートできる確率は簡単に求められます。これは最後に確認します。

\frac{30}{30} \times \frac{29}{30} \times \frac{28}{30} \times \cdots \times \frac{3}{30} \times \frac{2}{30} \times \frac{1}{30} = \frac{30!}{30^{30}} \approx 1.29 \times 10^{-29} \tag{3}

実装プログラム

今回も C# で作りました。各配列変数の意味は下記の通りです。

表１　配列変数の数学的意味
配列変数	数学的意味
`p[k]`	$P_n(k + 1)$
`q[k]`	$(1 - n / 30)(n / 30)^{k}$
`r[k]`	$P_{n + 1}(k + 1)$

PRECURE2.CS

using System;
using System.IO;
class PRECURE2 {
	static	int		Main() {
		var	p = new double[1000];
		var q = new double[1000];
		var	r = new double[1000];
		p[0] = 1.0;
		for(int n = 1; n < 30; n++) {
			q[0] = 1.0 - (double)n / 30;
			for(int i = 1; i < q.Length; i++)
				q[i] = q[i - 1] * (double)n / 30;
			for(int i = 0; i < r.Length; i++)
				r[i] = 0;
			for(int i = 0; i < p.Length - 1; i++)
				for(int j = 0; i + j < q.Length - 1; j++)
					r[i + j + 1] += p[i] * q[j];
			for(int i = 0; i < p.Length; i++)
				p[i] = r[i];
		}
		for(int i = 0; i < p.Length; i++)
			Console.Out.WriteLine("{0}\t{1}", i + 1, p[i]);
		return 0;
	}
}

ビルド用のバッチファイルは前回とほぼ同じです。

ビルド用のバッチファイル BUILD2.CMD はコチラ

※バッチファイルの文字コードは Shift-JIS で保存して下さい。

BUILD2.CMD

@echo off
setlocal
rem ----------------------------------------------------------------------------
rem ターゲットプログラム
rem ----------------------------------------------------------------------------
set TARGET=precure2

rem ----------------------------------------------------------------------------
rem csc.exe にパスが通っているかどうかの確認
rem ----------------------------------------------------------------------------
set DOTNETPATH=
for %%I in ( csc.exe ) do set DOTNETPATH=%%~$PATH:I
if defined DOTNETPATH goto LETS_COMPILE

rem ----------------------------------------------------------------------------
rem .NET Framework のインストールパスを探す
rem ----------------------------------------------------------------------------
for /F "usebackq tokens=1,2,*" %%I in (
  `reg query "HKLM\SOFTWARE\Microsoft\NET Framework Setup\NDP\v4\Full" /v InstallPath`
) do (
  if "%%I"=="InstallPath" set "DOTNETPATH=%%K"
)
if defined DOTNETPATH goto APPEND_PATH

echo .NET Framework のインストールパスが見つかりません！
pause
exit /b

:APPEND_PATH
echo パスに %DOTNETPATH% を追加します。
set "PATH=%PATH%;%DOTNETPATH%"

:LETS_COMPILE
rem ----------------------------------------------------------------------------
rem コンパイル
rem ----------------------------------------------------------------------------
echo %TARGET%.cs をコンパイルします。
if exist %TARGET%.exe del %TARGET%.exe
csc -nologo -o -w:4 %TARGET%.cs

if exist %TARGET%.exe goto SUCCESS
echo %TARGET% のコンパイルに失敗しました！
pause
exit /b

:SUCCESS
echo %TARGET% のコンパイルに成功しました。
pause
exit /b

計算結果

計算結果を以下に示します。今回の計算は一瞬で終わりました。なお，確率ピーク（最頻値）は前回のシミュレーションより僅かに増えて 102 箱となりましたが，中央値は 113 箱，平均値は 120 箱と前回と一致する結果になりました。

前回のシミュレーションでは確率ゼロの点が生じるため，対数スケール表示できませんでしたが，今回の計算では対数スケールで表示させることができます。

また 30 箱でコンプリートできる確率も求めることができ，それは $1.29 \times 10^{-12} $ となりました。これは式$(3)$での手計算と一致しています。

プリキュアチャレンジするしのけん氏が一兆人いれば，わずか 30 箱でコンプリートできる人が一人誕生するかもしれないのです。ただし，同程度の確率で 800 箱オーバーになってしまう気の毒な人も誕生してしまいます。

※ 30 箱未満は確率ゼロになるのでグラフにプロットしていません。

まとめ

大したプログラムでもないので最初から本記事のプログラムを出せば良かったようにも思いますが，この種の問題の最大の課題は正解が分からないということです。なので，前回の記事のシミュレーションが一つの解答（候補）であり，今回の計算とも＊ほぼ＊一致していることから＊まあまあ＊合っているんじゃないかなと思っています。

最頻値が僅かに変わった点については，本記事の計算に僅かな誤りがあった可能性も否定できませんが，どちらかというと前回の記事のシミュレーションにおける乱数の品質のほうが怪しいように思います。

なお，プリキュアチャレンジにおける確率分布曲線は統計学におけるガンマ分布あるいはアーラン分布に近いように思いますが，筆者の数学レベルが低過ぎてよく分かりません。

おまけ

このプログラムは，数値を書き換えるだけで他の問題にも応用が利きます。

同じ丸美屋の「ちいかわ」ふりかけのシール全 20 種に適用した場合，最頻値は 60 個，中央値は 67 個，平均値（期待値）は 72 個となりました。

一方，永谷園のお茶漬けの東海道五拾三次カード全 55 種に適用した場合，最頻値は 220 個，中央値は 241 個，平均値（期待値）は 253 個となりました。ちなみに五十三次なのに何故 55 種かというと，あくまで五十三次とは起点の日本橋と終点の三条大橋を除いた途中の宿場町のことだからなんですね。浮世絵カードにはしっかり日本橋と三条大橋も含まれているので全 55 種となります。

参考文献