エモクロアTRPGロールシャッハシンドロームのロスト率を求めてみた～導出編～

この記事は下記記事の導出編である。
https://note.com/terapotan/n/n593c80adbd7b

生存確率を求める

前提

　途中のループにおいて、PLの選択によってはどちらか一方が傷を負った状態のまま、あえてループを脱することも可能ではあるが、今回の確率計算においてそのような場合は考慮しない。両方無傷を引いたときか、もしくは共鳴レベルが10以上になったときのみ、ループを脱するものとする。
　このシナリオにおいて「ロストする」とは、逸脱時「胸、頭、全身」のいずれかであってなおかつ技能判定でファンブルを出した場合（つまり胸、頭、全身の部位を喪失した場合。さすがにこの部位を失って生存を続けるのは難しいだろう）とする。足や腕の喪失はロストに含まないものとする。
　なお、共鳴感情の一致や残響の効果は考慮しないものとし、進行は公式エモクロアルールブック、ロールシャッハシンドローム公式シナリオに則って行われるものとする（その卓独自の進行は一切考慮しない）。

　以下は記号の説明である。
$p_{dm}(d,j,s)$は$d$DM<=$j$を行い成功数が$s$となる確率を表す。なお$-d\leqq s \leqq 2d$である。$P_{dm}(d,j,s)$とは

$$ P_{dm}(d,j,s)=\sum_{i=s}^{2d}{p_{dm}(d,j,i)} $$
である。すなわち$d$DM<=$j$を行った成功数が$s$以上となる確率を表す。

$S_{1}$とは作中の共鳴判定を行って成功する確率を表している（sは共鳴を意味するSympathyの頭文字である。1は成功数を表している）。$S_{1}$は$P_{dm}$を使って
$$ S_1(k)=P_{dm}(k,9,1) $$
と書ける。ここで$k$とは〈∞共鳴〉の技能レベルを表している。$S_{3}$は作中の共鳴判定を行ってトリプル以上の成功を出す確率を表している。$S_{1}$と同様に$S_3$も$P_{dm}$を使って次のように書ける。
$$ S_3(k)=P_{dm}(k,9,3) $$

　$P_{ain}(L)$とは、痛み表のロールにおいて、共鳴者と(DLが担当する)NPC、少なくともどちらか一方が無傷でない出目を引いてしまう確率である。$L$とは無傷であるような出目の最大数を表す。以降無傷レベルと呼称する。痛みレベルが$L$であるとき、痛み表の$1,2,\cdots,L$がすべて無傷である。痛みレベル$L$はシナリオ記載の条件により、$L\leqq 5$である。

　$P_{ain}(L)$はどちらも無傷である確率から1を引いたものに等しい。したがって
$$ P_{ain}(L)=1-\frac{L^2}{64} $$
である。

　試行過程は次の通りである。

1. $k=1$,$L=1$からスタートする。
2. 共鳴判定を4回行う。詳細は以下。
      1. $S_1(k)$の確率でkに1加算する。
      2. 手順1に戻る。これを4回繰り返す。
3. 痛み表の判定を行う。手順の詳細は以下。
      1. $S_1(k)-S_3(k)$の確率で$k$にのみ1加算する。共鳴判定の成功。
      2. $S_3(k)$の確率で$k$と$L$にそれぞれ1加算する。トリプル以上の成功。
      3. $1-S_1(k)$の確率で$k$と$L$の値にそれぞれ0を加算する（つまり変動なし）。共鳴判定失敗。
         1. 以上の三手順は、共鳴判定を行って成功したら共鳴値を1上昇させ、そのうちトリプル以上の成功だった場合無傷レベルを1上昇させる、という手順を$S_1(k),S_3(k)$を使って書き表したものである。
      4. $1-P_{ain}(L)$の確率で試行を終了する。この形での試行終了の事象を$\alpha$とし、その確率を$p_{\alpha}$とする。$p_{\alpha}=1-P_{ain}(L)$でないことに注意せよ。詳細は後述するが、すべての$k$と$L$の場合においてその都度確率を求め合算しなくては真の$p_{\alpha}$は求められない。$\alpha$は双方がどちらも無傷を引いてループを脱出する、つまり無傷で生還する事象に対応する。
      5. $k\geqq10$のとき試行を終了する。この形での試行終了の事象を$\beta$とし、その確率を$p_{\beta}$とする。$\beta$は逸脱によるループの強制終了であって$\alpha$に該当しない事象に対応する。$\beta$において少なくとも一方が「胸、頭、全身、無傷」以外で終わるとき、この事象を$\gamma$とし、その確率を$p_{\gamma}$とする。$\beta$においてどちらも「胸、頭、全身」で終わるとき、この事象を$\delta$とし、その確率を$p_{\delta}$とする。
