【MHW:IB】装飾品の確率論

装飾品コンプの確率分布の最頻や平均(期待値)を計算していきます．

前提条件
報酬期待値
排出平均確率
多項分布
マルコフ連鎖
- モンテカルロ法

前提条件

滅日を周回するとして、封じられた珠、刻まれた珠、太古の珠の3種類から4スロットを0個の状態からコンプリートするということを考えます．
滅日を $m$ 周し、珠が累計 $\{a_r m \}=\{a_1 m,a_2 m,a_3 m\}$ 個でたときにコンプリートしたとします． $a_r$ は1周ででる珠 $r$ の個数期待値とします．
つまり、 $m$ 周すると、珠が $n=(a_1+a_2+a_3)m$ 個数得られます．

また、挑戦・回避のような装飾品は5個でコンプとするものを弱コンプ、7個でコンプとするものを強コンプと定義します．

報酬期待値

1周したときの報酬の期待値をまず計算していきます．滅日の珠に限っていえば、報酬内容確定確定枠数が1つ(封じられた珠3個)、報酬確定枠数が3つ、抽選確率22/32で報酬抽選枠が最大4つまで増えるようです．
つまり報酬枠の期待値というのは、

$1+3+(\frac{22}{32})+(\frac{22}{32})^2+(\frac{22}{32})^3+(\frac{22}{32})^4=5+\frac{46433}{65536}=5.70851...$

となっているようです．このうち、1枠分は封じられた珠3つで固定されていて、残りの4.70851...枠で各3種の珠が確率ででてきます．

封じられた珠

封じられた珠は1枠に対して、8%の確率で1つ、1%の確率で3つでるようです．すなわち、1周で得られる期待値は

$3\times 1\times 1 +(1\times 0.08 + 3\times 0.01)\times(4+\frac{46433}{65536})=3+\frac{3394347}{6553600}=3.51793...$

となります．

刻まれた珠

こちらは1枠に対して、41%の確率で1つでるようなので、

$1\times 0.41 \times(4+\frac{46433}{65536})=1+\frac{6098057}{6553600}=1.93049...$

となります．

太古の珠

こちらは1枠に対して、50%の確率で1つでるようで、

$1\times 0.50 \times(4+\frac{46433}{65536})=2+\frac{46433}{131072}=2.35426...$

となります．

したがって、

$a_1=3.51793...$
$a_2=1.93049...$
$a_3=2.35426...$

です．

排出平均確率

このように珠の種類が多く、確率もバラバラである場合計算が非常に煩雑になってしまうので、珠とそこからの排出確率を重みづけして平均化していきます．
珠 $r$ から装飾品 $j$ が排出される確率を $p_{j,r}$ と置きます．珠の排出期待値で重みづけ平均化し、それを装飾品 $j$ の排出確率 $p_j$ とします．

$p_j = \frac{a_1 p_{j,1} + a_2 p_{j,2} + a_3 p_{j,3}}{a_1 + a_2 + a_3}$

こうすることで、珠の種類数というのを除外して考えることができてモデルを単純化できます．(両辺に $a_1+a_2+a_3$ をかければ一周での獲得数期待値となる)

例えば火炎珠II【4】は封じられた珠から1.11%、刻まれた珠から0.75%、太古の珠から0.28%で排出されます．なので、平均化された仮想珠から、

$p_{\rm 火炎珠II【4】}=\frac{3.51793\times 0.0111 + 1.93049\times 0.0075 + 2.35426\times 0.0028}{3.51793 + 1.93049 + 2.35426}$
$=0.7705\%$

の確率で排出されるとみなせます．一周での平均化仮想珠の排出数期待値は $a_1+a_2+a_3=7.80268$ 個なので、一周あたりの火炎珠II獲得数期待値は $7.80268 \times 0.007705=0.0550$ 個となります．

多項分布

装飾品は $J$ 種類あり、ラストに揃う装飾品を $j$ とラベリングします． $n-1$ 個目の時点でラストの装飾品 $j$ 以外(装飾品 $k$ )はコンプしていて $j$ がラスト1個手前までコンプしている確率を $P^{\rm qc}_j(n)$ (quasi complete)とおき、 $n$ 個目でコンプするということを考えます．

