スクフェスの特技において「パターンC」という仕様が最近発見されました。これは、浦の星女学院パソコン部さんによって存在が提起されたスクフェスの仕様で、「○○秒間発動する」という内容の特技において、特技発動中に同じ特技の発動機会が2回以上あった場合、1回の発動機会しか考慮されない というものです。詳しくは、浦の星女学院パソコン部さんの以下のツイートにて解説されています。

マニアックですが超重要なネタです。高レベルのCOMBO FEVERを使うと画像のような状況が起きます。この場合、最も確からしい挙動はA〜Cのどれだと思いますか？実機のスコアと比較した結果、最も近いのはCでした。SIF Matchは現在Bなので、近い内にCに変更しようと思います。Aで計算すると壊れます😇 pic.twitter.com/oIL0mTvIbJ
— 浦の星女学院パソコン部 (@devsifmatch) 2019年4月22日

この現象は多くの方によって検証がなされており、例えば浦の星女学院パソコン部さんは実機のスコアと比較することによって検証を行なっていたり、穂夢さんとおくーやさんは特技の連続発動数をカウントすることによって検証 *1 *2 を行なっています。

この現象は、曲全体に対する特技発動時間の割合、つまりカバー率に影響を及ぼします。この時のカバー率を、浦の星女学院パソコン部さんは計算によって、 *3 おくーやさんはモンテカルロ法によって *4 求めています。

この記事の内容

この記事では、パターンCにおいて 何がどれくらい、なぜカバー率の低下に寄与しているのか を知ることを目的とします。そのために、モンテカルロ法によって求められた正確なカバー率の値をよく近似し、なおかつなるべく簡単な形式の理論式を導出します。具体的には、以下のことを行います。

マルコフ過程を用いて「スキップされる発動機会」と「有効な発動機会」などを表現し、系が定常状態にあるという近似を用いて、部員ステータスから簡単にカバー率を求めることのできる近似的な理論式を導きます。
結論として、パターンCの効果は、平均的な発動確率が部員のパラメータに依存した係数倍に縮められることに相当する ことを導きます。
このことから、どのような部員がパターンCによって強い補正を受けるか を分析します。

最終的に、曲のノーツ数と長さを $N_{\rm notes},\ T_{\rm music}$ 、特技が $N_{\rm trig}$ ノーツおきに $p$ の確率で $T_{\rm eff}$ 秒持続するとき、その曲全体における カバー率の期待値 $\hat r_{\rm cover}$ の近似式 として

$\displaystyle \hat r_{\rm cover} \approx \frac {p} {1 + (m-1)p + p^2q} \frac {N_{\rm maxtrigs} T_{\rm eff}} {T_{\rm music}}$

という式を導出します。ここで、 $N_{\rm maxtrigs}:= \left\lfloor \frac {N_{\rm notes}} {N_{\rm trig}} \right\rfloor$ は曲中の特技の発動機会の回数、 $m,q$ はそれぞれ特技発動に要する時間を $T_{\rm trig}$ とした時の $\frac {T_{\rm eff}} {T_{\rm trig}}$ の 整数部分および小数部分、 $\left\lfloor \cdot \right\rfloor$ は床関数を表しています。

2019.5.5 追記： 本記事の理論式を少しだけ改良したバージョンを作りました。改良された理論式では、予測誤差が4.0ポイントから1.6ポイント程度に改善されています。改良版の理論式の導出は、この記事で紹介している理論式の延長にあるため、導出について詳しく知りたい方は、現在の記事もご覧ください。

snowharinchan.hatenablog.com

2019.5.12 追記： 下記の記事の方法によって、上記の近似式をより簡単に導出することができます。

snowharinchan.hatenablog.com

この理論式はどれほど現実に当てはまっているのでしょうか。これを確かめるため、モンテカルロ法で求めたカバー率を真値とした時の、理論式を直接適用した時の場合の予測結果と、補正係数を掛けた場合の予測結果が以下の図になります。前者は予測誤差（RMSEで評価）が0.0401程度、後者は予測誤差が0.0210程度になりました。つまり、それぞれ カバー率を4ポイント、2ポイント程度の誤差で言い当てられている ことになります。

f:id:snowharinchan:20190503135530p:plain

ここで、真値のデータとして、おくーやさんのブログ記事中の6枚の画像に掲載されているデータを用いています。このような解析を行なって頂き、ありがとうございます。画像中の数値は、OCRによって文字起こしすることで利用しています。

また、この理論が面白いのは、 $\frac {T_{\rm eff}} {T_{\rm trig}}$ の整数部分および小数部分がカバー率計算の上で大きな役割を果たしている ことです。整数部分はマルコフ過程の状態数、小数部分はマルコフ過程の遷移確率に影響しています。同じ量から計算される整数部分と小数部分が、このようにモデル中の別々の量を定めている事も、注目すべき部分です。

また、近似的な理論式でありながら、比較的小さい予測誤差のモデルが導出されていることから、このモデルはこのパターンCという現象を説明する上である程度概要を捉えたモデルになっていることがわかります。