4. その定義から明らかに$p_{\alpha}+p_{\beta}=1$が成り立つ。

$p_{dm}(d,j,s)$をプログラムで求める

$n_{v,j}$を求める

　$p_{dm}(d,j,s)$を求める前に1個のサイコロを振ったときの各成功数の場合の数を求める。前者の値を求めるためには後者の値が必要である。以下1個のサイコロを振り判定数が$j$であったときに、成功数が$v$であるような場合の数を$n_{v,j}$と置く。
　$j$の値で場合分けを行って求める。

　$j=1$の場合、各出目における成功数は次の通りである。

　よって$n_{-1,1}=1,n_{0,1}=8,n_{1,1}=0,n_{2,1}=1$である。
　$j=10$の場合、各出目における成功数は次の通りである。

　よって$n_{-1,10}=0,n_{0,10}=1,n_{1,10}=8,n_{2,10}=1$である。

　証明は省略するが、$2\leqq j \leqq 9$のとき$n_{-1,j}=1,n_{0,j}=9-j,n_{1,j}=j-1,n_{2,j}=1$である。例えば$j=5$のとき

であるから$n_{-1,5}=1,n_{0,5}=4,n_{1,5}=4,n_{2,5}=1$であるがこれは先ほど挙げた式の値と一致する。

$p_{dm}(d,j,s)$をプログラムで求める

　$dp[d][s]$を$d$DM<=$j$において成功数$s$となる場合の数と定義する。ここで$p_{dm}(d,j,s)$は$dp[d][s]$を用いて次のように書ける。
$$  p_{dm}(d,j,s)=\frac{dp[d][s]}{10^d} $$
　このことを踏まえて2DM<=5において成功数4となる確率を求めてみよう。

　今ダイスを0個振ったとき、合計数が0になる場合の数は1である。このことから$d=0$の行の値は次のように埋められる。

　$d=1$の行を埋めていきたい。まず$dp[0][-2]$を考える。ここから$dp[1][-2]$に到達するにはダイスを1つ振って成功数0を出さなければならない。すなわち
$$ dp[1][-2]\leftarrow dp[1][-2] + dp[0][-2] \times n_{0,5} $$
である。ここで$\leftarrow$は代入を表す。同様にして$dp[0][-1]$から$dp[1][-2]$に到達するにはダイスを2つ振って成功数-1を出さなければならない。すなわち
$$ dp[1][-2]\leftarrow dp[1][-2] + dp[0][-1] \times n_{-1,5} $$
である。$dp[-1,1]$についても同じである。$dp[-2][0],dp[-1][0],dp[0][0]$からの到達を考えてそれぞれで場合の数を計算して足し合わせればよい。
　以上のことをまとめて式にしたのが下式である。$D$DM<=$j$のとき$dp[d][s]$の値を埋めるための公式として使える。
$$ dp[d][s]=\sum_{k=-1}^{2}v_{k,j}dp[d-1][s-k] $$
　ただし$s-k<-D$か$s-k>2D$のときは$dp[d][s-k]=0$とする。例えば$dp[1][-1]$は
$$ dp[1][-1]=v_{-1,5}dp[0][0]=1 $$
　である。$d=1,d=2$の列を埋めると最終的に次のような値になる。なお実際のプログラムでは遷移元を起点として遷移先の要素に対して場合の数を足していく形式になっているが、最終的に算出される値はどちらの方法でやっても変わらない。加算する順番が異なるだけである。

