2019-05-12

もっと簡単に例の近似式を求める方法

スクフェス

この記事の近似式をもっと簡単に出す方法です。

snowharinchan.hatenablog.com

導出

特技の持続時間を $T_{\rm eff}$ , 発動に要する時間を $T_{\rm trig}$ とし、 $T_{\rm eff}/T_{\rm trig}$ の整数部分を $m$ , 小数部分を $q$ とします。

まずは簡単のために $m=1$ の場合を考えます。後に $m \geq 2$ の場合を考えます。

今、曲が終わって、全部で $N_{\rm maxtrigs}$ 回の発動機会があったとします（特技が発動しなかった機会や、無効発動機会も含みます）。ここで、

有効発動機会において特技が発動しなかった回数を $N_{\rm Off}$
有効発動機会において特技が発動した回数を $N_{\rm On}$
無効発動機会の回数を $N_{\rm Skipped}$

とします。このように設定することで、

$\displaystyle N_{\rm Off} + N_{\rm On} + N_{\rm Skipped} = N_{\rm maxtrigs}$

となります。 $N_{\rm On}$ が何回であったのかを調べることで、 $T_{\rm eff} N_{\rm On}$ で合計の発動時間がわかり、 $T_{\rm eff} N_{\rm On} / T_{\rm music}$ （ただし $T_{\rm music}$ は曲全体の時間）でカバー率がわかるので、以下では $N_{\rm On}$ を求めます。

さて、有効発動機会だけに限れば、有効発動機会における特技の発動率は $p$ なので、平均的には

$\displaystyle N_{\rm Off} : N_{\rm On} = 1-p : p$

となります。

さて、この $N_{\rm On}$ の発動機会たちなのですが、この中には

自分をカバーする特技は無効発動機会を含んでいなかった
自分をカバーする特技は無効発動機会を含んでいた

という2つに分かれます。

そこで、これら全ての $N_{\rm On}$ の中から、後者の無効発動機会を含んでいたような発動機会たちを特定する必要があります。

そのために、ここで $N_{\rm On}$ の発動機会たち一人一人に次のような質問をしていきます。

「あなたの発動機会をカバーしている特技があると思いますが、その特技の更に１個後に連続で特技は発生しましたか？」

この質問に対して、平均的には 確率 $p$ で「はい」という答えが得られます。

そして、「はい」という答えが得られた人には更に質問を続けます。

「なるほど、それでは、その連続で発生した特技は、無効発動機会を一つ含んでいましたか？」

この質問に対して、平均的には 確率 $q$ で「はい」という答えが得られます。（理由は前回の記事をご参照ください。）

以上によって、 $N_{\rm On}$ の発動機会たちの中で、両方の質問に対して「はい」と答える発動機会の割合は $pq$ です。したがって、実は $N_{\rm Skipped} = pqN_{\rm On}$ である ことがわかります。

よって、

$\displaystyle N_{\rm On} : N_{\rm Skipped} = 1 : pq$

です。

これを先ほどの比と合わせると、

$\displaystyle N_{\rm Off} : N_{\rm On} : N_{\rm Skipped} = 1-p : p : p^2q$

となります。左辺は合計して $N_{\rm maxtrigs}$ になるので、

$\displaystyle \begin{aligned} N_{\rm On} &= \frac{p}{1-p + p + p^2q} N_{\rm maxtrigs} \\ &= \frac{p}{1 + p^2q} N_{\rm maxtrigs} \end{aligned}$

がわかりました。

整数部分が2以上の場合

$m \geq 2$ 以上の場合、特技が発動すると必ず $m-1$ 回余分に無効発動機会が発生します。しかし、特技の発生区間中に含まれ得る発動機会は $m$ 回か $m-1$ 回のいずれかであることは変わりありません。

したがって、 $(m-1)N_{\rm On}$ 回分余分に発生する無効発動機会を加えると、

$\displaystyle N_{\rm Skipped} = pqN_{\rm On} + (m-1)N_{\rm On}$

となるので、

$\displaystyle N_{\rm On} : N\_{\rm Skipped} = 1 : pq + (m-1)$

です。したがって、 $N_{\rm Off} : N_{\rm On} = 1-p:p$ であることは変わらないので、

$\displaystyle N_{\rm Off} : N_{\rm On} : N_{\rm Skipped} = 1-p : p : p^2q + (m-1)p$

となるので、

$\displaystyle \begin{aligned} N_{\rm On} &= \frac{p}{1-p + p + p^2q + (m-1)p} N_{\rm maxtrigs} \\ &= \frac{p}{1 + (m-1)p + p^2q} N_{\rm maxtrigs} \end{aligned}$

という式が得られました。

マルコフ連鎖の話との繋がり

$m=1$ の場合を考えます。上記の話で出てきた2つの式は、

$\displaystyle \begin{bmatrix} 1-p & 1-p \\ 0 & pq \end{bmatrix} \begin{bmatrix} N_{\rm Off} \\ N_{\rm On} \end{bmatrix} = \begin{bmatrix} N_{\rm Off} \\ N_{\rm Skipped} \end{bmatrix}$

のように表すことができます。

左辺に $N_{\rm Skipped}$ を、右辺に $N_{\rm On}$ を出すために変形すると、

$\displaystyle \begin{bmatrix} 1-p & 1-p & 0 \\ p & p & 0 \\ 0 & pq & 0 \end{bmatrix} \begin{bmatrix} N_{\rm Off} \\ N_{\rm On} \\ N_{\rm Skipped} \end{bmatrix} = \begin{bmatrix} N_{\rm Off} \\ N_{\rm On} \\ N_{\rm Skipped} \end{bmatrix}$

と書けます。

ここで $N_{\rm Skipped} = pqN_{\rm On}$ なので、同じものを引いて足すことで

$\displaystyle \begin{bmatrix} 1-p & 1-p & 0 \\ p & p(1-q) & 1 \\ 0 & pq & 0 \end{bmatrix} \begin{bmatrix} N_{\rm Off} \\ N_{\rm On} \\ N_{\rm Skipped} \end{bmatrix} = \begin{bmatrix} N_{\rm Off} \\ N_{\rm On} \\ N_{\rm Skipped} \end{bmatrix}$

です。

ここで、 $N_{\rm Off}/N_{\rm maxtrigs}$ は「ある発動機会が $\rm Off$ である確率」と読めるので、 $N_{\rm Off}/N_{\rm maxtrigs} \approx P_k(s={\rm Off})$ が言えます。したがって、両辺を $N_{\rm maxtrigs}$ で割ると

$\displaystyle \begin{bmatrix} 1-p & 1-p & 0 \\ p & p(1-q) & 1 \\ 0 & pq & 0 \end{bmatrix} \begin{bmatrix} P_k(s={\rm Off}) \\ P_k(s={\rm On}) \\ P_k(s={\rm Skipped}) \end{bmatrix} = \begin{bmatrix} P_k(s={\rm Off}) \\ P_k(s={\rm On}) \\ P_k(s={\rm Skipped}) \end{bmatrix}$

となります。 $\pi_s := \left[ P_k(s={\rm Off}) \ P_k(s={\rm On}) \ P_k(s={\rm Skipped}) \right]$ と書いて、両辺を転置すると

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

となります。

このようにして、今回の話は前回の記事のマルコフ連鎖の定常状態での話と繋がります。

今回の話は、この流れを逆に辿ることで導出しました。 $N_{\rm Skipped}$ の導出の話はややこしいのですが、マルコフ連鎖の話はその検算になっていると言えます。また、初期時刻付近での特技の発生確率についてはマルコフ連鎖のモデルの方が詳しく記述できているので、その影響の大きさを調べるためにはそちらのモデルの方が有効です。

2019-05-05

スクフェスの「パターンC」の理論式の改良

スクフェス

この記事では、以下の記事の理論式を改良します。

snowharinchan.hatenablog.com

パターンCにおけるカバー率の近似式 - 改良版

以下の記号を用います。

$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 rm,1} - T_{\rm rm,2}} {T_{\rm music}}$