また、上記は近似的な理論でしたが、この記事の最後にて、モンテカルロ法を用いて計算された真値と同様に精密な期待値の理論値を求める方針 について述べます。

以下では、この理論式を導出する過程を説明し、最後に理論の検証を行います。

導出

以下では、この理論式の導出方法について説明します。

マルコフ過程によるパターンCのモデル化

マルコフ過程を用いてパターンCをモデル化する際に最も重要となるのは、発動機会がスキップされるのがいつであるかを求めることです。この節では、まずはじめにその点について説明します。

発動機会がスキップされるのはいつか

発動機会が6秒ごと、発動時間が8.5秒の特技を考えます。これは例えば、500ノーツ120秒の曲（ノーツは均等に分布しているとします）において、25音ごとに発動する特技が120*25/500=6となり、これに該当します。

この特技が3回連続で発動した様子を表したのが下図です。

f:id:snowharinchan:20190503135533p:plain:w400

ここで、3回目の特技の発動時のみ、特技の発動時間中に2度発動機会が訪れてしまっているため、無効な発動機会が一度現れています。このような機会を無効発動機会と呼ぶことにします。

ここで上記の図から、特技発動時間間隔 $T_{\rm trig}$ と発動持続時間 $T_{\rm eff}$ の比 $r_{\rm eff,trig} := T_{\rm eff} / T_{\rm trig}$ の 小数部分 が主な原因になって無効発動機会が現れていることがわかります。この小数部分を $q:=r_{\rm eff,trig} - \left\lfloor r_{\rm eff,trig} \right\rfloor$ とおきます。例えば今 $T_{\rm eff} = 8.5\,{\rm sec},\ T_{\rm trig} = 6\,{\rm sec}$ のため $r_{\rm eff,trig} = 8.5/6 \approx 1.42$ となりますが、 $1.42$ を0倍、1倍、2倍、3倍、4倍、…としていくと $0,\ 1.43,\ 2.84,\ 4.26,\ 5.68,\cdots$ となり、整数部分のみに着目すると $0,\ 1,\ 2,\ 4,\ 5,\cdots$ となり、3がスキップされていることがわかります。つまり、発動持続時間 $T_{\rm eff}$ を0から順に整数倍していった時の整数部分に着目し、スキップされる整数が無効発動機会になっていることがわかります。ここで、発動機会は0番目から数えることに注意します。

発動開始時刻がランダムであると仮定した場合

さて、前節では、無効発動機会のパターンの鍵は $q$ が握っていることがわかりました。ノーツが均等に分布しているような曲では、特技が連続で発動した時の無効発動機会の現れるパターンは必ず一定です。しかし一定と言えどこのパターンは複雑で、先ほどのように $q=0.42$ の場合は、スキップされる数は $3,5,8,10,\cdots$ のように、この部分のみ切り出して見ると一見不規則になります。

ここで、このパターンを正確に表現するのではなく、近似的に表現するため、 無効発動機会はランダムに現れる とします。

また 簡単のため、以下ではしばらく発動持続時間中に最大でも2回のみ発動機会を含むような特技のみを考える ことにします。この仮定はやがて外し、後に一般的な理論を構築します。また、無効発動機会について考慮するため、特技発動時間間隔 $T_{\rm trig}$ が発動持続時間 $T_{\rm eff}$ 以下であるような特技を考えます。これらを合わせると、 $r_{\rm eff,trig} := T_{\rm eff} / T_{\rm trig}$ について、 $1 \leq r_{\rm eff,trig} \lt 2$ となります。

そのために、次のような状況を考えます。

数直線状の区間 $0 \leq x \leq 1$ 上の一様分布からランダムに一つ点 $x$ をサンプルする。 $x$ を端点とし、 $x$ より大きい点をもう一つの端点とする長さ $1+q$ の線分を考える。この線分が数直線上の整数点を2つ以上含む確率は $q$ で与えられる。

上記において、線分を特技の発動持続時間、数直線上の整数点を特技の発動機会と考えると、特技の開始時刻が完全にランダムであるとき、特技発動時間中に無効発動機会を含む確率は $q$ で与えられる ことがわかります。つまり $q$ が大きければ大きいほど長さ2の線分に近付くことで無効発動機会が生まれる確率が増え、逆に $q$ が小さいほど長さ1の線分に近付き無効発動機会が生まれる確率が減ります。

元のモデルでは、無効発動機会は小数部分 $q$ が整数倍されることで整数に繰り上がるタイミングで現れましたが、 $q$ が大きいほど繰り上がりが頻繁に起き、逆に $q$ が0に近ければ近いほど、整数がスキップされるまでに特技が連続で発動しなければならない回数が増えていきます。したがって、この元のモデルの状況と、開始時刻をランダムとした状況は近いことがわかります。

マルコフ過程による表現