$p_{\alpha},p_{\beta}.p_{\gamma},p_{\delta}$をプログラムで求める

　ここではプログラムの概略を記すのみとする。詳細についてはソースコードを参照のこと。

1～4週目における共鳴判定

　4週目の共鳴判定が終了した時点で$k$は$1,2,3,4,5$のいずれかの値を取る。ここで$k$がそれぞれの値になる確率を求めたい。
　方針としてはおおよそ次の通りである。学校の校庭に人を100人集める。そして学校の校庭を$k=1,\cdots,5$の五つに区分し、最初は$k=1$に100人配置しておく。このたとえ話における「一人」は本状況におけるプレイヤー一人に相当する。
　各々は事前に定めておいた確率に従って他の区分へと移動する。例えば1回目の試行では$100S_1(k)$人の人間が$k=2$へ移動し、$100-100S_1(k)$人の人間が$k=1$に留まる。移動が完了したら1回目の試行は終了である。この状況は1回目の共鳴判定終了時点で、$S_1(k)$の確率で$k=2$、$1-S_1(k)$の確率で$k=1$となる状況に相当する。

　2回目の試行も同様にして行う。1回目の試行終了時点で$k=1$にいる人間のうち、$100S_1(k)$%の人間は$k=2$へ移動し、$100(1-S_1(k))$%の人間は$k=1$に留まる。1回目の試行終了時点で$k=2$にいる人間のうち、$100S_1(k)$%の人間は$k=3$へ移動し$100(1-S_1(k))$%の人間は$k=2$に留まる。2回目の試行終了時点での最終的な人数は、$k=1$は$k=1$を動かずにいた人間の人数に等しく、$k=2$は$k=1$から$k=2$へ移動してきた人数と$k=2$を動かずにいた人数の和に等しく、$k=3$は$k=2$から$k=3$へ移動してきた人数に等しい。これで2回目の試行は終了である。

　後は同様の手順を4回繰り返せば、無事求めたい確率が求まる。適当に遷移確率を定めて、一度手動で上記アルゴリズムをやってみると理解が深まるかと思う。文章量に反してやっていることは意外にシンプルである。

5週目以降の痛み判定

　$p_{\alpha}$を求めるには、すべての$k$と$L$の組み合わせにおいて$\alpha$が生じる確率を求めればよい。なぜならいかなる過程を踏むにせよ、いずれかの$k$と$L$の組み合わせには到達するからである。これは$p_{\beta}$も同様である。以下dp[a][b]とは、$k=a$かつ$L=b$以外の状態から該当の状態に到達する確率を表すものと約束する。
　
　$k=a$かつ$L=b$であるとき、$k$と$L$の増減パターンは以下の三種類である（再掲）。

1. $S_1(k)-S_3(k)$の確率で$k$にのみ1加算する。共鳴判定成功。ここではパターン1と呼称する。
2. $S_3(k)$の確率で$k$と$L$にそれぞれ1加算する。トリプル以上の成功。ここではパターン2と呼称する。
      1. $L+1$が5より大きい場合、$L$に1は加算しない。
3. $1-S_1(k)$の確率で$k$と$L$の値にそれぞれ0を加算する（つまり変動なし）。共鳴判定失敗。ここではパターン3と呼称する。