装飾品 $j$ が $n$ 個目でコンプする確率 $P^{\rm com}_j (n)$ および最終的に求めたい装飾品が $n$ 個数目ではじめてコンプする確率 $P^{\rm com} (n)$ は、

$P^{\rm com}_j(n)=P^{\rm qc}_j(n) p_j$

$\displaystyle P^{\rm com}(n)=\sum_j P^{\rm com}_j(n)$

と表せます． $n$ に関する無限和をとれば全事象を表すので、 $\displaystyle c := c_j + \sum_{k} c_{k}$ とおいて、

$\displaystyle \sum_{n=c}^{\infty} P^{\rm com}(n) =1$

を満たします．

装飾品 $j$ 、 $k$ のコンプ個数を $c_j$ 、 $c_{k}$ とおくと、 $P^{\rm qc}_j(n)$ というのは、ラスト一種類の $j$ は $c_j -1$ 個ちょうど、それ以外の $k$ は少なくとも $c_k$ 個でて、残りの枠は $j$ が外れる確率といえます．
この少なくとも...という条件が組み合わせ数を爆発的に増やしています．これを計算するには余事象を考えると組み合わせ数はぐっと抑えることができるのですが、のちでやりたいことに合わせて、愚直に定式化していきます．

自然数(0を含む) $x_j \ (\leqq n - c)$ に対して、その順序付きの分割(結合) $X_j:=(x_1,x_2,\cdots ,x_{k \not = j} ,\cdots x_J)$ を考えます．成分 $x_k$ は $0$ でもよいとし、 $0$ 成分も含めた項数は $J-1$ となります．
$x_k$ は各装飾品 $k$ がコンプ個数から余剰に出た個数を表していて、それぞれ $c_{k} + x_{k}$ 個でたということです． $x_j$ は余剰個数の総和 $\displaystyle x_j :=\sum_{k} x_{k} \leqq n - c$ を表しています．
この順序付き分割全体を $C_n (x_j)$ とおきます．

求める確率というのは片端固定版の多項分布の総組み合わせの和です．
まず $\displaystyle p:=\sum_j p_j$ 、 $q:=1-p$ とおいておきます． $q$ というのが1~3スロ系の装飾品が出る確率に対応しています．そして、次の関数を定義します．

$\displaystyle f_{j}(X_j) = \prod_{k\not= j} \frac{ p_{k}^{c_k+x_{k}} }{ (c_k+x_{k})! }$

そうすると、

$\displaystyle P^{\rm qc}_j(n) = \frac{ (n-1)! }{ (c_j-1)! } p_j^{c_j-1} \sum_{x_j =0}^{n-c} \frac{q^{n-c-x_j}}{(n-c-x_j)!} \sum_{X_j \in C_n(x_j)} f_{j}(X_j)$

と書き表すことができます．

正規分布近似

ここで、 $n-1$ 個の珠を引いたとき、 $j$ が $c_j-1$ 個ちょうどでる確率 $P^{\rm sq}_j(n)$ (single quasi complete)を計算してみます．これはとても簡単で、

$\displaystyle P^{\rm sq}_j(n) = \frac{(n-1)!}{(c_j-1)!(n-c_j)!}p_j^{c_j-1}(1-p_j)^{n-c_j}$

とかけます．この単純なものから多項分布の性質を見出してみます．

二項分布である $P^{\rm sq}_j(n)$ は、 $n$ が十分大きければ正規分布に近似できるので、

$\displaystyle P^{\rm sq}_j(n) = \frac{1}{\sqrt{2\pi(n-1)p_j(1-p_j)}}\exp{\left(-\frac{(c_j-1-(n-1)p_j)^2}{2(n-1)p_j(1-p_j)}\right)}$
$\displaystyle \approx \frac{1}{\sqrt{2\pi n p_j(1-p_j)}}\exp{\left(-\frac{(c_j-np_j-1)^2}{2np_j(1-p_j)}\right)}$

として扱えます．

対象となる装飾品が2種類の場合(2次元多項分布)の簡易版のグラフをdesmosで描きました．パラメータを動かすことで確率の変動が視覚的にわかりやすいです．
www.desmos.com