上記の知見に基づいて、スクフェスの一つの楽曲中にわたる特技発動の様子をモデル化した、以下の図のようなマルコフ過程を考えます。

f:id:snowharinchan:20190503135539p:plain:h300

このマルコフ過程は以下のようにして構成されています。

Off, On, Skipped の3つの状態を持つ。
初期状態はOffの状態で、曲の開始時に初期状態からスタートする。
発動機会が一度訪れるごとに時刻を1つ進める。この時、この発動機会が有効発動機会であっても無効発動機会であっても、いずれの場合においても時刻を1進める。（つまり、前節の図において、真ん中の数直線の区切りの全てのタイミングにおいて時刻を1進めることとする。）
Offは特技が発動していない状態、Onは特技が発動している状態、Skippedは特技が発動しているが無効発動機会が訪れた状態を表すことにする。
On から On への遷移確率を $p(1-q)$ 、On から Skipped への遷移確率を $pq$ とする。 一方、Off から On への遷移確率は $p$ とする。これは、Off の状態から On になる際は（ $1 \leq r_{\rm eff,trig} \lt 2$ と仮定したことにより）無効発動機会が現れることはないが、特技が発動している状態で次の特技が予約された時、その開始時刻はランダムに定まるという近似を行い、したがって確率 $q$ で特技発動時間中に無効発動機会を含む様子を表している。
Skipped から On への遷移確率を $1$ とする。
Off 状態への遷移確率はSkippedを除くいずれの状態からも $1-p$ である。

このマルコフ過程には二つの特徴があります。

Skipped 状態を作ることで、余分に状態が消費される状況をモデル化している
On 状態を訪問した回数を数えることで、特技が発動した回数を数えることができる

以下では後者について説明します。

On 状態を訪問した回数を数えることで、特技が発動した回数を数えることができる

さて、私達が興味があるのはカバー率の期待値なので、楽曲中に特技が発動していた時間の期待値を知るのが目的になります。同じ特技が重複して発動することはないので、特技の発動回数の期待値を求めて $T_{\rm eff}$ 倍すれば良いことになります。（ここで、楽曲終了付近で発動した特技の端数時間はここでは考えません。）

ここで、上記のマルコフ過程において、「特技が発動した」という事を表すのは、太線で表した辺を通った際になります。 これらの辺は、それぞれ

A: Off → On: 元々特技が発動していない状態から、特技が発動した。
B: On → On: 特技発動中に、次の特技の発動を予約した。この時、次回の特技発動時は無効発動機会を含まないように予約した。
C: Skipped → On: この辺を通るからには、必ず On → Skipped の辺を通っている。 On → Skipped の辺は、「特技発動中に、次の特技の発動を予約した。この時、次回の特技発動時は無効発動機会を 含むように 予約した。」という状況を表している。

という状況を表しており、これらのいずれかの辺を通るということで、今回考えている特技発動の場合を全て言い尽くしていることになります。

さて、本節冒頭の図を見ると、前の時刻から現在時刻にかけてこれらA,B,Cの辺のいずれかを通ったということは、現在時刻において On 状態にいることの必要十分条件である ことがわかります。

以上のことから、曲中の特技の発動回数を数えるには、曲中にて On 状態を訪問した回数を求めれば良い ことがわかります。

状態訪問数の期待値

ここで、記号の整理と定義を行います。

$T_{\rm music}$ ：曲の長さ。
$N_{\rm notes}$ ：曲のノーツ数。
$T_{\rm eff}$ ：特技の発動持続時間。
$N_{\rm trig}$ ：特技の発動に必要なノーツ数。
$T_{\rm trig} := T_{\rm music}N_{\rm trig}/N_{\rm notes}$ ：特技の発動に必要な時間。
$N_{\rm maxtrigs} := \left\lfloor N_{\rm notes}/N_{\rm trig}\right\rfloor$ ：曲中の発動機会の総数。
$r_{\rm eff,trig} := T_{\rm eff} / T_{\rm trig}$ .
$r_{\rm eff,trig}$ の小数部分を $q:=r_{\rm eff,trig} - \left\lfloor r_{\rm eff,trig} \right\rfloor$ とおく。

前節のマルコフ過程の確率遷移行列 $A$ は以下の通りとなります。

$\displaystyle A = \begin{bmatrix} 1-p & p & 0 \\\ 1-p & p(1-q) & pq \\\ 0 & 1 & 0 \end{bmatrix}$

さて、楽曲開始時の初期状態は $x_0 := [1\ 0\ 0$ ]で与えられます。また、楽曲中には $N_{\rm maxtrigs}$ 回の発動機会が存在します。また、On 状態の状態ベクトルは $x_1 := [0\ 1\ 0$ ]で与えられます。

したがって、ある一つの時刻 $k$ $(0 \leq k \leq N_{\rm maxtrigs})$ に Off, On, Skipped 状態を訪問している確率をそれぞれ $P_k(s=0),\ P_k(s=1),\ P_k(s=2)$ と書くことにすると、時刻 $k$ において On 状態を訪問している確率は