　$k=a,L=b$に到達した上で、複数回ダイスを振っても、なお$k=a,L=b$に留まっている確率$p_c$を求める。同じ$k,L$に留まり、なおかつ$\alpha$にもならない（双方無傷とならない）から単純に考えれば
$$ p_c=dp[a][b] + dp[a][b](1 - S_1(k))P_{ain}(L) $$
　のようにすればいいように思う。しかしこの考えは誤りである。
　なぜなら、もう一度共鳴判定に失敗し$\alpha$にもならず同じ$k,L$に留まった確率、さらにもう一度続けて共鳴判定に失敗し$\alpha$にもならず同じ$k,L$に留まった確率……が考慮されていないためである。理論上無限回同じ場所に留まり続けられるのであって、確率を求めるにはその無限回の試行を考慮しなければならない。
　すなわち次のような式となる。
$$ \begin{align*} p_c&=dp[a][b] + dp[a][b](1 - S_1(k))P_{ain}(L)+dp[a][b]\{(1 - S_1(k))P_{ain}(L)\}^2+dp[a][b]\{(1 - S_1(k))P_{ain}(L)\}^3+\cdots\\ &=dp[a][b]\left(1+(1 - S_1(k))(P_{ain}(L)+\{(1 - S_1(k))P_{ain}(L)\}^2+\{(1 - S_1(k))P_{ain}(L)\}^3+\cdots\right) \end{align*} $$
　これは初項1、公比$(1 - S_1(k))P_{ain}(L)$の無限等比級数である。よって無限等比級数の和の公式から
$$ p_c=\frac{dp[a][b]}{1-(1 - S_1(k))P_{ain}(L)} $$
となる。

　以上を前提として踏まえた上で、$p_{\alpha},p_{\beta}$を求めることを考える。
　dp[a][b]の値が与えられているものとする。1～4週目における共鳴判定の手順があるから、値が与えられている要素があるのは自明である。
　$p_\alpha$から考える。ダイス結果のパターンは既に示した通り3通りしかないのだから、その3通りについて$\alpha$が起こる確率を考えて、足し合わせればよい。
　パターン1において$\alpha$が起こる確率$p_{\alpha1}$は、(ここに到達し複数回ダイスを振っても同じ場所に留まり続ける確率) X (共鳴判定に成功する確率) X (双方無傷を引く確率)　であるから
$$ p_{\alpha1}=p_cS_1(k)(1-P_{ain}(L)) $$
である。
　パターン2において$\alpha$が起こる確率$p_{\alpha2}$は、（ここに到達し複数回ダイスを振っても同じ場所に留まり続ける確率）X (トリプル以上の成功が共鳴判定で出る確率) X (双方無傷を引く確率)である。ただし、ここで双方無傷を引く確率は現在の無傷レベルから、一つ上がったところで計算を行わなくてはならない（$L+1$が5を超えていた場合はこの限りではない）。数式で表すと
$$ p_{\alpha2}=p_cS_3(k)(1-P_{ain}(L+1)) $$
である。
　パターン3において$\alpha$が起こる確率$p_{\alpha3}$は、(ここに到達し複数回ダイスを振っても同じ場所に留まり続ける確率) X (共鳴判定が失敗する確率) X (双方無傷を引く確率)　であるから
$$ p_{\alpha3}=p_c(1-S_1(k))(1-P_{ain}(L)) $$
　である。$k=a,L=b$で$\alpha$が起こる確率は、これら三つの値を足し合わせれば求められる。$p_{\alpha}$は到達可能なすべての$k,L$の場合において、上三つの値を足し合わせれば求められる。厳密にはこれ以外にも$k$が最大値に到達した時点での$\alpha$が生ずる確率を求め加算する必要があるが、その辺りの詳しいことに関してはソースコードを参照してほしい。

　$p_{\beta},p_{\gamma},p_{\delta}$についても$p_{\alpha}$と同じように愚直に確率を足し合わせていくだけであるから、説明を省略する。具体的な手順についてはソースコードを参照してほしい。
　なお$p_{\gamma},p_{\delta}$の導出方法は、$p_{\beta}$の導出とほぼ同じだが使用している関数が異なることに留意せよ。

　以上の処理が書かれたソースコードを下記に掲載する。

      def get_nv(j):
    """j (1～10) に応じた n_v の辞書を返す"""
    if j == 1:
        return {-1: 1, 0: 8, 1: 0, 2: 1}
    elif j == 10:
        return {-1: 0, 0: 1, 1: 8, 2: 1}
    else:  # 2 <= j <= 9
        return {-1: 1, 0: 9 - j, 1: j - 1, 2: 1}

def probability_of_dice(d, j, s):
    """
    d : サイコロの個数
    j : 成功判定の閾値 (1～10)
    s : 求める成功数の合計
    """
    if s < -d or s > 2*d:
        return 0.0
        
    nv = get_nv(j)

    #成功数にはマイナスも含む。配列の仕様上負のインデックスは取れない。
    #そこで成功数全体にdを足すことで、成功数が0始まりになるようにしている。
    offset = d
    size = 3 * d + 1
    dp = [0] * size
    #0個のサイコロを振って成功数が0となる場合の数は1通りである。これを初期値として入力する。
    dp[0 + offset] = 1

    for i in range(d):
        next_dp = [0] * size
        for t in range(-i, 2*i + 1):
            cur = dp[t + offset]
            if cur == 0:
                continue
            for v, cnt in nv.items():
                if cnt == 0:
                    continue
                new_t = t + v
                next_dp[new_t + offset] += cur * cnt
        dp = next_dp

    total_ways = dp[s + offset]
    return total_ways / (10 ** d)

def dice_probability_greater_successes(d,j,s):
    sum_probability = 0
    for i in range(s,2*d+1):
        sum_probability += probability_of_dice(d,j,i)
    return sum_probability

def return_problem_probability(max_k, P_dm_func):
    L_MAX = 5
    
    def pain_prob_func(L):
        """無傷レベルLにおける痛み表ロールで、で少なくとも一方が無傷を引いていない確率"""
        return (64 - L**2) / 64

    def pain_gamma_func(L):
        """無傷レベルLにおける痛み表ロールで少なくとも一方が「胸、頭、全身、無傷」以外の出目を引く確率"""
        return 1-(3+L)**2 /64
    
    def pain_delta_func(L):
        """無傷レベルLにおける痛み表ロールでどちらも「胸、頭、全身」の出目を引く確率"""
        return 9/64
    
    #max_k行、L_MAX列の配列を作成し、全要素0で初期化する
    dp = [[0.0] * (L_MAX + 1) for _ in range(max_k + 1)]

    # 4週目の共鳴判定終了時の、各共鳴レベルの確率を計算
    current_k_dist = {1: 1.0}
    for _ in range(4):
        next_k_dist = {}
        for k, prob in current_k_dist.items():
            s1 = P_dm_func(k, 1)
            k_succ = min(k + 1, max_k)
            next_k_dist[k_succ] = next_k_dist.get(k_succ, 0.0) + prob * s1
            next_k_dist[k] = next_k_dist.get(k, 0.0) + prob * (1 - s1)
        current_k_dist = next_k_dist

    for k, prob in current_k_dist.items():
        dp[k][1] = prob

    p_alpha = 0.0
    p_beta = 0.0
    p_gamma = 0.0
    p_delta = 0.0

    # 5週目以降の計算
    for k in range(1, max_k + 1):
        for L in range(1, L_MAX + 1):
            #dp[k][L]とは他のkとLからそのkとLに到達する確率である
            curr_p = dp[k][L]

            p1 = P_dm_func(k, 1)   # 成功数1以上の確率
            p3 = P_dm_func(k, 3)   # 成功数3以上の確率

            prob_nothing = 1.0 - p1          # kもLも増えない確率
            prob_k_plus_only = p1 - p3       # kのみ増加する確率
            prob_k_and_l_plus = p3           # kとL両方増加する確率

            #dp[a][b]に遷移した後、複数回共鳴判定を行ってもなお、同じk=a,L=bに留まり続けている確率
            w = curr_p / (1.0 - prob_nothing * pain_prob_func(L))

            # 共鳴レベルが10（最大値）に達している場合
            # このターンは逸脱するターンであり、痛み表のロールを行って（=p_alpha,p_betaの算出）計算を終了する。
            if k >= max_k:
                p_alpha += w * (1 - pain_prob_func(L))
                p_beta += w * pain_prob_func(L)
                p_gamma += w * pain_gamma_func(L)
                p_delta += w * pain_delta_func(L)

                continue
            else: #まだ共鳴レベルが9以下である場合

                p_alpha += w * (prob_nothing * (1 - pain_prob_func(L)))
                p_alpha += w * (prob_k_plus_only * (1 - pain_prob_func(L)))
                #L==L_MAXのときそれ以上Lが増加することはないから、L+1であるような場面であったとしても、Lにする必要がある。
                p_alpha += w * (prob_k_and_l_plus * (1 - (pain_prob_func(L+1) if L < L_MAX else pain_prob_func(L))))

                if k + 1 >= max_k:
                    p_beta += w * (prob_k_plus_only * pain_prob_func(L))
                    p_beta += w * (prob_k_and_l_plus * (pain_prob_func(L+1) if L < L_MAX else pain_prob_func(L)))
                    
                    p_gamma += w * (prob_k_plus_only * pain_gamma_func(L))
                    p_gamma += w * (prob_k_and_l_plus * (pain_gamma_func(L+1) if L < L_MAX else pain_gamma_func(L)))

                    p_delta += w * (prob_k_plus_only * pain_delta_func(L))
                    p_delta += w * (prob_k_and_l_plus * (pain_delta_func(L+1) if L < L_MAX else pain_delta_func(L)))                    
                    
                else:
                    dp[k+1][L] += w * (prob_k_plus_only * pain_prob_func(L))
                    if L < L_MAX:
                        dp[k+1][L+1] += w * (prob_k_and_l_plus * pain_prob_func(L+1))
                    else:
                        dp[k+1][L] += w * (prob_k_and_l_plus * pain_prob_func(L))

    return p_alpha, p_beta,p_gamma,p_delta

def sample_p(d, s):
    # サイコロd個、閾値9で成功数s以上
    return dice_probability_greater_successes(d,9,s)

MAXK = 10
a, b, c, d = return_problem_probability(MAXK, sample_p)

print(f"--- 設定: MAXK={MAXK}, L_MAX=5 ---")
print(f"パターンα (失敗) 確率: {a:.6%}")
print(f"パターンβ (成功) 確率: {b:.6%}")
print(f"合計: {a+b:.6%}")
print(f"パターンγ 確率: {c:.6%}")
print(f"パターンδ 確率: {d:.6%}")
    

シミュレーションとの整合

　無傷で生還できる確率が$p$であるとする。無傷で生還した場合1、そうでない場合に0となるような確率変数$X$を考えたとき、確率変数$X$はパラメータ$p$のベルヌーイ分布に従う。
　ここで$n$回試行したときの値の平均値を$\bar{X_n}$とする。中心極限定理から、$n$が十分大きいとき下式で表される確率変数$Z$は標準正規分布$\mathcal{N}(0,1)$に従う。
$$ Z=\frac{\sqrt{n}(\bar{X_n}-p)}{p(1-p)}\sim\mathcal{N}(0,1) $$
　真の値$p$からの誤差を$\epsilon$としたとき、$\bar{X_n}$が$p-\epsilon$から$p+\epsilon$に収まる確率$A$は、次のように求められる。
$$ A=\phi\left(\frac{\sqrt{n}((p+\epsilon)-p)}{p(1-p)}\right)-\phi\left(\frac{\sqrt{n}((p-\epsilon)-p)}{p(1-p)}\right) $$
　ただし$\phi(x)$は標準正規分布の累積分布関数で
$$ \phi(x)=\int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}}\exp(-y^2/2)dy $$
　である。