ここで、 $T_{\rm rm,1},\ T_{\rm rm,2}$ はそれぞれ特技の余剰発動時間で、

$\displaystyle \begin{aligned} T_{\rm rm,1} &= T_{eff} - \left( \frac {N_{\rm notes}} {N_{\rm trig}} - \left\lfloor \frac {N_{\rm notes}} {N_{\rm trig}} \right\rfloor \right) T_{trig} \\\ T_{\rm rm,2} &= \begin{cases} T_{\rm rm,1} - T_{trig} & T_{\rm rm,1} \geq T_{trig} \ \text{and} \ m=1 \\\ 0 & \text{otherwise} \end{cases} \end{aligned}$

です。

この理論式の大部分は、冒頭にて紹介した記事の導出に基づいています。そのため、詳しい導出についてはそちらをご覧ください。この記事では、前回の記事との差分の部分について紹介します。

導出

これは、「曲の終了間際に発動し、曲の終了後も続いた特技の発動時間は、カバー率に影響しない」という効果*1を考慮したものです。 $m \geq 1$ となる曲は、少なくとも一番最後の発動機会に発動した特技は、曲の終了後も効果が続いてしまいます。この時、曲の終了後も続いた時間を余剰時間と呼ぶことにします。前回の記事の理論値は、特技の発動回数×持続時間で全体の発動時間を計算していたのですが、この方法では曲の終了後も続いている時間も数えてしまっていました。この余剰時間は必ず発生するので、この時間を取り除くことでより良い理論値が得られることが期待されます。

一番最後の発動機会に特技が発動したときの余剰時間が $T_{\rm rm,1}$ 、一番最後から2番目の発動機会に特技が発動したときの余剰時間が $T_{\rm rm,2}$ です。曲のノーツ数を発動に必要なノーツ数で割った時の値が $34.4$ の場合、一番最後の発動機会では、発動機会 $0.4$ 回分の時間、つまり $0.4 T_{\rm trig}$ しか特技が持続しません。そのため、余剰の時間は $T_{\rm eff} - 0.4 T_{\rm trig}$ となります。

冒頭にて紹介した前回の記事の導出によると、任意の時間で特技が発動する確率は $\frac {p} {1 + (m-1) p + p^2 q}$ で与えられると近似できるので、余剰時間の期待値は $\frac {T_{\rm rm,1}} {1 + (m-1) p + p^2 q}$ となります。この値を前回導出した判定の持続時間から引き去ることで、余剰時間を取り除きます。

$T_{\rm rm,2}$ については、時刻が最後の発動機会に比べて発動機会1回分だけ後ろにずれているので、余剰時間が発生する場合は $T_{trig}$ 分だけ長くなります。ただし、 $T_{\rm rm,1} \lt T_{\rm trig}$ となり、特技が曲の時間中に収まる場合もあるため、この場合は $T_{\rm rm,2} = 0$ とします。

また、 $m=2$ 以上の場合2つ目の余剰時間は考慮していないのですが、これは結果的にその方が予測誤差がが下がったため、特に結果的に $m=2,3$ の場合が上手く説明できているためという理由です。 $m=2,3$ の場合の方が余剰時間は多く発生するので、一見これらの場合こそ $T_{\rm rm,2} = 0$ の効果を考えた方が良いように思えるのですが、なぜかこの項を取り除いた方が数値上は若干良い結果が得られているので、もしかしたら尤もらしい理由があると考えて、ここではこれらの項を取り除いています。

理論の検証

図2: T_{rm,2}をm=1の場合も入れた場合。 — 図2: $T_{rm,2}$ を $m=2,3$ の場合も入れた場合。

図1は理論式とモンテカルロ法によるシミュレーション値の比較です。予測誤差（RMSEで評価）は0.0158となっています。また、ここでは前回の記事にて紹介した外れ値について、 $m=1,q=1$ となるように計算しています。この場合を元の通りに計算すると予測誤差は 0.0198 と少しだけ増えますが、以前の理論値よりも減少しています。また、前回の記事の通りに計算した補正係数は0.996となり、全体として大きめの値を予測してしまうという前回の理論値における傾向が取り除かれていることがわかります。

図2は $T_{rm,2}$ を $m=2,3$ の場合も入れた場合の結果です。こちらは理論としてはより整合性があるのですが、予測誤差が0.0167と少し大きくなります。前回の記事の通りに計算した補正係数は0.998となり、同じく大きめに予測する傾向は取り除かれています。

結論

曲の終了付近に発動した特技は、曲の終了後も持続してしまうが、この時の余剰時間はカバー率に影響しない。このような余剰時間は必ず発生するが、前回の理論値は、この余剰時間を含めてしまっていた。
この余剰時間を取り除いた新たな理論値を導出したところ、前回の理論値に比べて、より小さな予測誤差で予測が可能になった。

免責

本稿における内容を利用した結果について、著者はいかなる責任も負いません。

*1:着想過程については https://twitter.com/siratama_z/status/1124553497502863361 の会話を参照ください。

2019-05-03

スクフェスの「パターンC」に関するマルコフ過程を用いた理論解析と特技性能の近似式

スクフェス

スクフェスの特技において「パターン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$ とおく。

$\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

訂正履歴

2019.5.4: 文中の ${\rm eff}$ と ${\rm trig}$ が入れ替わっているという指摘 *9 を受け、他の該当箇所も含め訂正しました。
2019.5.5: 改良版の理論式に関する情報を掲載。
2019.5.12: いくつかの誤植など*10を修正。また、より簡単な導出を行なった記事に関する情報を掲載。

免責

本稿における内容を利用した結果について、著者はいかなる責任も負いません。

*1:https://twitter.com/natsuowa_6000/status/1123273417870811136

*2:https://twitter.com/Octilleryo/status/1123513256545382400

*3:https://twitter.com/devsifmatch/status/1123843581507198977

*4:https://twitter.com/Octilleryo/status/1123513256545382400

*5:https://www.wolframalpha.com/input/?i=Eivenvector%5BTranspose%5B%7B%7B1-p,p,0%7D,%7B1-p,p*(1-q),p*q%7D,%7B0,1,0%7D%7D%5D%5D

*6:https://www.wolframalpha.com/input/?i=Eivenvector%5BTranspose%5B%7B%7B1-p,p,0,0%7D,%7B0,0,1,0%7D,%7B1-p,p*(1-q),0,p*q%7D,%7B0,1,0,0%7D%7D%5D%5D

*7:https://www.wolframalpha.com/input/?i=Eivenvector%5BTranspose%5B%7B%7B1-p,p,0,0,0%7D,%7B0,0,1,0,0%7D,%7B0,0,0,1,0%7D,%7B1-p,p*(1-q),0,0,p*q%7D,%7B0,1,0,0,0%7D%7D%5D%5D

*8:再びになりますが、このような解析を行なって頂き、ありがとうございます。

*9:https://twitter.com/siratama_z/status/1124217396938461184

*10:https://twitter.com/Darsein_LL/status/1125807874909835265

2019-03-04

スクフェスのスコア計算方法の予想をしていたら、浮動小数の深みを見た話 - 「float32では合うのに、float64では合わない計算」

概要

スクフェスのスコア計算方法の予想を多くの方としていた所*1、float32では合うのに、float64では合わない計算に直面しました。

一言で言えば浮動小数誤差なのですが、

ceil(x * 1.1) という計算において、float64において、特定の10の倍数のみ計算が合わない
float64では計算が合わない場合があるのに、float32では常に計算が合う
上記の現象が、スクフェスのスコア計算シミュレーションという実用的な場面で発生している

という面白い側面があります。

この現象の技術的な詳細は、別のブログ記事にて非常に細かく紹介しています。この現象は、計算結果の末尾たった1ビットがずれることで発生していることがわかりました。この1ビットがずれる場合とずれない場合が存在するのですが、その詳細について上記の記事で詳しく説明しています。

この記事では、この現象に辿り着くまでの議論の経緯を、スクフェス的な側面にスポットライトを当てて説明します。

背景