$\displaystyle P_k(s=1) = x_0 A^k x_1^\top$

で与えられることから、全時刻にわたって On 状態を訪問する回数 $N_{\rm trigs}$ （特技発動に要するノーツ数 $N_{\rm trig}$ とは異なります）の期待値 $\hat N_{\rm trigs}$ は

$\displaystyle \begin{aligned} \hat N_{\rm trigs} &:= \mathbb{E}\left[N_{\rm trigs}\right] \\\ &= \sum_{k=0}^{N_{\rm maxtrigs}} P_k(s=1) \\\ &= \sum_{k=0}^{N_{\rm maxtrigs}} x_0 A^k x_1^\top \end{aligned}$

で与えられます。

定常状態による近似

さて、上記の和を計算しても良いのですが、計算を簡単化するために、更にこのマルコフ過程が楽曲開始時から常に定常状態であったという近似を行います。

定常状態について説明するために、 $(p,q)=(0.81,0.42)$ の場合の、各時刻において各状態にいる確率の推移を見ると、以下のようになります。

f:id:snowharinchan:20190503135521p:plain

ここで、時刻5以降はほぼ同じ確率になっていることがわかります。この収束先の確率を定常状態と言います。この例では楽曲は時刻20まで流れているので（上記のグラフでは更に長い時刻がプロットしてありますが）、全体のほぼ3/4は定常状態にいることがわかります。つまり、 $x_0 A^k x_1^\top$ という量を計算していくと、比較的早い段階で計算結果が以降のどのような $k$ に対してもほぼ常に同じ値になるということが起きます。

この時の一定である確率の値を ${\pi}_s$ とおくと、定常状態では

$\displaystyle \begin{bmatrix} P_k(s=0) \\ P_k(s=1) \\ P_k(s=2) \end{bmatrix}^\top \approx {\pi}_s$

と書けます。

そこで、楽曲開始後（初期状態を除く）からずっと定常状態にいた という近似を行うと、

$\displaystyle \begin{aligned} \hat N_{\rm trigs} &= \sum_{k=0}^{N_{\rm maxtrigs}} x_0 A^k x_1^\top \\\ &\approx x_0 x_1^\top + \sum_{k=1}^{N_{\rm maxtrigs}} \pi_s x_1^\top \\\ &= N_{\rm maxtrigs} \pi_s x_1^\top \\\ \end{aligned}$

のように表すことができます。ここで、和の変形には初期状態では On 状態にいなかったことを用いました。

さて、 ${\pi}_s$ の具体的な値を求めます。定常状態 ${\pi}_s$ では次の時刻における確率が前の時刻と変わらないことから、

$\displaystyle \pi_s A = \pi_s$

と書けます。これは $A$ の固有値1に対する左固有ベクトル（つまり $A^\top$ の右固有ベクトル）であるため、これを求める *5 と

$\displaystyle x_{\rm eig} = \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & 1 \end{bmatrix}$

となります。これを要素の総和が1となるように正規化することで

$\displaystyle \pi_s =\frac 1 {\left(\frac{1-p}{p^2 q} + \frac 1 {pq} + 1\right)} \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & 1 \end{bmatrix}$

となります。

更に、最終的に $\pi_s$ は $x_1$ 倍してしまうので、 $s=1$ の要素のみ計算すればよいことになります。よって最終的に

$\displaystyle \begin{aligned} \hat N_{\rm trigs} &\approx N_{\rm maxtrigs} \pi_s x_1^\top \\ &= N_{\rm maxtrigs} \frac 1 {\left(\frac{1-p}{p^2 q} + \frac 1 {pq} + 1\right)} \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & 1 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \\ &= \frac {p} {1 + p^2q} N_{\rm maxtrigs}\\ \end{aligned}$

が得られます。

これに特技持続時間をかけて楽曲の長さで割ったのがカバー率なので、

$\displaystyle \begin{aligned} \hat r_{\rm cover} &:= {\mathbb E} \left[ r_{\rm cover} \right] \\\ &\approx \frac{\hat N_{\rm trigs} T_{\rm eff}}{T_{\rm music}}\\\ &\approx \frac {p} {1 + p^2q} \frac {N_{\rm maxtrigs} T_{\rm eff}} {T_{\rm music}} \end{aligned}$

のように、カバー率を部員や楽曲のパラメータによって表すことができました。

整数部分が2より大きい場合

現時点での理論式を用いて、実際の値の予測を行うと以下のようになります。

f:id:snowharinchan:20190503135523p:plain

ここで、図中の Integer part は $r_{\rm eff,trig}$ の整数部分を表しています。このように、整数部分が1の場合は説明できているのですが、整数部分が2、3の場合は非常にずれが大きくなっています。

そこで、この場合の理論を考えてみましょう。以下では、

$m:= \left\lfloor r_{\rm eff,trig} \right\rfloor$ ： $r_{\rm eff,trig}$ の整数部分

と定義します。