　今理論値は$p=0.78111257$である。また100万回のシミュレーション結果である$\bar{X}$は$0.781218$である。このとき$\epsilon=0.001$としてシミュレーション結果が$p±\epsilon$に収まる確率は99.999999%と求まる。所定の条件下でシミュレーションを行うと$\bar{X}$はほぼ確実にこの範囲内に収まるという事実を示しており、実際$\bar{X}$はこの範囲内に収まっている。
　この事実は理論値とシミュレーション結果の正しさを直接立証するものではないが、少なくとも現時点では理論値とシミュレーション結果に矛盾はないことを示している。

　シミュレーションに用いたコードは次の通りである。

      import random
import collections

#痛み表のうち、unHarmeds個が無傷であったとき
#1d8を振り、「どちらも」出た目が無傷であった場合はTrue,それ以外の場合はFalseを返す。
#1<=unHarmeds<=5以外の範囲が指定された場合、RuntimeErrorを発生させる。
def decideAllPlayersUnharmed(unHarmeds,playerDice,npcDice):
    #範囲外の場合
    if not(1<=unHarmeds<=5):
        raise RuntimeError
    
    if(playerDice <= unHarmeds) and (npcDice <= unHarmeds):
        return True
    else: 
        return False
    

#LevelDM<=Powerの成功値を求める。
#乱数を用いるため、結果は毎回異なる。
#0をペナルティダイス、1をボーナスダイスとする。
#Powerは1以上9以下、Levelは1以上である必要がある。これ以外の範囲の整数を指定した場合例外を発生させる。
def calcSuccessfulValue(Level,Power):
    #範囲外の場合
    if not(1<=Power<=9) or Level<=0:
        raise RuntimeError
    
    diceValueList = [random.randint(1,10) for i in range(Level)]
    successfulValue = 0
    
    
    for i in diceValueList:
        #ペナルティダイス
        if i == 10:
            successfulValue -=1 
            continue
        
        #ボーナスダイス
        if i == 1:
            successfulValue += 1
        
        if i <= Power:
            successfulValue += 1
    
    return successfulValue