2019年2月20日のスクフェスのアップデートにおいて様々な変更が行われました。その大きなものの一つに、回復特技と判定特技によってタップ時のスコアが増加するようになった、というものがあります。

@magu_maki さんが、下記のtwitterのスレッドにおいて非常に詳細なデータと共に、回復特技のタップスコア計算方法の詳細なデータと共に、スコア計算方法の予想を載せています。

回復の新機能について検証しました。

体力が41,43,44の3つのユニット(属性一致、グループ不一致、回復特技のみ)で試したところ、それぞれ2.04%,2.18%,2.25%×ゲージがMAXになった回数分ユニットの属性値が増えるとして計算すると、実測値と全て一致しました。

⇒上昇率は0.5+(体力-19)×0.07と予想！ pic.twitter.com/j7xTUzBOdb
— まぐお＠ﾌﾗｲﾝｸﾞﾊﾞｰﾄﾞ=ｴﾘｺﾇｽⅡ世 (@magu_maki) 2019年2月20日

このデータと予想は非常に精密で、

データの数が非常に多い。
それぞれ体力ゲージのMAX回数やユニット属性値、Perfect/Great判定など様々な状況での観測である。
それにも関わらず、それらの数多くの実測値を全て言い当て続けている。

という特徴を持っていました。

しかし、この驚くべき精度で次々と実測値が的中される中、たった一つだけ実測値と予想値が合わない観測が発生していました。

応援タップSCOREアップの計算方法はまだ解が見つかっていないけど、とりあえず回復のパラメータアップの方の検証を進め、体力41～54まで判明。なお、応援タップSCOREアップは、SCOREの繰り下げ前の小数にアップ率を掛けて繰り上げると良さそうだけど、1つだけ合わない実測値ががががが… pic.twitter.com/9bsTnL18Nm
— まぐお＠ﾌﾗｲﾝｸﾞﾊﾞｰﾄﾞ=ｴﾘｺﾇｽⅡ世 (@magu_maki) 2019年2月22日

上記のツイートの3枚目の画像の下側にある、赤枠で囲まれた「-1」のデータがその観測値です。

このデータでは、以前のデータとは違う点が一つありました。それは、タップスコア計算において、「タップスコアアップ」というメドレーフェスティバルのイベント中に発生する「応援」の効果を受けて、通常よりタップスコアが上昇するというものでした。

タップスコアアップは、タップスコアが1.1倍になるという効果を持っているということが広く知られています。@magu_maki さんの計算によると、

通常のタップスコア計算後、つまり切り捨て処理などを行った後に、得られたスコアを1.1倍して、その後切り上げをする

という計算を行うと、実測値とほぼ一致することがわかりました。しかし、1つのデータのみこの計算では説明ができていませんでした。

唯一説明できなかったデータにおける実際のスコア画面、ユニット編成、タップスコアアップの様子。ユニット属性値から予想されるスコアは902だったが、スコア実測値は903だった。画像は@magu_makiさんよりご提供頂きました。 — 唯一説明できなかったデータにおける実際のスコア画面、ユニット編成、タップスコアアップの様子。ユニット属性値から予想されるスコアは902だったが、スコア実測値は903だった。画像は@magu_makiさんよりご提供*2頂きました。

f:id:snowharinchan:20190304204009j:plain — 唯一説明できなかったデータにおける実際のスコア画面、ユニット編成、タップスコアアップの様子。ユニット属性値から予想されるスコアは902だったが、スコア実測値は903だった。画像は@magu_makiさんよりご提供*2頂きました。

試みられた計算: 10の倍数が関係している？

この現象を説明するために、様々な計算が試みられました。

下記のツイートでは、@magu_makiさんによって、切り上げや切り捨てをするタイミングや、行う対象の数を変更する事が試みられました。

+1とは？
他に2通りの計算方法を試したのですがどちらも合わずでした。誤差が出ている箇所のGreatの切捨後のSCOREが820(赤の部分)で1桁目が0なので、その0.1倍が82と整数になり繰り上がらないため902になり、実測値の903と合わないので、その辺り何かヒントにならないですかね？ pic.twitter.com/sfMIc2sCiQ
— まぐお＠ﾌﾗｲﾝｸﾞﾊﾞｰﾄﾞ=ｴﾘｺﾇｽⅡ世 (@magu_maki) 2019年2月22日

この時点で、

「誤差が出ている箇所のGreatの切捨後のSCOREが820(赤の部分)で1桁目が0」で、
「その0.1倍が82と整数になり繰り上がらないため902に」なるために「実測値の903と合わない」のではないか？

という事が@magu_makiさんによって予想されます。

深まる謎: 丸め込み誤差のはずなのに、10の倍数でも正しく丸め込まれる場合がある？

また、@siratama_zさんによって、丸め込み誤差を考慮したような下記の方法が試みられました。

すみません、説明が足らなかった上に言ってることが間違っていました。
「タップSCOREアップ×1.1」から「繰上」の計算を、ROUNDUP(x,0)ではなくROUNDDOWN(x+1,0)でやるとうまくいくかなと思ったのですが、ゲージ2の101-200の「97」は実測値とずれていないことを見落としていました。 pic.twitter.com/duvsx9VePe
— siratama (@siratama_z) 2019年2月22日

このような計算を行うことで、10の倍数の時のみ通常の切り上げよりも1だけ大きい出力を得ることができ、820*1.1 → 903という計算が説明できるようになりました。

しかしここで、重要な事が指摘されます。それは、同じ10の倍数であるはずの他のデータでは、元々切り上げが正しく計算できているため、この計算方法ではそちらのデータについては説明ができなくなってしまうことが指摘されます。具体的には、実測値は

820*1.1 は820+82+1 = 903と計算されている
970*1.1 は970+97 = 1067と計算されている

という不可解な事が起きているのです。

この点は、@magu_makiさんによっても問題として明確に挙げられています。

まさにこの970があることが厄介で、820と同じように1の位が0なのにこちらは合ってるんですよね。何か変えるとこっちが違ってきちゃって…うーむ…🤔
— まぐお＠ﾌﾗｲﾝｸﾞﾊﾞｰﾄﾞ=ｴﾘｺﾇｽⅡ世 (@magu_maki) 2019年2月22日

この時点で、

1.1倍という計算は小数を含む計算なので、入力が10の倍数の時、切り上げ前の結果がちょうど整数になる関係で誤差が発生するのではないか
820では切り上げが正しく行われないが、970では正しく行われるようなメカニズムは何か？

という事が明確に意識されていました。つまりこの問題の最も難しい点は、10の倍数の中でも、正しく切り上げが計算される場合がある事を説明するという点でした。

打開策: 「float32では計算が合うのに、float64で計算すると計算が合わない」事がある

ここで@snowharinchan*3によって、下記の事が指摘されます。

わかったかもしれません！実はスクフェスは内部的にfloat64を使っていて、10の倍数であるはずの820を1.1倍した後にもfloat誤差が乗っていて、それを切り上げて誤差が発生するのではないでしょうか？
— レン@凛ちゃん推し (@SnowhaRinchan) 2019年2月23日

このように仮定して、float32の場合とfloat64の場合、そしてまぐおさんの表にある、1.1倍適用前のスコアが820、970の場合で計算すると、なんとfloat64で計算する場合、820だとceilがずれるが970だとずれない、という驚くべき結果が得られましたw pic.twitter.com/Jfkpc4ybrU
— レン@凛ちゃん推し (@SnowhaRinchan) 2019年2月23日

実際、1.1倍適用前のスコアが0から1999までの場合で、float32とfloat64でのceilの結果がずれる場合を全て列挙するとこのようになり、基本的に10の倍数なのですが、940から1280の間に空白があったり、他にも470-640、60-80などにも空白があります。 pic.twitter.com/4QNFvdfTKU
— レン@凛ちゃん推し (@SnowhaRinchan) 2019年2月23日

なので、1.1倍適用前のスコア実測値がここに列挙されているスコアになり、かつ、1.1倍応援が適用された場合、float誤差によってceil関数の計算がずれる、と予測することができます。
— レン@凛ちゃん推し (@SnowhaRinchan) 2019年2月23日

