クラップス「Make 'Em All」の確率を計算する

カジノゲーム「クラップス」のサイドベットの一つである「Make 'Em All」という少し変わったルールの確率を、「確率モデルによる近似」と「動的計画法（DP）」の二通りの方針で計算する方法を紹介します。

## クラップスの「Make 'Em All」ベットとは？

今回計算するのは、「Make 'Em All」というベットです。ルールは以下の通りです。

- **勝ち条件**： シューター（サイコロを2つ投げる人）が出目の合計値`7`を出す前に、`2, 3, 4, 5, 6, 8, 9, 10, 11, 12`の全ての目を出す
- **負け条件**： 上記の全ての目を出す前に、`7`が出てしまう

このルールが複雑なのは、「特定の順番で出す必要はない」が「`7`が出る前に全て出す」必要がある点、そして対象となる目が複数ある点です。

## 確率計算へのアプローチ

この種の複雑な条件を持つ確率を計算する一つのアプローチとして、ここでは指数分布を用いた連続的な確率モデルと積分による近似計算を試みます。

### 確率モデルによる近似

計算を簡略化するため、以下のようなモデルを考えます。

1. 各サイコロの目の合計値（2から12）が初めて出現するまでにかかる試行回数 `T_i` が、それぞれ独立に指数分布に従うと仮定します。この指数分布のパラメータ `λ_i` は、それぞれの出目 `i` が1回の試行で出る確率（例：2の目が出る確率は `1/36`、3の目は `2/36` など）に比例すると考えます。具体的には、`λ_i`をその出目が出る確率そのものとします

1. これらの指数分布に従う時間 `T_i` はすべて互いに独立であると仮定します

この仮定は、本来離散的なサイコロの試行を連続的な時間で近似するものです。

### 勝利条件の数式化

「Make 'Em All」ベットに勝つ条件は、「7の出目 `T_7`が初めて出現するまでにかかる時間が、7以外の指定された全ての出目（`{2,3,4,5,6,8,9,10,11,12}`）がそれぞれ初めて出現するまでにかかる時間のうち最も遅いものよりも大きい」と言い換えられます。これを数式で表すと次のようになります。

```math
max(T_2,T_3,T_4,T_5,T_6,T_8,T_9,T_{10},T_{11},T_{12}) < T_7
```

私たちが計算したいのは、この条件が満たされる確率です。

```math
P(max(T_2,T_3,T_4,T_5,T_6,T_8,T_9,T_{10},T_{11},T_{12}) < T_7)
```

### 確率計算

それぞれの出目`i`が出る確率を`λ_i`とします。`T_i`は指数分布に従うと仮定したので、確率密度関数を用いて以下のように式変形します。

```math
P(max(T_2, \dots, T_{12}) < T_7) = \int_{0}^{ \infty } P(max(T_2, \dots, T_{12}) < t) λ_7 e^{ -λ_7 t} dt
```

また、`T_i`は互いに独立であると仮定したため、`P(max(T_2,…,T_12) < t)`は独立な事象の同時確率と同じです。

```math
P(max(T_2, \dots, T_{12}) < t) = P(T_2 < t) \times P(T_3 < t) \times \dots \times P(T_{12} < t)
```

```math
= (1 - λ_2 e^{-λ_2 t}) (1 - λ_3 e^{-λ_3 t}) \dots (1 - λ_{12} e^{-λ_{12} t})
```

これを元の式に代入して

```math
P(max(T_2, \dots, T_{12}) < T_7) = \int_{0}^{ \infty } (1 - λ_2 e^{-λ_2 t})  \dots (1 - λ_{12} e^{-λ_{12} t}) λ_7 e^{ -λ_7 t} dt
```