まず、 $m=2$ である時を考えます。先ほどの数直線の上の線分の図において、数直線上で長さ $2+q$ の線分を動かすことを考えると、この場合 発動機会の消費数は2回か3回に限られる ことがわかります。また、Off状態から特技が発動すると、必ず無効発動機会を1回消費して特技を発動させることがわかり、更に特技を予約で発動させた場合、確率 $q$ で更にもう1回余分に無効発動機会を消費することがわかります。

以上のことを反映すると、On状態の前に一回無効発動機会を挟む などを行うと良いことがわかります。そこで、以下のようなマルコフ過程を考えます。

f:id:snowharinchan:20190503135542p:plain:h400

$m=1$ の時との差分は、

Off状態がOnではなくSkipped_1に繋がっている
Skipped状態がOnではなくSkipped_1に繋がっている
On状態の自己ループだった辺がOnではなくSkipped_1に繋がっている

というものです。Skipped_1状態は無効発動機会を表しているので、Skipped状態と同様、On状態への遷移確率を1に設定します。

さて、 $m=1$ の場合と同様にして、この場合はSkipped_1の状態の訪問数を数えることで、特技の発動数を数えることができます。その他の状態は、無効発動機会か、特技が発動していない状態を表しています。

同様にして、 $m=3$ の場合は以下のようなマルコフ過程を考えます。

f:id:snowharinchan:20190503135546p:plain:h375

この状態遷移図は $m=2$ の時とほぼ同様なのですが、

Skipped_1がOn状態ではなくSkipped_2状態に繋がっていて、Skipped_2がOn状態に繋がっている

という違いがあります。また、この場合もSkipped_1状態という一つの状態への訪問数を数えることで、楽曲中の特技発動の回数を数えることができます。

これらの状態遷移図の確率遷移行列は以下のように表すことができます。

$\displaystyle A_{m=2} = \begin{bmatrix} 1-p & p & 0 & 0 \\\ 0 & 0 & 1 & 0 \\\ 1-p & p(1-q) & 0 & pq \\\ 0 & 1 & 0 & 0 \end{bmatrix}$

$\displaystyle A_{m=3} = \begin{bmatrix} 1-p & p & 0 & 0 & 0 \\\ 0 & 0 & 1 & 0 & 0 \\\ 0 & 0 & 0 & 1 & 0 \\\ 1-p & p(1-q) & 0 & 0 & pq \\\ 0 & 1 & 0 & 0 & 0 \end{bmatrix}$

これらの行列の固有値1に対する左固有ベクトルを求める *6 *7 と、

$\displaystyle x_{{\rm eig},m=2} = \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & \frac 1 {pq} & 1 \end{bmatrix}$

$\displaystyle x_{{\rm eig},m=3} = \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & \frac 1 {pq} & \frac 1 {pq} & 1 \end{bmatrix}$

となります。比較のため $m=1$ の場合を再掲すると、

$\displaystyle x_{{\rm eig},m=1} = \begin{bmatrix} \frac{1-p}{p^2 q} & \frac 1 {pq} & 1 \end{bmatrix}$

となります。

さて、 $m=1,2,3$ の全ての場合において、特技の発動回数を数えるにあたって訪問数を数えるべき状態の成分は $1/(pq)$ になっています。 $m=1,2,3$ のそれぞれの場合において $1/(pq)$ が $1,2,3$ 回現れているので、次のような 一般式 を考えることができます。

2019.5.12追記： 上記のことは、数学的帰納法によって簡単に示すことができます。

$\displaystyle \begin{aligned} \hat N_{\rm trigs} &\approx N_{\rm maxtrigs} \pi_s x_1^\top \\\ &= N_{\rm maxtrigs} \frac {\frac 1 {pq}} {\left(\frac{1-p}{p^2 q} + m \frac 1 {pq} + 1\right)} \\\ &= N_{\rm maxtrigs} \frac {p} {\left(1-p + m p + p^2 q\right)} \\\ &= \frac {p} {1 + (m-1) p + p^2 q}N_{\rm maxtrigs} \\\ \end{aligned}$

以前との違いは、分母の $\frac 1 {pq}$ に $m$ が掛かって $m \frac 1 {pq}$ になったことです。これをカバー率の期待値に代入すると、

$\displaystyle \begin{aligned} \hat r_{\rm cover} &:= {\mathbb E} \left[ r_{\rm cover} \right] \\\ &\approx \frac{\hat N_{\rm trigs} T_{\rm eff}}{T_{\rm music}}\\\ &= \frac {p} {1 + (m-1) p + p^2 q} \frac {N_{\rm maxtrigs} T_{\rm eff}} {T_{\rm music}} \end{aligned}$

が得られます。

$m=1$ を代入すると、分母の第二項が0になり、以前求めた結果と一致することがわかります。

最終的な理論式

以上の結果を整理します。