この指摘には二つの大きなポイントがあります。

１つ目は、float64で計算すると、同じ10の倍数でも、820では切り上げの計算が合わないのに、970では合う現象を再現できる事を指摘した点です。つまり、当初予想されていた通りこの現象は浮動小数計算由来の切り上げ誤差だったが、詳しく計算すると、中には計算が合う場合も存在する事を指摘した点です。

しかし、@magu_makiさんのExcelによる検証でも、1.1倍の処理に浮動小数が使われています。にも関わらず、なぜExcelのシートでは現象が再現されなかったのでしょうか？

その答えとなる２つ目のポイントは、float64ではこの現象は発生するのに、float32ではこの現象は発生せず常に切り上げが正しく計算されるという事を指摘した点です。同じ計算をfloat32で行うと同じ不動小数の仲間でも、常に正しく計算が行われました。そして、どうやらExcelではfloat32が使われているらしく（ExcelのVBA環境にて WorkSheetFunction.Ceiling を用いると、float32における計算結果が再現されることなどから推定）、float32を用いているExcelでは現象が再現できなかったという事がわかりました。

調べてみると、VBAではDouble型を用いてfloat64による計算が可能であることがわかります。VBAで作った関数はExcelから呼べるため、これを使ってfloat64で1.1倍した後に切り上げを行う関数を作成した所、820、970のデータを含め、他のデータも説明が可能になりました。*4

発生する理由：末尾1ビットのみがずれる場合がある

この現象が発生する理由を詳しく追っていくと、float64で計算した場合、末尾の1ビットのみ1になり、切り上げ時に誤差が発生する場合があることが発生原因であることがわかります。

末尾のビットがずれるだけなら、小数計算を厳密に扱えない計算機では普通の事の様に思えますが、面白いのは10の倍数の中でも末尾のビットがずれない場合があるということと、float64ではこれが発生するのに、float32では発生しないということです。

この詳しい原理については、冒頭でも紹介した別の記事にて詳しく説明しています。

ゲーム的なコメント

float64によって切り上げの計算が合わない場合は、全て最終的なスコアの値が本来の値よりも1だけ大きく計算されます。スクフェス的には、ボーナスがタップスコア「アップ」とのみ謳っており、「1.1倍される」とは言っていません。この現象が発生している時も発生していない時も、いずれもボーナスが適用される前よりもスコアが上がっているという事実は変わらないことから、「タップスコアアップ」というボーナスの名称には反していない、と個人的に考えています。

結論、感想

スクフェスのスコア計算において、「1.1倍した後に切り上げ」という計算において切り上げ誤差が発生していた。誤差が発生した数は常に10の倍数であったが、10の倍数の中でも、誤差が発生しない数が存在し、謎を呼んだ。
計算をfloat64で行うことで現象が再現できることがわかった。
同時に、float32で行うと再現できないことがわかった。
float64で計算を行うことで、唯一説明できなかったデータが説明可能になった。

技術的な詳細は、冒頭でも紹介した下記の記事にて記していますが、この切り上げの問題は末尾のたった1ビットの誤差によって発生しています。21世紀である2019年の現在、このようなビットレベルの現象が、実際にiOSやAndroidで配布されているスクフェスというゲームのスコア予想という実用的な場面で発生している事にワクワクを覚えました。

snowharinchan.hatenablog.com

*1:スクフェスの2019年2月のアップデートで変更された新たなスコア計算方法の解析を @magu_maki さんが行っていた所、一つだけ説明できないデータ点が浮上しました。様々な計算方法が試みられましたが、これに対して、私が計算がずれる原因は、浮動小数点由来の誤差ではないかと予想しました。その結果、浮動小数の誤差まで考慮すると、実際に観測されたデータが全て説明できるようになりました。

*2:https://twitter.com/magu_maki/status/1099710995600072704

*3:この記事の著者です

*4:https://twitter.com/magu_maki/status/1099301573279109120

2019-02-24

float32では正しく計算されるが、float64では正しく計算されない場合のある式：ceil(1.1*x) とそのラブライブ！スクフェスのスコア計算における影響

概要

色々な計算を行っていた所*1、次のようなfloat32とfloat64における一般的な現象に直面しました。それは、

ceil(1.1 * x)

という計算は、 float32だと計算が合うが、float64だと計算が合わない場合がある という現象です。つまり、 計算精度が高いはずのfloat64の方が計算が合いにくい 、という現象が起きました。

一言で言えば浮動小数の誤差で、1.1倍している関係から、基本的に計算が合わないのは10の倍数です。しかし、 10の倍数でも、中にはfloat64で計算の合う場合もある ことがわかりました。

例えば、x が

230の時、float32では ceil(1.1 * x) が253 となるのに対し、float64では 254 と正しく計算されない
しかし540の時、 float32でも 594 、float64でも 594 と、正しく計算される
この現象は、Python（numpy使用）、C、Javaで確認される

試しに0から1500未満の整数でこの現象が起きる（つまり、float64で計算が合わない）ような数を列挙すると、

50 90 100 110 170 180 190 200 210 220 230 330 340 350 360 370 380 390 400 410 420 430 440 450 460 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390 1400 1410 1420 1430 1440 1450 1460 1470 1480 1490

のように、一見不規則な列になります。

もっとも驚くべきことは、このような現象が 実用上の場面 で起きているということです。具体的には、後に紹介する通り、@magu_maki さんの観測データを用いて、スクフェスの点数計算を手元でシミュレーションしようとした際、float32とfloat64で計算結果が異なり、float64で計算しないと挙動が再現できなかった、という点です。

ここでは、

なぜfloat64で計算が合わない場合があるのか
計算が合わないのは基本的に入力が10の倍数の時であるように見えるが、なぜ10の倍数の中でも、計算が合う時があるのか

について、float型のデータ表現方法や2進数展開などについて詳しく見ながら分析します。

概要
- 目次
計算が合わない場合の列挙
1.1のfloat32とfloat64における表現
1.1を10倍した時
1.1を一般の10の倍数倍した時
四捨五入でも計算が合わない場合
いつ誤差が発生するのか
float32の場合
実応用：ラブライブ！スクフェスにおけるスコア計算
結論

計算が合わない場合の列挙

まずは、冒頭の計算がfloat32とfloat64の場合で異なるような正の整数xを列挙します。試しに、0から1500未満の整数について調べます。

以下のコードを使います。

public class DoubleFloat {
    public static void main(String[] args) {
        float a_32 = 1.1f, b_32;
        double a_64 = 1.1d, b_64;
        for (int i = 0; i < 1500; i++) {
            b_32 = i;
            b_64 = i;
            if (Math.ceil(a_32*b_32) != Math.ceil(a_64*b_64)) {
                System.out.printf("%d %f %f\n", i, Math.ceil(a_32*b_32), Math.ceil(a_64*b_64));
            }        
        }
    }
}

すると、以下の結果を得ます。

50 90 100 110 170 180 190 200 210 220 230 330 340 350 360 370 380 390 400 410 420 430 440 450 460 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390 1400 1410 1420 1430 1440 1450 1460 1470 1480 1490

Math.ceil の結果も出力するように、9行目を

System.out.printf("%d %f %f\n", i, Math.ceil(a_32*b_32), Math.ceil(a_64*b_64));

のように変更すると、以下のような結果になります。

50 55.000000 56.000000
90 99.000000 100.000000
100 110.000000 111.000000
110 121.000000 122.000000
170 187.000000 188.000000
180 198.000000 199.000000
190 209.000000 210.000000
200 220.000000 221.000000
210 231.000000 232.000000
（中略）
1450 1595.000000 1596.000000
1460 1606.000000 1607.000000
1470 1617.000000 1618.000000
1480 1628.000000 1629.000000
1490 1639.000000 1640.000000

同様の結果が、CおよびPythonでも、以下のコードを使うことで再現できます。

#include <stdio.h>
#include <math.h>

