プリキュアカレーのシール全30種コンプするには何箱必要か？その２～状態遷移モデルで考える

# はじめに

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

[動画](https://www.youtube.com/watch?v=BZxgchEPXis)の中で<b>しのけん</b>氏自身が述べていましたが，この種のコレクターズアイテムは最後のほうにはひたすら当たりを待つ忍耐ゲームとなります。全 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}$ となります。


```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontSize 16
  fontColor black
  fontStyle normal
}
</style>
title
<img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/45038e41-2267-47a1-96fc-79d9de2f795a.png>

図１　状態遷移確率
end title
@enduml
```

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

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

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

```math
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 箱だけでコンプリートできる確率は簡単に求められます。これは最後に確認します。

```math
\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# で作りました。各配列変数の意味は下記の通りです。

<table>
<caption>表１　配列変数の数学的意味</caption>
<thead>
<tr><th style="text-align:center">配列変数</th>
    <th style="text-align:center">数学的意味</th></tr>
</thead>
<tbody>
<tr></tr>
<tr><td style="text-align:center"><code>p[k]</code></td>
    <td style="text-align:center">$P_n(k + 1)$</td></tr>
<tr></tr>
<tr><td style="text-align:center"><code>q[k]</code></td>
    <td style="text-align:center">$(1 - n / 30)(n / 30)^{k}$</td></tr>
<tr></tr>
<tr><td style="text-align:center"><code>r[k]</code></td>
    <td style="text-align:center">$P_{n + 1}(k + 1)$</td></tr>
<tr><th></th><th></th></tr>
</tbody>
</table>

```C#: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;
	}
}
```

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

<details><summary><b>ビルド用のバッチファイル BUILD2.CMD はコチラ</b></summary>

※バッチファイルの文字コードは 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
```

</details>

# 計算結果

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

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontSize 16
  fontColor black
  fontStyle normal
}
</style>
title
<img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/1d53a349-7146-9471-378e-dfaff5045704.png>

図２　シール全 30 種をコンプするまでに食する箱数とその確率（線形スケール）
end title
@enduml
```

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

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

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

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontSize 16
  fontColor black
  fontStyle normal
}
</style>
title
<img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/85fd29a8-2ce6-a9db-2757-c95cf7298d2a.png>

図３　シール全 30 種をコンプするまでに食する箱数とその確率（対数スケール）
end title
@enduml
```

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

# まとめ

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

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

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

# おまけ

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

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

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontSize 16
  fontColor black
  fontStyle normal
}
</style>
title
<img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/b651f4ff-c2c0-8fb0-2480-7c5981abb4d1.png>

図４　丸美屋の「ちいかわ」ふりかけのシール全20種の場合
end title
@enduml
```

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

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontSize 16
  fontColor black
  fontStyle normal
}
</style>
title
<img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/fe80e883-4af8-a2eb-3b82-11548d38a0d3.png>

図５　永谷園のお茶漬けの東海道五拾三次カード全55種の場合
end title
@enduml
```


# 参考文献

- [プリキュアカレーのシールコンプするまで限界食いするも一生ダブり続けて絶望する漢 - しのけんの別腹 <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/ee21aa0e-add2-222e-dadc-d8ffad3be7a7.png">](https://www.youtube.com/watch?v=BZxgchEPXis)

- [プリキュアシリーズ - 丸美屋](https://www.marumiya.co.jp/product/series/character/precure/)

- [ちいかわシリーズ - 丸美屋](https://www.marumiya.co.jp/product/series/character/chiikawa/)

- [商品情報 - 永谷園](https://www.nagatanien.co.jp/product/)

- [アーラン分布 - wikipedia](https://ja.wikipedia.org/wiki/アーラン分布)

- [待ち行列理論 その1 - Qiita](https://qiita.com/ogata-k/items/ed53406a8e6010fe4a43)

- [ガンマ関数の形が階乗の一般化として現れることをポアソン過程から導く - Qiita](https://qiita.com/Ken-ichi_Hironaka/items/64c19f8651f7d0b6f429)