$T_{\rm music}$ ：曲の長さ。
$N_{\rm notes}$ ：曲のノーツ数。
$T_{\rm eff}$ ：特技の発動持続時間。
$N_{\rm trig}$ ：特技の発動に必要なノーツ数。
$T_{\rm trig} := T_{\rm music}N_{\rm trig}/N_{\rm notes}$ ：特技の発動に必要な時間。
$N_{\rm maxtrigs} := \left\lfloor N_{\rm notes}/N_{\rm trig}\right\rfloor$ ：曲中の発動機会の総数。
$r_{\rm eff,trig} := T_{\rm eff} / T_{\rm trig}$ .
$r_{\rm eff,trig}$ の小数部分を $q:=r_{\rm eff,trig} - \left\lfloor r_{\rm eff,trig} \right\rfloor$ とおく。
$r_{\rm eff,trig}$ の整数部分を $m:=\left\lfloor r_{\rm eff,trig} \right\rfloor$ とおく。

すると、全楽曲にわたって特技が発動している時間の割合、すなわちカバー率 $r_{\rm cover}$ の期待値 $\hat r_{\rm cover}$ の近似値は、以下の式で与えられます。

$\displaystyle \hat r_{\rm cover} \approx \frac {p} {1 + (m-1) p + p^2 q} \frac {N_{\rm maxtrigs} T_{\rm eff}} {T_{\rm music}}$

2019.5.5 追記： 「この記事の内容」にて触れた通り、この理論式を少しだけ改良したバージョンを作りました。改良された理論式では、予測誤差が4.0ポイントから1.6ポイント程度に改善されています。

snowharinchan.hatenablog.com

理論の検証

f:id:snowharinchan:20190503135530p:plain

f:id:snowharinchan:20190503135527p:plain

さて、具体的な式を導出したからには、実際のデータに当てはまるかを確認する必要があります。

上記のグラフは冒頭で掲載したものと同様のものです。ここでは、真値のデータとして、おくーやさんのブログ記事中の6枚の画像に掲載されている値をOCRによって文字起こししたものを用いています。 *8

１枚目の図は、理論値そのものを真値と比較した様子、２枚目の図は、理論値に補正係数を掛けたものを真値と比較した様子になっています。

まず１枚目の図において、真値と予測値の間のRMSEを計算すると、0.0401程度になります。これは冒頭でも紹介した通り、理論式の平均の予測誤差がおよそ4ポイント程になっている事を表しています。

まず、これを「整数部分が2より大きい場合」の節の冒頭にて掲載したグラフと比較すると、 $m=2$ と $m=3$ の場合において見事に予測ができていることがわかります。これによって、先ほどの理論は $m=2,3$ の場合をうまく説明していると言えます。

また、１枚目の図を見ると、予測値が真値に比べて大きめの値を予測している傾向があることがわかります。そこで、予測値と真値の間の切片0とした時の回帰直線を最小二乗法によって求めると、 $0.9592954$ という値の係数が得られます。この係数を予測値に掛けたものを新たな予測値として、真値と比較したものが2枚目の図になります。ここで、通常このようなパラメータ決定を行う際は、まずデータを学習データとテストデータに分割し、学習データでパラメータを決定した上で、テストデータにおいて予測モデルを評価しなければならないのですが、簡単のため全てのデータを用いています。その事を承知の上でRMSEを計算すると（これは、訓練誤差を測っていることに相当します）、0.0210という値を得ます。

このような補正係数によって予測値が改善されるということから、現実と理論において何かしらのギャップが存在することがわかります。本稿の最後では、より精密な理論式を導出するための方針について述べます。一方で、そのような補正係数を入れない場合においても予測誤差がある程度小さいことから、このギャップは十分小さく、この理論式はパターンCという現象を概ね捉えていると言えます。

外れ値

ここで、これらの図において、 $m=2$ の場合において1つだけ外れ値があることがわかります。このデータを調べると、ノーツ数600、秒数120、20ノーツおきに64%の確率で8秒発動する特技であることがわかります。この特技の $m$ と $q$ を計算すると、 $m=2,\ q=0$ となり、小数部分がちょうど $0$ であることがわかります。補正係数を用いない場合の予測誤差を見ると-0.088になっており、絶対値で見るとデータの中では最も予測誤差が大きくなっています。

そのため、 $q=0$ が原因となっているかについて考えます。

まずこのような状況の場合、特技終了の瞬間が必ず発動機会と重なるため、その発動機会を無効発動機会とみなすか、有効な発動機会とみなすかが実装に依存する ためであると考えられます。この事を考慮して、この場合のみ特別に $m=1,\ q=1$ として計算すると、予測誤差が0.0396となり、他のデータと同程度に説明できるようになりました。 この変更は、 $q=1$ としていることから、必ず無効発動機会を表す状態に遷移する ような実装を採用していることに対応しますが、この変更によってデータが説明可能になったということは、元のモンテカルロ法による実験もそのような実装が用いられていたと推測することができます。