int main(int argc, char const *argv[])
{
    float a_32 = 1.1, b_32;
    double a_64 = 1.1, b_64;

    for (int i = 0; i < 1500; i++) {
        b_32 = i;
        b_64 = i;
        if (ceil(a_32*b_32) != ceil(a_64*b_64)) {
            printf("%d %f %f\n", i, ceil(a_32*b_32), ceil(a_64*b_64));
        }
    }

    return 0;
}

import numpy as np
for x in range(1499):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)

    a_64 = np.array(1.1, dtype=np.float64)
    b_64 = np.array(x, dtype=np.float64)
    ceil_64 = np.ceil(a_64*b_64)
    
    if ceil_32 != ceil_64:
        print(x, ceil_32, ceil_64)

さて、先程の実行結果を見ると、基本的に10の倍数において、float32とfloat64で結果が1違う、ということがわかります。しかし、 10の倍数の中でも、計算結果が合う場合がある ことがわかります。

以下では、その原因を分析していきます。

1.1のfloat32とfloat64における表現

まず、float32とfloat64における1.1の、浮動小数型としての実際の2進数表現を出力します。

以下のJavaのコードを使います。

public class Sif {

    public static String float2string (float x) {
        return Integer.toBinaryString(Float.floatToRawIntBits(x));
    }

    public static String double2string (double x) {
        return Long.toBinaryString(Double.doubleToRawLongBits(x));
    }

    public static void main(String[] args) {
        float a_32 = 1.1f;
        double a_64 = 1.1d;
        System.out.printf("float32 (binary):\n");
        System.out.printf("1.1      : 00%s\n", float2string(a_32));
        System.out.printf("           ^~~~~~~~~=======================\n");
        System.out.printf("float64 (binary):\n");
        System.out.printf("1.1      : 00%s\n", double2string(a_64));
        System.out.printf("           ^~~~~~~~~~~~====================================================\n");
        System.out.printf("           ^: 符号部\n");
        System.out.printf("           ~: 指数部\n");
        System.out.printf("           =: 仮数部\n");
    }
}

すると、以下のような結果を得ます。

float32 (binary):
1.1      : 00111111100011001100110011001101
           ^~~~~~~~~=======================
float64 (binary):
1.1      : 0011111111110001100110011001100110011001100110011001100110011010
           ^~~~~~~~~~~~====================================================
           ^: 符号部
           ~: 指数部
           =: 仮数部

ここで、浮動小数の規格*2を見ると、float32、float64では、それぞれ後ろ23ビット、52ビットが仮数部ビットになっていることがわかります。つまり、float32, float64における仮数部ビットは、それぞれ

00011001100110011001101
0001100110011001100110011001100110011001100110011010

です。

ここで、浮動小数では、 仮数部の先頭に1をつけたものが、実際の数値としての仮数部になります。 これは、2進数では0以外の数の先頭の桁は必ず1になることから、その桁を省略しても必ず数を復元することができることに由来しています。そのため、実際の仮数部は

100011001100110011001101
10001100110011001100110011001100110011001100110011010

になります。

続いて、これらをわかりやすく扱うため10進数の整数に直します。そのために、以下のPythonコードを使います。

print(int("100011001100110011001101", 2))
print(int("10001100110011001100110011001100110011001100110011010", 2))

すると、以下の結果を得ます。

9227469
4953959590107546

一方、指数部の方を見てみると、それぞれ 01111111 と 01111111111 になっています。これらをunsigned intとして解釈して10進数に直すと127、1023になります。しかし、浮動小数の規格では、数値としての指数部の値に バイアス を加えたものを指数部ビットに入れています。float32、float64ではバイアスはそれぞれ127、1023なので、これらを差し引くと、数値としての指数部はいずれもちょうど0になることがわかります。

つまり、float32、float64では、それぞれ

$1.1 \approx 9227469 * 2^{-23} = 1.10000002384185791015625$

$1.1 \approx 4953959590107546 * 2^{-52} = 1.100000000000000088817841970012523233890533447265625$

という近似を使って1.1を表していることがわかります。ここで、 $2^{-23}$ 、 $2^{-52}$ 倍するのは、仮数部がそれぞれ23ビット、52ビットで表現されているためです。

このような誤差が出る原因は、1.1（10進表記）は10進数では循環小数にならないが、2進数では循環小数になることにあります。実際、1.1の仮数部を見てみると、00011001100…のように、000の後に1100が続いていることがわかります。浮動小数における1.1は、この循環小数を途中で打ち切ったものとして表現されているのが、この誤差を生む原因になっています。

1.1を10倍した時

さて、冒頭の計算が合わないのは、計算した範囲だと全て入力が10の倍数である場合でした。そこで、先程計算した仮数部が10倍された時の結果を見ます。そのために、仮数部の2倍、8倍、10倍を求めます。

以下のコードを使います。

a_32 = 9227469
print("a_32:   {0:30b}".format(a_32))
print()
print("a_32*2: {0:30b}".format(a_32*2))
print("a_32*8: {0:30b}".format(a_32*8))
print("a_32*10:{0:30b}".format(a_32*10))
print(" "*(11+24) + "^")
print(" "*11 + "".join(map(str, range(10)))*3)
print(" "*11 + (" "*9).join(map(str, range(3))))
print()

a_64 = 4953959590107546
print("a_64:   {0:60b}".format(a_64))
print()
print("a_64*2: {0:60b}".format(a_64*2))
print("a_64*8: {0:60b}".format(a_64*8))
print("a_64*10:{0:60b}".format(a_64*10))
print(" "*(12+53) + "^")
print(" "*12 + "".join(map(str, range(10)))*6)
print(" "*12 + (" "*9).join(map(str, range(6))))

すると、以下のようになります。

a_32:         100011001100110011001101

a_32*2:      1000110011001100110011010
a_32*8:    100011001100110011001101000
a_32*10:   101100000000000000000000010
                                   ^
           012345678901234567890123456789
           0         1         2

a_64:          10001100110011001100110011001100110011001100110011010

a_64*2:       100011001100110011001100110011001100110011001100110100
a_64*8:     10001100110011001100110011001100110011001100110011010000
a_64*10:    10110000000000000000000000000000000000000000000000000100
                                                                 ^
            012345678901234567890123456789012345678901234567890123456789
            0         1         2         3         4         5

上記において、 a_32 と a_64 が、それぞれ1.1のfloat32、float64における仮数部を整数に直したものを表しています。

まず、2進数において2倍するという操作は0を末尾に1つ付け加える事に対応し、8倍するという操作は0を末尾に3つ付け加えるという操作に対応します。そして、10倍するという操作はそれらの結果を足す事に相当します。

それらを踏まえた上で、上記のように1.1のfloat32、float64における仮数部を10倍すると、循環部分の1100が綺麗に打ち消し合い、末尾付近の桁のみ誤差が残ることがわかります。上記の結果において、 ^ で示した位置が、最終的に仮数部ビットに含まれない最後の桁になっています。末尾の誤差の位と仮数部ビットに含まれない位の位置関係を見ると、 float64の場合の方が、float32に比べて、誤差が仮数部ビットに含まれない桁に近い ことがわかります。

これについて、float64の場合により詳しい計算を行っていきます。

1.1を一般の10の倍数倍した時

さて、1.1を一般の10の倍数倍した時に、float64の仮数部の途中計算の様子を見るために、1.1のfloat64における仮数部を10の倍数倍した結果の2進数表記を見てみましょう。

以下のコードを使います。

import numpy as np
a = 4953959590107546

def ceil32_and_64_are_different(x):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)

    a_64 = np.array(1.1, dtype=np.float64)
    b_64 = np.array(x, dtype=np.float64)
    ceil_64 = np.ceil(a_64*b_64)
    
    return ceil_32 != ceil_64

for i in [10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,310,320,330,340,630,640,650,660,920,930,940,950,1270,1280,1290,1300]:
    o = ("x" if ceil32_and_64_are_different(i) else " ") + " "
    binstr = "{0:b}".format(a*i)
    print(o + "{0:4d}:{1} {2}".format(i, binstr[:54], binstr[54:]))