最長時間近似

※最長時間近似という言葉は造語です．
上のdesmosでいろいろパラメータをいじると見えてくるのが、 $p_1 < < p_2$ だったり、 $c_1 >> c_2$ だったりすると、最終的な確率分布に2の方は殆ど寄与せず、単純な二項分布として近似できるということがわかります．すなわち、 $n$ が十分大きいとき、相対的に排出確率が高く、コンプ個数が少ないものは外れ枠に圧縮し、考えるべき変数を減らして近似できるというわけです．2種類以上ある場合には、各一種類のみに着目したときのコンプまでの時間期待値 $T^{\rm com}_j:=E[ n ]$ で比較するのがよさそうです．

$\displaystyle T^{\rm com}_j=\sum_{n=c_j}^{\infty} n P^{\rm com}_j(n)$
$\displaystyle =\sum_{n=c_j}^{\infty}\frac{n!}{(c_j-1)!(n-c_j)!}p_j^{c_j}(1-p_j)^{n-c_j}$
$\displaystyle =\frac{c_j}{p_j}$

近似の仕方としては、まずコンプ時間期待値の最大値を求めます．その最大となる装飾品を $j$ とラベリングします．ある $\delta$ を定めて、 $\displaystyle \frac{c_k p_j}{p_k c_j} \geqq \delta$ を満たす装飾品 $k$ を準最長時間としてグルーピングします．逆に、それを満たさないもの(一時的にそれらを $l$ とラベリングします)は外れ枠として圧縮されます．すなわち、 $l$ は1～3スロ系のものと同一視して $\displaystyle q\rightarrow q+\sum_l p_l$ と再定義されます．

$\delta$ を小さくすればするほど近似精度が上がりますが準最長時間グループに分類される種類数が増えて計算コストが上がります．今回の装飾品の例では $\delta=0.5$ で十分ぽそうです．

一般化クーポンコレクター問題

最長・準最長時間に分類される種類数とそれらのコンプ個数が極めて小さければ、後述する遷移確率行列を考えて計算するのも現実的ですが例えば強コンプの場合で具体的に計算してみると、最長時間が5329.13となり、それらは挑戦・回避、挑戦・体術、全開・回避、全開・体術、底力・回避、底力・体術の6つです． $\delta = 0.5$ として準最長時間に分類されるものは(最長も含め)102種類となります．準最長時間を無視して( $\delta=0.8$ 程度に相当)最長時間のみで計算しようとしても、たったの6種類だけで状態数は $(7+1)^6 = 262144$ 個となります．
弱コンプの場合でも、最長時間が3806.52で、挑戦・回避などの22種類、 $\delta = 0.5$ の準最長まで含めると114種類あります．最長時間だけで状態数 $(5+1)^{22}=131621703842267136$ 個です．とはいえ、最長時間だけをみれば対称性があることが多いので工夫すれば計算できる可能性はあります．

等確率・等コンプ個数のものが複数種類あるだけのような状況であれば、それは一般化クーポンコレクター問題(GCCP)のSSC(the Single Slippage Configuration)モデルというものに相当します．
これはある程度研究されている内容ですが、確率分布を計算できているものは見つかりませんでした．

GCCPのSSCを考えるものとして、「当たり」が $m$ 種類あって、排出確率がそれぞれ $p$ 、コンプ個数が $c$ とします．そうすると、コンプ確率というのは、

$x\ (\leqq n-mc)$ の0を含む順序付き分割 $X:=(x_1,x_2 ,\cdots ,x_{m-1})$ とその全体 $C_n(x)$ 、外れ確率 $q:=1-mp$ をおいて、

$\displaystyle f(X) = \frac{ p^{(m-1)c+x} }{ (c+x_1) ! \cdots (c+x_{m-1})! }$

とすると、

$\displaystyle P^{\rm com}(n) = \frac{ m(n-1)! }{ (c-1)! } p^c \sum_{x =0}^{n-mc} \frac{q^{n-mc-x}}{(n-mc-x)!} \sum_{X \in C_n(x)} f(X)$
$\displaystyle = \frac{ m(n-1)! }{ (c-1)! } p^{mc} \sum_{x =0}^{n-mc} \frac{p^{x} q^{n-mc-x}}{(n-mc-x)!} \sum_{X \in C_n(x)} \frac{ 1 }{ (c+x_1) ! \cdots (c+x_{m-1})! }$