また、全データを見ると、 $q=0$ となっているデータがもう一点だけありました。こちらはノーツ数500、秒数120、32ノーツおきに81%の確率で8秒間発動する特技で、 $m=1,\ q=0$ となっています。しかし、こちらのデータの予測誤差は0.0299となっていて、試しに $m=0,\ q=1$ と修正すると、逆に予測誤差が0.177となり、予測がずれてしまいました。これは次のように解釈できます。 $m=1$ の場合、特技発動の開始地点と終了地点は必ず（時間を $T_{\rm trig}$ 単位で測った時に）整数点上になり、予約で特技が発動しても同様です。したがって、特技が発動していたとしても次の発動機会は必ず特技終了の瞬間になるので、これが無効発動機会とみなされない実装がモンテカルロ法において用いられていたことが考えられます。

議論

以下では、

近似の解釈
分散の評価
どのような部員がパターンCの影響を受けるか
精密な理論式の構築の方針

について考えます。

近似の解釈

まず、定常状態を仮定したことで、最終的に特技発動回数の分布を生起確率 $p/({1 + (m-1)p + p^2q})$ 、回数 $N_{\rm maxtrigs}$ 回の 二項分布で近似 したことになります。また、

$\displaystyle \hat N_{\rm trigs} \approx \frac {p} {1 + (m-1)p + p^2q} N_{\rm maxtrigs}\\$

と表されていますが、これは発動機会の総数に発動確率 $p$ を掛ける、というナイーブな予測

$\displaystyle \hat N_{\rm trigs,naive} := pN_{\rm maxtrigs}\\$

において、発動確率 $p$ が補正係数 $1/(1+(m-1)p+p^2q)$ 倍に縮められているものであるとみなしている、と考えることができます。

つまり パターンCが発生しない場合と比べた時のパターンCの効果は、平均的な発動確率が $1/(1+(m-1)p+p^2q)$ 倍に縮められることに相当する と考えることができます。したがって、この補正係数の大きさを分析することで、部員がパターンCから受ける影響を定量的に解析することができます。

2019.5.12追記： $N_{\rm trigs,naive}$ は曲全体の演奏時間よりも長くなり得るため、この方法ではカバー率が1以上に計算されることになります。この補正係数はそのような計算方法に対するパターンCのカバー率の補正率を表しているため、「カバー率が1以上にならない場合に対してパターンCがどれだけカバー率が小さくなるか」という、本来知りたい値よりも強い補正が計算結果として現れてしまいます。

分散の評価

この形式で表すことの利点の一つは、分散の評価ができることです。定常状態を仮定したことで、に特技発動回数の分布を二項分布で近似したことになるため、この二項分布の分散を求めることで、

$\displaystyle \begin{aligned} {\mathbb V} \left[ r_{\rm cover} \right] &\approx {\mathbb V} \left[ \frac{N_{\rm trigs} T_{\rm eff}}{T_{\rm music}} \right] \\\ &\approx {\mathbb V} \left[ N_{\rm trigs} \right] \frac {T_{\rm eff}^2} {T_{\rm music}^2} \\\ &= N_{\rm maxtrigs} \frac {p} {1 + (m-1)p + p^2q} \\\ &\ \ \ \ \ \times \left(1 - \frac {p} {1 + (m-1)p + p^2q}\right) \frac {T_{\rm eff}^2} {T_{\rm music}^2} \\\ &= N_{\rm maxtrigs} \frac {p \left(1+(m-1)p + p^2q - p\right)} {(1 + (m-1)p +p^2q)^2} \frac {T_{\rm eff}^2} {T_{\rm music}^2} \\\ &= N_{\rm maxtrigs} \frac {p \left(1+(m-2)p + p^2q\right)} {(1 + (m-1)p +p^2q)^2} \frac {T_{\rm eff}^2} {T_{\rm music}^2} \\\ \end{aligned}$

となります。

500ノーツ120秒の楽曲に対して先ほどの $N_{\rm trig}=25,\ p=0.81,\ T_{\rm trig}=8.5,\ m=1,\ q=0.42$ の特技を当てはめると、分散は $0.0232$ 程度になります。標準偏差に直すと $0.152$ になるため、カバー率は15%程度変動すると見込むことができます。

どのような部員がどの程度パターンCの影響を受けるか

ここでは、補正係数 $1/(1+(m-1)p+p^2q)$ に関する議論を行います。補正係数を

$\alpha := 1/(1+(m-1)p+p^2q)$ ：補正係数

と定義しておきます。

まずわかるのは、補正係数の分母が大きいほど大きな補正を受けるということです。したがって、分母が大きくなる要因を考えていきます。

この理論式において最も特徴的なのは、 $q$ が大きい部員ほど強い補正を受ける ということです。これは、以前考えた数直線上の線分の長さを考えると、 $q$ が大きいほどスキップされる発動機会が多いということなので、この点を考えると自然なことですが、これが理論式にはっきりと現れていること、そして その $q$ の値をどれほど心配しなければならないのかという定量的な分析 ができることが重要です。