print(" "*(7+53) + "^")
print(" "*7 + ("".join(map(str, range(10)))*6)[:54])
print(" "*7 + (" "*9).join(map(str, range(6))))

すると、以下のような結果を得ます。

    10:101100000000000000000000000000000000000000000000000001 00
    20:101100000000000000000000000000000000000000000000000001 000
    30:100001000000000000000000000000000000000000000000000000 1100
    40:101100000000000000000000000000000000000000000000000001 0000
x   50:110111000000000000000000000000000000000000000000000001 0100
    60:100001000000000000000000000000000000000000000000000000 11000
    70:100110100000000000000000000000000000000000000000000000 11100
    80:101100000000000000000000000000000000000000000000000001 00000
x   90:110001100000000000000000000000000000000000000000000001 00100
x  100:110111000000000000000000000000000000000000000000000001 01000
x  110:111100100000000000000000000000000000000000000000000001 01100
   120:100001000000000000000000000000000000000000000000000000 110000
   130:100011110000000000000000000000000000000000000000000000 110100
   140:100110100000000000000000000000000000000000000000000000 111000
   150:101001010000000000000000000000000000000000000000000000 111100
   160:101100000000000000000000000000000000000000000000000001 000000
x  170:101110110000000000000000000000000000000000000000000001 000100
x  180:110001100000000000000000000000000000000000000000000001 001000
   310:101010101000000000000000000000000000000000000000000000 1111100
   320:101100000000000000000000000000000000000000000000000001 0000000
x  330:101101011000000000000000000000000000000000000000000001 0000100
x  340:101110110000000000000000000000000000000000000000000001 0001000
   630:101011010100000000000000000000000000000000000000000000 11111100
   640:101100000000000000000000000000000000000000000000000001 00000000
x  650:101100101100000000000000000000000000000000000000000001 00000100
x  660:101101011000000000000000000000000000000000000000000001 00001000
x  920:111111010000000000000000000000000000000000000000000001 01110000
x  930:111111111100000000000000000000000000000000000000000001 01110100
   940:100000010100000000000000000000000000000000000000000000 101111000
   950:100000101010000000000000000000000000000000000000000000 101111100
  1270:101011101010000000000000000000000000000000000000000000 111111100
  1280:101100000000000000000000000000000000000000000000000001 000000000
x 1290:101100010110000000000000000000000000000000000000000001 000000100
x 1300:101100101100000000000000000000000000000000000000000001 000001000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

ここで、 x はfloat32とfloat64で結果が合わない場合、 ^ はfloat64において仮数部に含まれない最初の桁、つまり仮数部ビットの53ビット目に相当する位を表しています。浮動小数の規格では仮数部の先頭の位の1は省略されるので、実際の位としては上から数えて54桁目に相当します。読みやすさのため、53ビット目以降をスペースで区切っています。

さて、実際の乗算の計算では、こうして得られた仮数部をfloat64表現に直すために、最後にこの仮数部の丸め込み処理を行います。仮数部は52ビットあるので、 ^ のある53ビット目で丸め込み をすることで52ビット目を決定します。（正確には、基本は四捨五入を行うが、丁度中間の値のみ切り下げを行う「最近接丸め」を使っているようです。これについては、次の節で説明します。また、仮数部の冒頭の1を省略した状態での53ビット目を表しています。）つまり、 53ビット目が1ならば52ビット目が1に、53ビット目が0ならば52ビット目が0 になるはずです。

試しに930、940の場合にどのように丸め込まれるかを見てみると、

x  930:111111111100000000000000000000000000000000000000000001 01110100
x  930:11111111110000000000000000000000000000000000000000001

   940:100000010100000000000000000000000000000000000000000000 101111000
   940:10000001010000000000000000000000000000000000000000000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

のようになります。すると、 930の場合は仮数部の一番低い位に誤差1が乗っているのに対し、940の場合は仮数部には誤差が乗っていない ことがわかります。これが、 float64において ceil(1.1 * x) の計算が合わない場合の原因 、および float32とfloat64で計算が合う場合と合わない場合の違い になっています。

四捨五入でも計算が合わない場合

さて、^ のある53ビット目で四捨五入 を行って52ビット目を決定すると書きましたが、これだけでは説明できない場合が存在します。それは、10,20,40,80,160,320,640,1280の場合で、これらの場合では53ビット目が1になっているにもかかわらず、ceilにおいて正しく結果が計算されていることがわかります。

これは、おそらく浮動小数の規格において、丸め方として 最近接丸め を使っているためであると予想されます。Wikipediaの浮動小数の規格の記事における「浮動小数点数の丸め」の項目を見ると、最近接丸めについて以下のように説明されています。

最も近くの表現できる値へ丸める。表現可能な2つの値の中間の値であったら、一番低い仮数ビット（桁）が0になるほうを採用する。これは二進での標準動作かつ十進でも推奨となっている。

ここで、これが「二進での標準動作かつ十進でも推奨となっている」と説明されていることから、浮動小数の掛け算の結果においても同様の処理が行われていると考えることができます。

さて、上記の結果において10,20,40,80,160,320,640,1280の場合のみを見てみると、仮数部の 54桁目以降が全て0 になっていることがわかります。

    10:101100000000000000000000000000000000000000000000000001 00
    20:101100000000000000000000000000000000000000000000000001 000
    40:101100000000000000000000000000000000000000000000000001 0000
    80:101100000000000000000000000000000000000000000000000001 00000
   160:101100000000000000000000000000000000000000000000000001 000000
   320:101100000000000000000000000000000000000000000000000001 0000000
   640:101100000000000000000000000000000000000000000000000001 00000000
  1280:101100000000000000000000000000000000000000000000000001 000000000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

四捨五入を行う桁である53桁目が1で、それ以降の54桁目以降が全て0ということは、これらの場合において、仮数部は 丁度四捨五入の中間の位置 にあることになります。最近接丸みの説明では、この場合は 仮数部が0となるように丸める 、つまり 切り捨てを行う と書いてあるので、53桁目が1であるにも関わらず、切り捨てを行うことで最終的な結果を得ます。

その結果、 これらの場合では仮数部の最終桁は0になる ことから、浮動小数由来の誤差が乗らず、 ceilが正しく計算される ことになります。

つまり、丸め込みによる仮数部の52桁目の決定方法は、

53桁目が0ならば、52桁目は0
53桁目が1の場合、
- 54桁目以降が全て0ならば、52桁目は0
- そうでない場合は、通常の四捨五入の通り、52桁目は1

という処理を行うことで決定できると予想されます。

この方法を用いると、上記の0から1499の場合の ceil(1.1 * x) の計算結果を全て説明することができます。

いつ誤差が発生するのか

float64における1.1の仮数部の表現を10倍した時の数の表現を見ると、

a_64*10:    10110000000000000000000000000000000000000000000000000100
                                                                 ^
            012345678901234567890123456789012345678901234567890123456789
            0         1         2         3         4         5

のように、誤差は1桁のみ存在することがわかります。そのため、1.1 * m * 10 という数の誤差の桁数は、mの桁数そのものになることがわかります。

また、float32とfloat64において結果が合わない場合の一覧を見ると、計算結果の整数部の桁が繰り上がった時に、 誤差の桁が52ビット目から押し出されることで計算が合うようになっている ことがわかります。

以上と先程の丸め込みの方法を踏まえると、float64で ceil(1.1 * x) が正しく計算される場合は、以下のようにして概ね説明できます。

x = 2ⁿ * 10 （10進表現。n=0,1,2,…）の場合は、計算が合う。
x = 2ⁿ * 10 + 10m （10進表現。m=1,2,…）を計算して行く。
この計算結果の整数部を２進数表示した時の桁数を見ていく。桁数が繰り上がるまでは、仮数部に誤差が乗っていることから ceil(1.1 * x) の計算は合わない。桁が繰り上がった後は、繰り上がったことによって浮動小数による誤差が52ビット目から押し出され、丸め込みによって誤差が無くなり、 ceil(1.1 * x) の計算が正しく行われるようになる。
2ⁿ * 10 + 10m == 2ⁿ⁺¹ * 10 （10進表現）に到達したら、1.に戻る。