#脱出できたら"Success"を、脱出できなかったら"Failure"を返す。（どちらも文字列）
def decideExit():
    resonanceValue = 1

    #痛み表更新までの間のループ（1週から4週）
    for i in range(4):
        #共鳴判定に成功した場合、共鳴レベルを1上昇させる
        if calcSuccessfulValue(resonanceValue,9) >= 1:
            resonanceValue += 1

    unHarmedValue = 1
    exitFlag = False #trueなら脱出成功、falseなら脱出失敗
    #4週以降
    while resonanceValue <= 9:#逸脱するまでループを続行
        #共鳴判定を行う。ここでトリプル以上の成功が出た場合、痛み表のうち一つを無傷に変える
        successfulValue = calcSuccessfulValue(resonanceValue,9)
        #成功の場合
        if successfulValue >= 1:
            resonanceValue += 1
            #トリプル以上の成功の場合、痛み表のうち一つを無傷に変える
            if successfulValue >= 3 and unHarmedValue < 5:
                unHarmedValue += 1

        playerDice = random.randint(1,8)
        npcDice = random.randint(1,8)
        if decideAllPlayersUnharmed(unHarmedValue,playerDice,npcDice):
            exitFlag = True
            break

    if exitFlag :
        return "Mukizu"
    elif (playerDice >= 6) and (npcDice >= 6):
        return "playerAndNpcYabai"
    elif (playerDice >= 6):
        return "playerYabai"
    elif (npcDice >= 6):
        return "npcYabai"
    else:
        return "Failure"