と書き表すことができます．

最後の式で最も大変なのは

$\displaystyle \sum_{X \in C_n(x)} \frac{ 1 }{ (c+x_1) ! \cdots (c+x_{m-1})! }$

の部分です．簡単のため、 $x$ は $m-1$ の倍数、つまり $x=(m-1)s$ であるとき(等分配)を考えてみます．
すべて $x_k=s$ のとき、分母の部分は最小となるはずです．

$\displaystyle (c+s)! \cdots (c+s)! = (m-1)(c+s)!$

ここから次に大きいのは1個分ずらしたもので、

$\displaystyle (c+s+1)! (c+s-1)! (c+s)! \cdots (c+s)! = \frac{c+s+1}{c+s}(m-1)(c+s)!$

これは $m-1$ 通りあります．

次は2個分ずらすものが考えられます．

$\displaystyle (c+s+2)! (c+s-2)! (c+s)! \cdots (c+s)! = \frac{(c+s+2)(c+s+1)}{(c+s)(c+s-1)}(m-1)(c+s)!$

$\displaystyle (c+s+2)! (c+s-1)! (c+s-1)! (c+s)! \cdots (c+s)! = \frac{(c+s+2)(c+s+1)}{(c+s)^2}(m-1)(c+s)!$

これらは $\displaystyle {}_{m-1}{\rm C}_{2}\times (m-3)$ 通り、

$\displaystyle (c+s+1)! (c+s+1)! (c+s-2)! (c+s)! \cdots (c+s)! = \frac{(c+s+1)^2}{(c+s)(c+s-1)}(m-1)(c+s)!$

$\displaystyle (c+s+1)! (c+s+1)! (c+s-1)!(c+s-1)! (c+s)! \cdots (c+s)! = \frac{(c+s+1)^2}{(c+s)^2}(m-1)(c+s)!$

これらは $\displaystyle {}_{m-1}{\rm C}_{2}\times {}_{m-3}{\rm C}_{2}$ 通りあります．

これを計算していくのは骨が折れるので、 $1/(m-1)(c+s)!$ に、項数 $m-1$ で順序付きの非負整数の分割数をかけたものを近似値とするとよさそうです．
$x$ を $0$ 成分を含んで $m-1$ 個の項数に分割する方法は、 $x$ 個のものと $(m-1)-1$ 個の仕切りを一列に並べる重複組み合わせ、すなわち $x+m-2$ 個のものから $m-2$ 個を選ぶ組み合わせなので、 $\displaystyle {}_{x+m-2}{\rm C}_{m-2}=\frac{(x+m-2)!}{x!(m-2)!}$ 通りです．したがって、

$\displaystyle P^{\rm com}(n) \approx \frac{ m(n-1)! }{(m-1)! (c-1)!} p^{mc} \sum_{x =0}^{n-mc} \frac{p^{x} q^{n-mc-x}(x+m-2)!}{x!(n-mc-x)!\Gamma(c+\frac{x}{m-1}+1)}$

と近似することにします．対数つまり情報量に直しておきます．

$\displaystyle I^{\rm com}(n) = -\log P^{\rm com}(n)$
$\displaystyle = -\log m -\log (n-1)! +\log (m-1)! +\log (c-1)! -mc\log p -\log \sum_{x =0}^{n-mc} \frac{p^{x} q^{n-mc-x}(x+m-2)!}{x!(n-mc-x)!\Gamma(c+\frac{x}{m-1}+1)}$

スターリングの近似式を使って、

$\displaystyle I^{\rm com}(n) = -\log m -(n-1)\log (n-1) +(m-1)\log (m-1) +(c-1)\log (c-1) -mc\log p + n - m - c -\log \sum_{x =0}^{n-mc} \frac{p^{x} q^{n-mc-x}(x+m-2)!}{x!(n-mc-x)!\Gamma(c+\frac{x}{m-1}+1)}$