float32の場合

最後に、同様の計算をfloat32において行った場合について見ます。

以下のコードを使います。

import numpy as np
a = 9227469

for i in [10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,310,320,330,340,630,640,650,660,920,930,940,950,1270,1280,1290,1300]:
    binstr = "{0:b}".format(a*i)
    print("  {0:4d}:{1} {2}".format(i, binstr[:25], binstr[25:]))
print(" "*(7+24) + "^")
print(" "*7 + ("".join(map(str, range(10)))*6)[:25])
print(" "*7 + (" "*9).join(map(str, range(3))))

すると、以下のような結果を得ます。

    10:1011000000000000000000000 10
    20:1011000000000000000000000 100
    30:1000010000000000000000000 0110
    40:1011000000000000000000000 1000
    50:1101110000000000000000000 1010
    60:1000010000000000000000000 01100
    70:1001101000000000000000000 01110
    80:1011000000000000000000000 10000
    90:1100011000000000000000000 10010
   100:1101110000000000000000000 10100
   110:1111001000000000000000000 10110
   120:1000010000000000000000000 011000
   130:1000111100000000000000000 011010
   140:1001101000000000000000000 011100
   150:1010010100000000000000000 011110
   160:1011000000000000000000000 100000
   170:1011101100000000000000000 100010
   180:1100011000000000000000000 100100
   310:1010101010000000000000000 0111110
   320:1011000000000000000000000 1000000
   330:1011010110000000000000000 1000010
   340:1011101100000000000000000 1000100
   630:1010110101000000000000000 01111110
   640:1011000000000000000000000 10000000
   650:1011001011000000000000000 10000010
   660:1011010110000000000000000 10000100
   920:1111110100000000000000000 10111000
   930:1111111111000000000000000 10111010
   940:1000000101000000000000000 010111100
   950:1000001010100000000000000 010111110
  1270:1010111010100000000000000 011111110
  1280:1011000000000000000000000 100000000
  1290:1011000101100000000000000 100000010
  1300:1011001011000000000000000 100000100
                               ^
       0123456789012345678901234
       0         1         2

float32では仮数部は23ビットあるので、 ^ で示した24ビット目（冒頭の1を省略して数えた場合）に対して丸め込みを行った結果が最終的な計算結果となります。しかし、上記の場合では 23ビット目はいずれも0 なので、 これらの場合では常にceilの計算が正しく行われる ことがわかります。

実際、他の場合でもこれが成り立つか調べるため、

import numpy as np
for x in range(200000):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)
    ceil = x + x // 10 + (1 if x % 10 > 0 else 0)

    if ceil_32 != ceil:
        print(x, ceil_32, ceil)

のようなコードを実行すると、0から200000未満の全ての x において、float32では ceil(1.1 * x) が正しく計算されていることがわかります。

これは、この記事の冒頭付近の「1.1のfloat32とfloat64における表現」の節で見た通り、float32における1.1の表現では、循環小数を丸め上げたことによる仮数部の誤差が、最終的に丸め上げる位に対して小さくなっているため起きていると考えられます。

実応用：ラブライブ！スクフェスにおけるスコア計算

この現象はスクフェスにおけるスコア計算という実用的な場面において発生しています。

応援タップSCOREアップの計算方法はまだ解が見つかっていないけど、とりあえず回復のパラメータアップの方の検証を進め、体力41～54まで判明。なお、応援タップSCOREアップは、SCOREの繰り下げ前の小数にアップ率を掛けて繰り上げると良さそうだけど、1つだけ合わない実測値ががががが… pic.twitter.com/9bsTnL18Nm
— まぐお＠ﾌﾗｲﾝｸﾞﾊﾞｰﾄﾞ=ｴﾘｺﾇｽⅡ世 (@magu_maki) 2019年2月22日

この@magu_makiさんの解析の概要を説明すると、

スクフェスにおいて、１回タップした時のスコアの理論式の解析を行っている
理論式の中に、ある整数（ベースとなるタップスコア）を1.1倍した上で切り上げをする という操作がある
この操作が、入力が820の時は計算が合わない（上記ツイートの3枚目の画像の赤枠の数値）
しかし、同じ10の倍数である970では計算が合っている（同画像の赤枠の1つ右、2つ下の枠の数値）

分析の結果、以下の事がわかりました。

冒頭で列挙した、 ceil(1.1 * x) がfloat32とfloat64で計算が合わない場合を見ると、820では計算が合わないが、970では合うことがわかる
計算が行われたExcelではfloat32で計算が行われている
Excelのマクロ、つまりVBAではDouble型、つまりfloat64での計算が使用可能
VBAを使って切り上げ処理を実装し、切り上げ処理のみfloat64で計算するようにした所、今まで計算が合っていなかった場所を含め、他の場所でも計算が合うようになった

結果的に、スクフェスのゲーム中における実際のスコアの測定値を説明するにはfloat64で切り上げ計算を行う必要があることがわかりました。

このようなビット単位で発生している現象が、2019年の今、ファミコンなどではなく、実際にiOSやAndroidで配布されているアプリの挙動を再現するという実用的な場面で起きているということに、今回は非常に驚かされました。

結論

ceil(1.1 * x) の計算が、float32では正しく計算されるが、float64では正しく計算されない場合があることを確認した。この現象がPython、Java、Cで発生することを確認した。
この現象は基本的に入力が10の倍数の時に発生するようである、と予想できることを見た。しかし、10の倍数の中でも正しく計算される場合が存在することを確認した。
10の倍数の中でも正しく計算される場合がなぜ存在するかについて、float32とfloat64における 1.1 のデータ表現における仮数部の計算を分析することで、詳しい説明を行った。特に、float64で x が10,20,40,80,160,320,640,1280の場合は、浮動小数の最近接丸めを考慮する必要がある可能性があることを指摘した。
float32の場合で同様の計算を行うと、0から200000未満の整数の入力では、浮動小数由来の誤差が発生せず、 ceil(1.1 * x) が正しく計算されることを確認した。
以上のような浮動小数の誤差を考慮して計算を行った結果、@magu_maki によるスクフェスのスコア計算の解析において唯一説明できなかったデータ点が説明できるようになり、実際に観測されたデータが全て説明できるようになった。

*2:https://ja.wikipedia.org/wiki/IEEE_754 を参照

2019-02-02

うでまえクラスのこと

スクフェスを初めて通算ログイン日数実に157日目、ついにロイヤルエキスパートの称号をゲットした。

エキスパート+で初めて課題に詰まった時は、当時の現状で全く行き詰まったと思った。そこからどれだけ練習すれば先に進めるか全くわからなかった。でもtwitterを見ると、エキスパート++、そしてロイヤルエキスパートも、実際にクリアしている人がいた。もしそうならば、自分もできるはずだと思った。

そこで、他の人の動画を前より注意深く見ると、人差し指押しにしていたり、10速にしていたりと、うでまえクラスを始める前は深く注目していなかった部分に気付いた。全く行き詰まっていたので、「これを試したら上手くなるかな…？」と思って色々試した。

すると、それまで全く行き詰まっていたと思っていたのに、新たな課題がクリアできるようになることがあった。勿論失敗した試みもあったけど、実りがある時は一気に実力が上がった。うでまえクラスをやっていた中でもこれはかなり新鮮で面白い体験だった。これは面白い体験であると同時に、そのような小さな工夫だけで実力が一気に上がるくらい、未熟さが大きいということだった。自分の実力が俯瞰できない時がある、ということでもあるかもしれない。

そして、特に週に1回きりのみ挑戦可能なフルコンボ課題の本番では、運動神経がかなり研ぎ澄まされ緊張感に溢れた。曲の終盤は特に緊張して、そしてやりきった後は物凄い達成感と開放感に溢れた。これ程のスリルを味わわせてくれるものはなかなか無いくらいの面白い体験だった。