successValueList = []

SIZE=1000000
for i in range(SIZE):
    successValueList.append(decideExit())

countDict=collections.Counter(successValueList)
allCount=len(successValueList)

for key in list(countDict.keys()):
    print(str(key) + ":" + str(100*countDict[key]/allCount)+"%")

#ループ終了後、ファンブル処理
pfc=0.25#PCがファンブルを出す確率
pfd=0.25#DLがファンブルを出す確率

for i in range(SIZE):
    if successValueList[i] == "playerYabai":
        #ファンブルを引いた場合
        if (True if random.random() <= pfc else False):
            #PCは死亡する
            successValueList[i] = "SIBOU"
    elif successValueList[i] == "npcYabai":
        #ファンブルを引いた場合
        if (True if random.random() <= pfd else False):
            #PCは死亡する
            successValueList[i] = "SIBOU"  
    elif successValueList[i] == "playerAndNpcYabai":
        #プレイヤーかDLがファンブルを引いた場合
        if (True if random.random() <= pfc else False) or (True if random.random() <= pfd else False):
            #PCかDLは死亡する
            successValueList[i] = "SIBOU"  

countDict=collections.Counter(successValueList)
allCount=len(successValueList)
print()
for key in list(countDict.keys()):
    print(str(key) + ":" + str(100*countDict[key]/allCount)+"%")

    

      import math