しかし、最後の項はどうしようもありません．なにかいい現実的にできる計算方法があればいいのですが．

逐次合成

先ほどのdesmosをヒントに、1種のものから同種のものが2つあるものへの合成(変換)ができると、それを逐次的に繰り返すことで2の累乗数の確率分布が得られそうです．
www.desmos.com

この青い点線のもの(以下 $P_1(n)$ )から赤い実線のグラフ( $P_2(n)$ )を $P_1$ の形で近似的でもいいので作ることが目標です．

$\displaystyle P^{\rm com}_1(n) = \sqrt{\frac{p}{2\pi n (1-p)}} \exp\left(-\frac{(c-np-1)^2}{2np(1-p)}\right)$

$\displaystyle P^{\rm com}_2(n) = \frac{1}{\pi n\sqrt{1-2p}} \sum_{x=0}^{\infty} \exp\left(-\frac{1-p}{2np(1-2p)} \left( 2(1-p)(c-np+x)(c-np-1)+(1+x)^2 \right) \right)$

それぞれの極大値や幅を求めて比較というのができると良かったのですが、 $P_2(n)$ のほうがあまりにも厳しすぎたのでdesmosで性質を探ってみると、おおよそ、青のコンプの時間期待値と赤の極大値が一致しています．完全に一致しているというわけではなく、同じようなずれ方をしているということで、オーダーまでずれるということはなさそうです．
数学としてはこんなやり方ではだめですが、先に進めずに結果が得られないというのも面白くないのでとりあえずこれを利用してみます．
$P_1(n)$ は極大となる $n$ のほうは(微分などで)簡単に求められて、

$\displaystyle n=\frac{1}{2}\left(1-\frac{1}{p}\right) + \sqrt{\frac{1}{4}\left(1-\frac{1}{p}\right)^2 + \left(\frac{c-1}{p}\right)^2}$
$\displaystyle \approx \frac{\frac{1}{2}(\sqrt{4c^2-8c+5}+1)}{p}$

です． $p$ は小さいので-1次近似しておきました． $P_2$ が $P_1$ と同じ関数形で近似できているとすると、ここの点が $P_1$ のほうのコンプ期待値 $\frac{c_1}{p_1}$ の点と近しいというわけです．
関数をみると、 $p$ のほうをずらすと幅が変わってしまうので、グラフ的に幅自体はそこまで変わってないことをみると、 $c$ のほうを書き直せるとよさそうなので、逐次的に

$c_2=\sqrt{c_1(c_1+1)} + 1$

としていくと近しい関数の概形が得られることが期待できます．
最長時間の装飾品が $2^m$ 種あったとき、上述のように $c$ を逐次的に計算し、 $m$ 回目の値 $c_m$ による確率分布を近似値として最終的に用いていきます．
つまり、関数

$f(c):=\sqrt{c(c+1)} + 1$

を定め、 $m$ 回反復合成関数 $f^m(c)$ によって、

$\displaystyle P^{\rm com}(n) = \sqrt{\frac{p}{2\pi n (1-p)}} \exp\left(-\frac{(f^m(c)-np-1)^2}{2np(1-p)}\right)$

とかけます．

確率分布の概形

弱コンプで考えていきます．データを入力し、装飾品の種類数は $293$ 種で、全コンプ個数は $871$ 個(つまりコンプ時間の理論値)、単一の最長時間は排出平均確率 $0.131 \%$ でコンプ個数5個の $3806.5$ 個目で、このようなものは $22$ 種類ありました．
$f^4(5)=10.9$ 、 $f^5(5)=12.4$ で、 $22$ は $2^4=16$ と $2^5=32$ の間なので $11$ としておきます．
コンプの期待値は $8374.3$ 個、すなわち $1073$ 周、最頻値は $7242.4$ 個の $928$ 周です．様々な雑な近似により実際は結構ずれると思いますが、ざっくり1000周前後でコンプリートするということがわかりました．
www.desmos.com