```math
= \frac {6} {36} \int_{0}^{ \infty } e^{ - \frac {6} {36} t} (1 - e^{ - \frac {1} {36} t})^2 (1 - e^{ - \frac {2} {36} t})^2 (1 - e^{ - \frac {3} {36} t})^2 (1 - e^{ - \frac {4} {36} t})^2 (1 - e^{ - \frac {5} {36} t})^2 dt
```

```math
= \frac {1} {6} \int_{0}^{ \infty } e^{ - \frac {1} {6} t} (1 - e^{ - \frac {1} {36} t})^2 (1 - e^{ - \frac {1} {18} t})^2 (1 - e^{ - \frac {1} {12} t})^2 (1 - e^{ - \frac {1} {9} t})^2 (1 - e^{ - \frac {5} {36} t})^2 dt
```

この積分を数値計算すると、約 0.0052577 が得られます。これは確率としておよそ 1/190 に相当します。

## 離散モデルによる正確な計算（動的計画法）

サイコロの目は離散的なので、離散的なモデルで考える方がより直接的で正確な値を求められます。ここで強力なツールとなるのが動的計画法（DP）です。DPは、大きな問題を小さな部分問題に分割し、それぞれの結果を記録・再利用することで効率的に解を求める手法です。

### 状態の定義

DPの状態として、「どの対象の目が出たか」という情報と、「その状態（まだ7が出ていない条件で、特定の目の組み合わせが達成されている状態）に至る確率」を管理します。

「どの目が出たか」は、11個の数字（2,3,4,5,6,7,8,9,10,11,12）に対応するビット列（ビットマスク）で表現できます。例えば、0番目のビットが2、1番目のビットが3、…、5番目のビットが7、…、10番目のビットが12に対応するとします。

DPテーブルの要素 `dp[mask]` は、「それまでのロールで7が出ることなく、`mask` で示される目の組み合わせが達成された確率」を表します。（実装では、ロール回数のパリティを使ってDPテーブルを使いまわします）

### 状態遷移

現在の状態（あるマスク `current_mask` が確率 `p_current` で達成されている）から、次のサイコロの目 `roll` が出た場合（その出目が出る確率を `P(roll)` とします）、新しい状態のマスク `next_mask = current_mask | (1 << roll_bit)` に遷移します。その遷移後の状態に至る確率には `p_current * P(roll)` が加算されます。

注意すべきポイント：

1. 7の処理
    - もし現在の状態 `current_mask` が既に7のビットを含んでいたら、その経路は既に負けなので、それ以上計算を進めません
    - 現在のロールで7が出た場合、その経路も負けとなります。その確率を「勝ちにつながる次の状態」に加算してはいけません
1. 勝利条件
    - `mask_completed` を、7以外の全ての必須の目が揃った状態（例：`0b11111011111`）と定義します
    - あるロールの結果、状態が `mask_completed` になったら、その分の確率を最終的な答え `answer` に加算します
1. 既に勝利した状態の処理
    - 一度 `mask_completed` に到達した状態（その確率が `answer` に加算された後）は、それ以上ロールを続けても「新たに勝利する」わけではないので、次のDP遷移の計算からは除外します。これにより、確率を二重に加算することを防ぎます

### 実装例

以下はC++を用いた実装の一例です。