また、 $m$ が大きい部員ほど強い補正を受けることもわかります。 $m$ は１回の特技発動時間に含まれる発動機会の数を表しているので、 $m-1$ は 必ずスキップされる状態の数 を表しています。これが大きいほど補正が大きくなるのは自然です。実際理論式と真値の比較のグラフを見ると、全体として $m=2,3$ の場合の方が $m=1$ よりもカバー率が小さいことが見て取れます。

また、そもそも発動確率が小さい部員は受ける影響が少ないことがわかります。

さて、 $q$ の大小によって受ける補正の大きさが違うことはわかったのですが、難しいのは $q$ は部員のパラメータだけでなく、曲のノーツ数と時間に依存して変化する ことです。また、この解析ではノーツが等間隔に分布していると仮定しましたが、ノーツ密度が高い地帯では $q$ が上昇します。ここで、 $q$ は何かに単純に比例するわけではなく、発動持続時間と発動機会の時間間隔の比の 小数部分で定まる ことから、直感的な予測が難しいのですが、ノーツ密度が少しだけ上昇した時、 $r_{\rm eff,trig}$ の整数部分は変化せず、小数部分のみが上昇することから、微小なノーツ密度の変化の範囲では、概ね $q$ はノーツ密度に比例すると考えられます。

f:id:snowharinchan:20190503135536p:plain

上図は、データ中の全部員について補正係数 $\alpha$ を計算した時の分布です。この図より、パターンCによって概ね0.8倍程度にカバー率が制限されていることがわかります。また、 $m=2,3$ の部員は、 $m=1$ の部員に比べて補正が強い傾向にあることがわかります。

f:id:snowharinchan:20190503135550p:plain

上図は、 $q$ の値と補正係数の関係のプロットです。上図より、ほぼ全ての部員について、 $q$ （と $m$ ）の値を見れば概ね補正係数が特定できることがわかります。具体的には、

$m=1$ の場合は、補正係数が $1$ から $0.7$ の間で分布しており、 $q$ に概ね比例してこの間を推移する。すなわち、 $\alpha \approx 1-0.3q$ と表せる。
$m=2$ の場合は、補正係数はほぼ $0.6$ で一定である。
$m=3$ の場合は、データが少ないので定かではないが、唯一のデータについては補正係数は $0.5$ となっている。

精密な理論式の構築の方針

本稿では、マルコフ過程によるモデルを構築した際に様々な近似を行いましたが、ここでは、モンテカルロ法を用いて計算された真値と同様に精密な期待値の理論値が求められると期待できる方針について述べます。

本稿においてマルコフ過程のモデルを構築した際、予約によって発生する無効発動機会はランダムに発生するという近似を行いました。しかし、この現象は以下のようにして取り払うことで、現象をより精密に捉えることができます。

全ての発動機会の数 $N_{\rm maxtrigs}$ 分だけ状態を用意する。これに特技が発動していない Off 状態を加えた、合計 $N_{\rm maxtrigs} + 1$ 状態を用意する。
$N_{\rm maxtrigs}$ 回連続で特技が発生した時の、無効発動機会が発生するタイミングを数える。これらの無効発動機会に対応した状態を作る。

これらの操作によって、具体的には本稿にて $m=3$ の時考えた状態遷移図におけるSkippedからSkipped_1へのループを取り払い、後退方向の辺をOff状態への辺のみに制限しながら、状態を直線状に伸ばして繋いでいった状態遷移図が出来上がります。この時の確率遷移行列のサイズは $(N_{\rm maxtrigs} + 1) \times (N_{\rm maxtrigs} + 1)$ となります。

また、定常状態を仮定するという近似を取り払うには、

定常状態を仮定するのではなく、 $x_0 A^k x_{\rm eff}$ の総和によって計算する

ということが考えられます。

以上の計算は、主に $(N_{\rm maxtrigs} + 1) \times (N_{\rm maxtrigs} + 1)$ 行列の最高で $(N_{\rm maxtrigs} + 1)$ 乗の計算を $N_{\rm maxtrigs}$ 回行うことでできますが、 $N_{\rm maxtrigs}$ は多くても高々40程度なので、比較的少ない計算量で期待値を計算することができるはずです。

結論

マルコフ過程を用いて「スキップされる発動機会」と「有効な発動機会」などを表現し、系が定常状態にあるという近似を用いて、部員ステータスから簡単にカバー率を求めることのできる近似的な理論式を導いた。
パターンCの効果は、平均的な発動確率が部員のパラメータに依存した係数倍（ $1/(1+(m-1)p+p^2q)$ ）に縮められることに相当する ことを導いた。
このことから、特技の持続時間を特技発動機会の時間間隔で割った値の小数部分が大きい部員ほど、パターンCによって強い補正を受ける ことを見た。
上記は近似的な理論であったが、モンテカルロ法を用いて計算された真値と同様に精密な期待値の理論値を求める方針 について述べた。

snowharinchan.hatenablog.com