from scipy.stats import norm


p=0.78111257
sigma=p*(1-p)
barX=0.781218
n=1000000
epsilon=0.001

z_lower=(math.sqrt(n)*((p-epsilon)-p))/sigma
z_upper=(math.sqrt(n)*((p+epsilon)-p))/sigma

lower=norm.cdf(z_lower)
upper=norm.cdf(z_upper)

print(upper - lower)
    

ロスト確率を求める

理論値

　PCもしくはDLが操作するNPC、少なくともどちらか一方がロストする確率を求める。
　まず逸脱時に少なくともどちらか一方（これは両方そうであった場合も含まれる）が「胸、頭、全身」の出目を引いて終わる確率は、$p_{\beta}-p_{\gamma}$である。このうち、どちらも「胸、頭、全身」の出目を引いて終わる確率は$p_{\delta}$であるから、どちらか一方だけが「胸、頭、全身」の出目を引いて終わる確率（これは両方そうであった場合は含まれない）は、$(p_{\beta}-p_{\gamma})-p_{\delta}$である。
　PCの運動系技能あるいは生存系技能でファンブルを出す確率を$p_{fc}$とし、DLが操作するNPCの運動系技能あるいは生存系技能でファンブルを出す確率を$p_{fd}$とする。少なくともどちか一方が「胸、頭、全身」の出目を引いて終わる場合とは、次の3パターンを指す。

1. PC:「胸、頭、全身」以外の出目　DLNPC:「胸、頭、全身」の出目
2. PC:「胸、頭、全身」の出目　DLNPC:「胸、頭、全身」以外の出目
3. PC:「胸、頭、全身」の出目　DLNPC:「胸、頭、全身」の出目

　これらの場合、それぞれにおいて、PCもしくはDLNPCの少なくともどちらか一方がファンブルを引く確率を求め、最後に求めた値を足し合わせればよい。算出式は次の通りである。
$$ \frac{1}{2}((p_{\beta}-p_{\gamma})-p_{\delta})(p_{fc}+p_{fd})+p_{\delta}(1-(1-p_{fc})(1-p_{fd})) $$
　
$p_{\alpha}=0.78111257,p_{\beta}=0.21888743,p_{\gamma}=0.00150232,p_{\delta}=0.05030506$、そして$p_{fc}=p_{fd}=0.25$を仮定するとロスト確率は$0.06377847625$、およそ6.4%と求まる。

シミュレーションとの整合

　生存確率を求めたときと同様の方法でシミュレーション値との整合を確認する。100万回のシミュレーション結果である$\bar{X}$は$0.063917$である。このとき$\epsilon=0.002$としてシミュレーション結果が$p±\epsilon$に収まる確率は99.9%と求まる。$\bar{X}$はこの範囲内に収まっており、シミュレーション結果と理論値に矛盾はないことが分かる。

      import math

import numpy as np
from scipy.stats import norm


p=0.06377847625
sigma=p*(1-p)
barX=0.063917
n=1000000
epsilon=0.0002

z_lower=(math.sqrt(n)*((p-epsilon)-p))/sigma
z_upper=(math.sqrt(n)*((p+epsilon)-p))/sigma

lower=norm.cdf(z_lower)
upper=norm.cdf(z_upper)

print(upper - lower)
    

参考文献

https://kotobank.jp/word/%E7%AD%89%E6%AF%94%E7%B4%9A%E6%95%B0-104058#E7.99.BE.E7.A7.91.E4.BA.8B.E5.85.B8.E3.83.9E.E3.82.A4.E3.83.9A.E3.83.87.E3.82.A3.E3.82.A2