```c++
#include <iostream>
#include <iomanip>

using namespace std;

int main() {
    // dp[ロールのパリティ(0か1)][出た目のビットマスク]
    // この配列は、各状態（特定の目の組み合わせが出た状態）になる確率を保持します。
    // ただし、7がまだ出ていないという条件下での確率です。
    double dp[2][1 << 11] = {}; // 11ビット: 2,3,4,5,6, (7), 8,9,10,11,12 に対応

    // 各目が出る確率 (2から12まで)
    // dice_prob[0] は P(2), dice_prob[5] は P(7), dice_prob[10] は P(12)
    double dice_prob[] = {
        1.0 / 36.0, // 2
        2.0 / 36.0, // 3
        3.0 / 36.0, // 4
        4.0 / 36.0, // 5
        5.0 / 36.0, // 6
        6.0 / 36.0, // 7 (インデックス 5)
        5.0 / 36.0, // 8
        4.0 / 36.0, // 9
        3.0 / 36.0, // 10
        2.0 / 36.0, // 11
        1.0 / 36.0, // 12
    };

    const int SEVEN_INDEX = 5; // 7に対応するdice_prob配列のインデックス
    const int SEVEN_BIT = (1 << SEVEN_INDEX); // 7に対応するビットマスク中のビット

    // 勝利条件となるマスク: 2-12の全ての目(7を除く)が出た状態
    // (1 << 11) - 1 は全ての11ビットが立った状態 (0b11111111111)
    // SEVEN_BIT と XOR することで、7のビットだけを反転させる (この場合は0にする)
    int mask_completed = ((1 << 11) - 1) ^ SEVEN_BIT; // 例: 0b11111011111
   
    // 初期状態: まだ何もサイコロを振っていないので、出た目はなし(マスク0)。この確率は1。
    dp[0][0] = 1.0;

    double answer = 0.0; // 最終的に求める「Make 'Em All」が成功する確率
    int N = 1000; // 最大ロール回数。この回数までに決着がつかない確率は非常に小さいと仮定。

    for (int i = 0; i < N; i++) { // i を現在のロール回数-1 とする (0からN-1まで)
        int current_dp_layer = i % 2;     // 現在の確率分布が格納されているDP層
        int next_dp_layer = (i + 1) % 2;  // 次のロール後の確率分布を計算するDP層

        // 次のDP層を初期化
        for (int state_idx = 0; state_idx < (1 << 11); state_idx++) {
            dp[next_dp_layer][state_idx] = 0.0;
        }

        // 現在の各状態 (state) から遷移を計算
        for (int state = 0; state < (1 << 11); state++) {
            // もし現在の状態 'state' に既に7のビットが含まれていたら、
            // このパスは7が出て負けが確定しているので、これ以上計算しない。
            if ((state & SEVEN_BIT) != 0) {
                continue;
            }

            // もし現在の状態が既に勝利条件を満たしていたら、
            // その確率は前のロールで answer に加算済みのはずなので、ここでは何もしない。
            // (このパスはここで終わり、新たな勝利は生まない)
            if (state == mask_completed) {
                continue;
            }
            
            // 次のロールで出る可能性のある全ての目 (roll_idx: 0から10) について計算
            for (int roll_idx = 0; roll_idx < 11; roll_idx++) {
                double prob_of_this_roll = dice_prob[roll_idx]; // 今回振った目が出る確率
                int bit_for_this_roll = (1 << roll_idx);       // 今回振った目に対応するビット

                // もし今回振った目が7だったら (roll_idx が SEVEN_INDEX と同じ)
                if (roll_idx == SEVEN_INDEX) {
                    // 7が出たので、このパスは負け。
                    // dp[next_dp_layer] には何も加算しない (勝ちにはつながらない)。
                    // この dp[current_dp_layer][state] * prob_of_this_roll の確率質量はここで消滅する。
                } else {
                    // 7以外の目が出た場合
                    int new_state = state | bit_for_this_roll; // 新しい出目状態のマスク

                    // 7以外の目が出て、状態が更新される確率を次層のDPに加算
                    dp[next_dp_layer][new_state] += dp[current_dp_layer][state] * prob_of_this_roll;
                }
            }
        }
        // このロール(i+1回目)の結果、新たに勝利条件(mask_completed)を達成した確率をanswerに加算
        // dp[next_dp_layer][mask_completed] には、このロールで7以外の目を出し、
        // 結果としてmask_completedに到達した経路の確率が集まっている。
        answer += dp[next_dp_layer][mask_completed];
    }

    cout << fixed << setprecision(7) << answer << endl;

    return 0;
}
```

上記のC++コードを実行すると、近似計算と同様に約 0.0052577 という結果が得られます。