言い換えれば「平均的には」、約1000周すれば22種類ある攻撃IIと同レアのものうち、どれか1つは11個ほど揃うくらいになっていて、なおかつ滅日から排出される4スロット装飾品はすべて揃っているということになります．戦闘時間2分、帰還まで1分、ロードや出発準備で1分とすると、67時間はかかります．
上位1%くらいの人は約500周までにコンプ可能で、下位1%くらいの人はコンプするのに約2000周以上を要します．※あくまで滅日のみから得られるものを考慮した場合で、他クエストに行っている場合は当然その分緩和されます．

計算に使用したcsvファイルはこちら(2のほうが弱コンプです)．入力はMinaiさんに手伝ってもらいました．
https://drive.google.com/drive/folders/1LmCwuds7NVm_p7rVB9BVp0cs8SXe6oRJ?usp=drive_link

実は強コンプのほうで同じ計算をしてみると、こちらのほうが早く揃いやすいという結果になってしまっています．これは最長時間のみを考慮した場合、強コンプはたったの7種類しかなく、準最長の寄与がわりとあるのにも関わらず無視してしまっているからだと思われます(つまり $\delta=1$ としているのが問題)．

以下おまけ

マルコフ連鎖

別の考え方として、各装飾品の所持状況(各コンプ個数が最大値)を確率過程とするとマルコフ連鎖とみなせます．
全装飾品を設定したコンプリート個数に到達した状態が吸収状態となっています．
装飾品 $j\ (1,2,\cdots , J)$ のコンプ個数を $c_j$ とすると、全状態の総数というのは所持数0個ということも考慮すると、 $N = (c_1+1)(c_2+1)\cdots (c_J +1)$ 通りになります．
$N\times N$ の遷移確率行列 $P$ を考えることで、珠 $n$ 個目で初期状態から吸収状態 $N$ になる確率を計算できます．

計算できます、とはいいましたが293種類の装飾品で一般的なコンプ個数を設定した場合、状態数は $9.58\times 10^{210}$ 通りにもなるので $(9.58\times 10^{210})\times (9.58\times 10^{210})$ サイズの行列を計算する必要になり、人間はもちろんコンピュータでも不可能でしょう．
そこで先ほどの最長時間近似が効いてきて、この再定義によって種類数が2~3種類程度に圧縮されればぐっと行列サイズを下げて現実的に計算することができます．
しかし、5種類以上にもなってくると行列サイズが数千以上になっていき計算時間がとても長くなるか、強制終了してしまいます．

とはいえ遷移確率行列の形だけでも記していこうとは思います．
まず、 $\displaystyle i=\sum_j x_j$ を満たす状態 $X_i=(x_1,x_2,\cdots,x_J)$ を考えます(さっきまでとは異なる定義であることに注意してください)．
同じ $i$ でも内容が異なる状態 $X'_i$ への遷移する確率はいずれかの装飾品を減らす必要があって、そのような確率は0(装飾品が勝手に消えることはない)です．
もちろん、 $i-1$ 以下のいかなる状態への遷移も装飾品が減ることがないことから0です．
一方、状態がかわらない、つまり4スロ系の装飾品の所持状況が変わらないことはあり、その自己遷移の確率はコンプ個数に達している装飾品の排出確率の和

$\displaystyle q(X_i) := \sum_{\substack{j \\ x_j=c_j}} p_j$

を用いて、 $q+q(X_i)$ と書けます．そして、 $X_i$ から装飾品 $j\ (x_j < c_j)$ がひとつ増えた状態 $X_{i+1}^{j}=(x_1,x_2,\cdots,x_j+1,\cdots,x_J)$ への遷移確率というのが $p_j$ なわけです．
遷移確率行列の $(X_i,X_i)$ 成分は $q+q(X_i)$ 、 $(X_i,X_{i+1}^{j})$ 成分は $p_j$ というわけです．
ある行を取り出すとこんな感じです(行和は確率の和なので1になる)．

コンプした状態 $(c_1,c_2,\cdots,c_J)$ は吸収状態なので自己遷移確率が $1$ (ここを $(N,N)$ 成分とします)で、他成分は0です．
このような行を $N$ 個並べることで遷移確率行列 $P$ ができます．その遷移確率行列を $n$ 乗して初期状態の行の最後の列 $N$ の値が珠 $n$ 個目で初期状態から吸収状態になる確率 ${P^n}_{0N}$ になります．