「リズムに合わせて画面をタッチ！」だけでルールは説明できるのに、大きな目標を目指そうと思った途端、ここまで総合的な身体能力が求められて、緊張を味わわせてくれて、そして目標を達成した時に感動させてくれるのは凄いと思った。

もちろんまだ達成できていない事はかなり多い。例えば日替わり楽曲のMasterのフルコンはまだできていない。そこに辿り着くまでにもまた色々な面白さがあると思う。

2018-09-30

ラブライブ！1期7話の穂乃果の数学のテストはかなり作り込まれている！

ラブライブ！1期7話では、μ’ｓのメンバーはラブライブへの出場の許可をもらうため、学校の理事長に相談します。その際、出場は許可するが、μ’ｓのメンバーが全員次の期末試験で赤点を取らないという条件を課され、μ’ｓのメンバーは数学の試験のために奮闘します。

この話の最後では、穂乃果の受けた実際の数学のテストを見ることが出来ます。この数学のテストと穂乃果の解答が非常に上手く作られていて（制作の方々に敬意を贈ります！）、この時の穂乃果の解答から、穂乃果の絶妙な数学の理解具合や試験中の悩み具合を見ることが出来ます。

こちらの記事で、ラブライブ！1期7話の穂乃果の数学のテストが文字起こしされています：

ameblo.jp

以下、この記事では、問題文および穂乃果の答案は、こちらの記事で文字起こしされた問題および穂乃果の答案をベースに分析をします。

問1-(5): 合っているけど…！穂乃果の深い数学への理解、そして奮闘

まず面白いのは、問1-(5)です。ここからは、穂乃果の数学の理解具合や、試験中の悩み具合を伺うことができます。

問1-(5)は次のような問題になっています：

$ \frac{1}{2- \sqrt{3}} $の整数部分を$ a $、[小]¹数の部分を$ b $とするとき

(a) $ a = \fbox{ ? } $ (b) $ b = \fbox{ ? } $ (c) $ {a^2 + 4b + b^2 = \fbox{ ? }} $

この問題に対する、穂乃果の解答は以下のようになっています。

(a): $3$ （正解）
(b): $732$ （不正解）
(c) (空欄) （不正解）

この解答において、(b)の「732」が穂乃果の絶妙な数学の理解具合、そして(c)から穂乃果の悩み具合を読み取ることができるんです！

(a)の解答

まずは、(a)の答えを導いてみましょう。

$ \frac 1 {2 - \sqrt{3}} $ の分母を有理化すると、

$ {\displaystyle \frac 1 {2 - \sqrt{3}} \ = \frac {2 + \sqrt 3} {2^2 - \sqrt{3}^2} \ = \frac {2 + \sqrt 3} {4-3} \ = 2 + \sqrt 3 }$

となります。これを小数表現に直すと

$ 2 + \sqrt 3 = 2 + 1.73205\cdots = 3.73205\cdots $

となるので、まず(a)は3で穂乃果の正解となります。

(b)の解答は合っている!?

さて、(b)は小数部分を答える問題ですが、よく見ると四捨五入さえすれば答えは732で合っています！！

それにも関わらず、不正解になっています。これは、この問の設問は

[小] ²数の部分を$ b $とするとき

となっていることが理由であると考えられます。このように小数点以下第何位で四捨五入するかについて指定がない時は、$ \sqrt3 - 1 $ のような形式が正解になります³。このような理由から不正解となっていたと考えるのが最も有力でしょう。

ここからわかるのは、穂乃果は$ \sqrt 3 = 1.73205\cdots $であることを確実に理解していて、 更に $ 3.732\cdots $ という点まで試験中に辿り着いていた ということなんです!!更に、そもそもこの式に辿り着くには、(a)を導出するための変形が行える必要があります。つまり、穂乃果は

分母の有理化
ルートの計算
$ \sqrt 3 = 1.73205\cdots $であるということ

を理解していたということになります。この時点で、かなり本質的にこの問題を解くための数学を理解していると言えるでしょう。

問6

この問題は放物線の平行移動に関する問題です。問題文にある放物線の式は、一見正解とは全く異なるものとなっています。実際、xとyで-2と3平行移動、入れ替えて3と-2平行移動させても、一致しません。

「これって、どんな放物線の形になるんだろう…」と思い、Google先生にグラフを描画してもらうと…………!!

f:id:snowharinchan:20181001000353p:plain

あああああーーーッ！！！これ！！x=-2と3を通る放物線じゃん！！

穂乃果は、$ (x+2)(x-3) $ という導出を行ったに違いない。

では、どうしてそのような答えに辿り着いたか？これは、２次方程式が $ (x-\alpha)(x-\beta)=0$ の形に変形できることを知っているからこそできる技なのです!!

問5: メタなネタ

問5は穂乃果のことに関してわかるわけでは無いのですが、メタなネタになっています。

その問題文はこうなっています：

職場の制作をいくつかの班に分けるのに、1班3人にすると、12人が入れない。また、1班7人にすると、最後の班だけは7人に満たないという。班の数および制作の人数を求めよ。

班の数　　　班，　進行の人数　　　人

アニメ制作のネタが入っています。

この問題は、穂乃果はいずれも空欄で答えて点を貰っていません。そのことと問題文を合わせると、これは一見解けない問題のように見えますが、実はちゃんと解くことができます。

その解答は以下のようになります。

班の数を$x$、制作の人数を$y$と置きます。 1班3人にすると12人が入れないことから、

$3x+12=y$

が言えます。また、1班7人にした場合、最後の班で不足している人数を$z$と置くと、

$7x-z=y, 0<z<7$

が言えます。これらを連立させると、

$x = z/4+3, y=3(z+28)/4, 0<z<7$

が言えます。$x,y,z$は整数であることから$z=4$が言えて、したがって$x=4, y=24$となります。

$z$を導入せずより高校数学風に解くならば、2番目の式を$ 7(x-1) < y < 7x $ として、$x$が整数であるという条件を使うのが良いでしょう。⁴

最後に

以上、ラブライブ！における数学の解答に関する考察でした！

冒頭で取り上げた記事で文字起こしされた結果では「少数」となっていました。↩
[1]と同様。↩
詳しくは、例えば整数部分と小数部分の意味と例題 | 高校数学の美しい物語などを参照。このような時、小数部分の0.732…は無限小数のため、「732」では省略してしまっていて正確に表すことができません。$ \sqrt 3 = 1.732\cdots $ なので、ここから$1$を引くと$ \sqrt 3 - 1 = 0.732\cdots $のように、小数部分だけを得ることができるので、これが正解となるのです。↩
参考: http://www.wolframalpha.com/input/?i=7(x-1)%3C3x%2B12%3C7x ↩

導出

整数部分が2以上の場合

マルコフ連鎖の話との繋がり

パターンCにおけるカバー率の近似式 - 改良版

導出

理論の検証

結論

免責

この記事の内容

目次

導出

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

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

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

マルコフ過程による表現

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

状態訪問数の期待値

定常状態による近似

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

最終的な理論式

理論の検証

外れ値

議論

近似の解釈

分散の評価

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

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

結論

訂正履歴

免責

概要

背景

試みられた計算: 10の倍数が関係している？

深まる謎: 丸め込み誤差のはずなのに、10の倍数でも正しく丸め込まれる場合がある？

打開策: 「float32では計算が合うのに、float64で計算すると計算が合わない」事がある

発生する理由：末尾1ビットのみがずれる場合がある

ゲーム的なコメント

結論、感想

概要

目次

計算が合わない場合の列挙

1.1のfloat32とfloat64における表現

1.1を10倍した時

1.1を一般の10の倍数倍した時

四捨五入でも計算が合わない場合

いつ誤差が発生するのか

float32の場合

実応用：ラブライブ！スクフェスにおけるスコア計算

結論

問1-(5): 合っているけど…！穂乃果の深い数学への理解、そして奮闘

(a)の解答

(b)の解答は合っている!?

問6

問5: メタなネタ

最後に