DanceDanceRevolution1 は KONAMI が開発している音楽ゲームのひとつです。DanceDanceRevolution には譜面ごとに難易度が振られていて2,その譜面をプレイするのがどれほど難しいかを示しています。
それとは別に,グルーブレーダーという仕組みがあって,その譜面の傾向を示します。各要素は以下の 5 通りです。
- STREAM
- 平均密度。曲中のオブジェ数が多いほど高くなる。
- VOLTAGE
- 最高密度。最もオブジェが多い 4 拍のオブジェ数が多いほど高くなる。
- AIR
- 跳躍頻度。同時踏みや踏んではいけないオブジェが多いほど高くなる。
- FREEZE
- 拘束度。どこかのパネルを踏み続ける拍が多いほど高くなる。
- CHAOS
- 変則度。細かいリズムや変速が多いほど高くなる。
グルーブレーダーの数値は譜面自体から厳密に計算されて求めることができます。計算式は公開されていませんが,有志のプレイヤーによってかなりの精度で明らかにされています。
一方で難易度の数値は制作側によって人為的に決定されています。そのため,バージョンアップのタイミングなどで難易度の見直しなどが行われる場合があります。
では,グルーブレーダーの数値から難易度を推定することは可能でしょうか。やってみましょう。
環境
- Python3 + JupyterLab
- Matplotlib
- NumPy
- Pandas
- PyCM
- SciPy
- Seaborn
前準備
インポートなどをしておきます。
import math
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from pycm import ConfusionMatrix
from scipy.optimize import minimize, differential_evolution, Bounds
def ustd_coefficient(n):
try:
return math.sqrt(n / 2) * math.gamma((n - 1) / 2) / math.gamma(n / 2)
except OverflowError:
return math.sqrt(n / (n - 1.5))
def std_u(a):
return np.std(a) * ustd_coefficient(len(a))
oo = math.inf
各譜面のデータを読み込みます。データとして,2020-03-19 当時の BEMANI wiki から旧曲と DanceDanceRevolution A20 の新曲のデータを CSV にしておきました。こちらに配置しておきます。今回は旧曲をフィッティングに使用する訓練データ,新曲を評価データとします。では,読み込んで DataFrame にしていきます。
old_csv = Path('./old.csv')
new_csv = Path('./new.csv')
train_df = pd.read_csv(old_csv)
test_df = pd.read_csv(new_csv)
display(train_df)
display(test_df)
| VERSION | MUSIC | SEQUENCE | LEVEL | STREAM | VOLTAGE | AIR | FREEZE | CHAOS | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | DanceDanceRevolution A | 愛言葉 | BEGINNER | 3 | 21 | 22 | 7 | 26 | 0 |
| 1 | DanceDanceRevolution A | 愛言葉 | BASIC | 5 | 34 | 22 | 18 | 26 | 0 |
| 2 | DanceDanceRevolution A | 愛言葉 | DIFFICULT | 7 | 43 | 34 | 23 | 26 | 7 |
| 3 | DanceDanceRevolution A | 愛言葉 | EXPERT | 11 | 63 | 45 | 21 | 25 | 28 |
| 4 | DanceDanceRevolution A | 天ノ弱 | BEGINNER | 3 | 20 | 25 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3390 | DanceDanceRevolution 1st | PARANOiA | EXPERT | 11 | 67 | 52 | 25 | 0 | 17 |
| 3391 | DanceDanceRevolution 1st | TRIP MACHINE | BEGINNER | 3 | 25 | 26 | 5 | 0 | 0 |
| 3392 | DanceDanceRevolution 1st | TRIP MACHINE | BASIC | 8 | 47 | 40 | 14 | 0 | 4 |
| 3393 | DanceDanceRevolution 1st | TRIP MACHINE | DIFFICULT | 9 | 52 | 40 | 30 | 0 | 7 |
| 3394 | DanceDanceRevolution 1st | TRIP MACHINE | EXPERT | 10 | 56 | 53 | 36 | 0 | 12 |
3395 rows × 9 columns
| VERSION | MUSIC | SEQUENCE | LEVEL | STREAM | VOLTAGE | AIR | FREEZE | CHAOS | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | DanceDanceRevolution A20 | おーまい！らぶりー！すうぃーてぃ！だーりん！ | BEGINNER | 3 | 18 | 21 | 5 | 16 | 0 |
| 1 | DanceDanceRevolution A20 | おーまい！らぶりー！すうぃーてぃ！だーりん！ | BASIC | 7 | 37 | 28 | 18 | 39 | 0 |
| 2 | DanceDanceRevolution A20 | おーまい！らぶりー！すうぃーてぃ！だーりん！ | DIFFICULT | 12 | 60 | 56 | 54 | 55 | 21 |
| 3 | DanceDanceRevolution A20 | おーまい！らぶりー！すうぃーてぃ！だーりん！ | EXPERT | 15 | 95 | 99 | 30 | 25 | 100 |
| 4 | DanceDanceRevolution A20 | 革命パッショネイト | BEGINNER | 3 | 16 | 16 | 1 | 35 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 380 | DanceDanceRevolution A20 | 50th Memorial Songs -The BEMANI History- | EXPERT | 13 | 63 | 79 | 14 | 62 | 63 |
| 381 | DanceDanceRevolution A20 | 50th Memorial Songs -二人の時 ～under the cherry bl... | BEGINNER | 3 | 17 | 20 | 3 | 46 | 0 |
| 382 | DanceDanceRevolution A20 | 50th Memorial Songs -二人の時 ～under the cherry bl... | BASIC | 7 | 40 | 33 | 36 | 29 | 0 |
| 383 | DanceDanceRevolution A20 | 50th Memorial Songs -二人の時 ～under the cherry bl... | DIFFICULT | 9 | 50 | 46 | 47 | 3 | 6 |
| 384 | DanceDanceRevolution A20 | 50th Memorial Songs -二人の時 ～under the cherry bl... | EXPERT | 12 | 73 | 60 | 60 | 15 | 32 |
385 rows × 9 columns
さらに,各グルーブレーダーの数値を標準化していきます。訓練データの平均が 0 で標準偏差が 1 になるようにし,同様の操作を評価データにも行います。
grs = ['STREAM', 'VOLTAGE', 'AIR', 'FREEZE', 'CHAOS']
sgrs = ['S_{}'.format(gr) for gr in grs]
m = {}
s = {}
for gr, sgr in zip(grs, sgrs):
v = train_df.loc[:, gr].values
v_t = test_df.loc[:, gr].values
m[gr] = np.mean(v)
s[gr] = std_u(v)
train_df[sgr] = (v - m[gr]) / s[gr]
test_df[sgr] = (v_t - m[gr]) / s[gr]
display(train_df.loc[:, sgrs])
display(test_df.loc[:, sgrs])
| S_STREAM | S_VOLTAGE | S_AIR | S_FREEZE | S_CHAOS | |
|---|---|---|---|---|---|
| 0 | -0.981448 | -0.838977 | -0.636332 | 0.056063 | -0.661167 |
| 1 | -0.534364 | -0.838977 | -0.160513 | 0.056063 | -0.661167 |
| 2 | -0.224844 | -0.405051 | 0.055768 | 0.056063 | -0.441192 |
| 3 | 0.462978 | -0.007285 | -0.030744 | 0.014296 | 0.218735 |
| 4 | -1.015839 | -0.730495 | -0.939125 | -1.029883 | -0.661167 |
| ... | ... | ... | ... | ... | ... |
| 3390 | 0.600542 | 0.245838 | 0.142280 | -1.029883 | -0.126941 |
| 3391 | -0.843883 | -0.694335 | -0.722844 | -1.029883 | -0.661167 |
| 3392 | -0.087279 | -0.188088 | -0.333538 | -1.029883 | -0.535467 |
| 3393 | 0.084676 | -0.188088 | 0.358562 | -1.029883 | -0.441192 |
| 3394 | 0.222240 | 0.281999 | 0.618099 | -1.029883 | -0.284066 |
3395 rows × 5 columns
| S_STREAM | S_VOLTAGE | S_AIR | S_FREEZE | S_CHAOS | |
|---|---|---|---|---|---|
| 0 | -1.08462 | -0.87514 | -0.72284 | -0.36161 | -0.66117 |
| 1 | -0.43119 | -0.62201 | -0.16051 | 0.599036 | -0.66117 |
| 2 | 0.359805 | 0.39048 | 1.396711 | 1.26731 | -0.00124 |
| 3 | 1.563493 | 1.945381 | 0.358562 | 0.014296 | 2.481343 |
| 4 | -1.1534 | -1.05594 | -0.89587 | 0.431967 | -0.66117 |
| ... | ... | ... | ... | ... | ... |
| 380 | 0.462978 | 1.222171 | -0.33354 | 1.55968 | 1.318614 |
| 381 | -1.11901 | -0.9113 | -0.80936 | 0.891406 | -0.66117 |
| 382 | -0.32802 | -0.44121 | 0.618099 | 0.181364 | -0.66117 |
| 383 | 0.015894 | 0.028875 | 1.093917 | -0.90458 | -0.47262 |
| 384 | 0.806889 | 0.535122 | 1.656248 | -0.40338 | 0.344436 |
385 rows × 5 columns
そして,各譜面のグルーブレーダーを示す 2 階テンソルと,各譜面の難易度を示す 1 階テンソルを抜き出しておきます。
train_sgr_arr = train_df.loc[:, sgrs].values
test_sgr_arr = test_df.loc[:, sgrs].values
train_level_arr = train_df.loc[:, 'LEVEL'].values
test_level_arr = test_df.loc[:, 'LEVEL'].values
最小二乗法による重回帰分析
重回帰分析は以下のような考え方によります。
説明変数群 $x_n$ と目的変数 $y$ があります。今回の場合説明変数とはグルーブレーダーの各値のことです。目的変数とは難易度です。このとき,
$$
y' = k_0 + k_1x_1 + k_2x_2 + \cdots + k_nx_n
$$
なる係数群 $k_n$ と $y’$ を考えたとき,$m$ 個のデータの二乗誤差 $e^2 := \sum_{i = 1}^m\left(y'_i - y_i\right)^2$ が最も小さくなるような $k_n$ を探索していきます。今回は,このような最適化問題を SciPy を使用した差分進化法で求めていきます。
まず,最小化したい関数を定義します。先程の $e^2$ です。
def hadprosum(a, b):
return (a * b).sum(axis=1)
def estimate(x, sgr_arr):
x_const = x[0]
x_coef = x[1:]
return hadprosum(sgr_arr, x_coef) + x_const
def sqerr(x):
est = estimate(x, train_sgr_arr)
return ((est - train_level_arr) ** 2).sum()
これを SciPy のdifferential_evolution関数に与えます。探索範囲に関しては色々試しながら十分と思われる範囲をズドンで与えています。
bounds = Bounds([0.] * 6, [10.] * 6)
result = differential_evolution(sqerr, bounds, seed=300)
print(result)
fun: 5170.056057917698
jac: array([-0.00236469, 0.14933903, 0.15834303, 0.07094059, 0.01737135,
0.1551598 ])
message: 'Optimization terminated successfully.'
nfev: 3546
nit: 37
success: True
x: array([8.04447683, 2.64586828, 0.58686288, 0.42785461, 0.45934494,
0.4635763 ])
この結果をみると,最も影響を与えているのは STREAM で,そのあとに VOLTAGE,CHAOS,FREEZE,AIR と続いているようです。
では,実際に得られたパラメータを使って評価を行ってみます。
まずは予測用の関数を定義します。この関数にパラメータとグルーブレーダーのテンソルを与えると予測された難易度のテンソルが返ります。
def pred1(x, sgr_arr):
est = estimate(x, sgr_arr).clip(1., 19.)
return np.round(est).astype(int)
この関数の返り値と実際の難易度を,PyCM のConfusionMatrixに与えることで,混同行列オブジェクトが生成されます。これのプロパティにアクセスし,正解率とマクロ F 値を求めていきます。
train_pred1_arr = pred1(result.x, train_sgr_arr)
test_pred1_arr = pred1(result.x, test_sgr_arr)
train_cm1 = ConfusionMatrix(train_level_arr, train_pred1_arr)
test_cm1 = ConfusionMatrix(test_level_arr, test_pred1_arr)
print('====================')
print('Train Score')
print(' Accuracy: {}'.format(train_cm1.Overall_ACC))
print(' Fmeasure: {}'.format(train_cm1.F1_Macro))
print('====================')
print('Test Score')
print(' Accuracy: {}'.format(test_cm1.Overall_ACC))
print(' Fmeasure: {}'.format(test_cm1.F1_Macro))
print('====================')
====================
Train Score
Accuracy: 0.33431516936671574
Fmeasure: 0.2785969345790368
====================
Test Score
Accuracy: 0.3142857142857143
Fmeasure: 0.24415916194348555
====================
正解率は 31.4% となりました。これは期待していたよりもだいぶ低い結果です。混同行列を Seaborn でヒートマップ化して見てみます。
plt.figure(figsize=(10, 10), dpi=200)
sns.heatmap(pd.DataFrame(test_cm1.table), annot=True, square=True, cmap='Blues')
plt.show()
横軸が実際の難易度,縦軸が予測された難易度です。これを見ると,低い難易度の曲は高く,一定以上の難易度に関しては低く評価され,さらに高い難易度になると一転して過大評価しているようです。ヒートマップは直線ではなく内側に曲がって弓なりの形状を描いています。
最尤推定法による (順序) ロジスティック回帰分析
続いてロジスティック回帰分析を試してみます。ロジスティック回帰分析はなんらかの確率を目的変数にとる分析に向いています。以下のような式を考えます。
$$
y' = \frac{1}{1 + \exp\left[-\left(k_0 + k_1x_1 + k_2x_2 + \cdots + k_nx_n\right)\right]}
$$
この $y'$ は 0 から 1 の値をとります。$m$ 個のデータで目的変数 $y$ が 0 か 1 かで定まっているとき,尤度 $l = \prod_{i = 1}^m y_iy_i' + (1 - y_i)(1 - y_i')$ を最大化することを考えます。数式で書くと分かりづらいですが,$y'$ は 陽性である確率であり,$y$ が 1 すなわち陽性であればそのまま $y'$ を,$y$ が 0 すなわち陰性であれば陰性の確率つまり 1 から $y'$ を除いた値を使用するということです。$(0, 1)$ の区間となる $y’$ をどんどんかけていけばその値は小さくなりすぎて計算機では扱いづらくなるため,対数尤度 $\log l = \sum_{i = 1}^m \log \left[y_iy_i' + (1 - y_i)(1 - y_i')\right]$ を最大化するのが普通です。
今回は陰性か陽性かではなく,順序のあるクラスのどこに属するかを扱います。こういう場合定数項 $k_0$ のみが異なる複数のロジスティック曲線を想定し,それぞれ “難易度 2 以上である確率”, “難易度 3 以上である確率”, ... と考えます。“難易度 2 である確率” は “難易度 2 以上である確率” から “難易度 3 以上である確率” を引いたものですから,これによって尤度を求めることができます。
まず,計算のために難易度をワンホット形式の 2 階テンソルに変換しておきます。
train_level_onehot_arr = np.zeros(train_level_arr.shape + (19,))
for i, l in np.ndenumerate(train_level_arr):
train_level_onehot_arr[i, l - 1] = 1.
つづいて最小化すべき関数を与えます。最小化なので上記の対数尤度にマイナスをつけたものを定義します。
def upperscore(x, sgr_arr):
x_const = np.append(np.append(oo, x[:18].copy()), -oo) # 1以上の確率1と20以上の確率0のために両端に無限大を挿入
x_coef = x[18:]
var = np.asarray([hadprosum(sgr_arr, x_coef)]).T
cons = np.asarray([x_const])
return 1 / (1 + np.exp(-(var + cons)))
def score(x, sgr_arr):
us = upperscore(x, sgr_arr)
us_2 = np.roll(us, -1)
return np.delete(us - us_2, -1, axis=1) # ずらして引き,末尾を削除して各難易度の確率を得る
def mloglh(x):
sc = score(x, train_sgr_arr)
ret = -(np.log((sc * train_level_onehot_arr).sum(axis=1).clip(1e-323, oo)).sum())
return ret
探索を行います。先程に比べてかなり時間がかかるので注意です。
bounds = Bounds([-60.] * 18 + [0] * 5, [20.] * 18 + [10] * 5)
result = differential_evolution(mloglh, bounds, seed=300)
print(result)
fun: 4116.792196474322
jac: array([ 0.00272848, 0.00636646, -0.00090949, 0.00327418, -0.00563887,
-0.00291038, -0.00509317, 0.00045475, 0.00800355, 0.00536602,
-0.00673026, 0.00536602, 0.00782165, -0.01209628, 0.00154614,
-0.0003638 , 0.00218279, 0.00582077, 0.04783942, 0.03237801,
0.01400622, 0.00682121, 0.03601599])
message: 'Optimization terminated successfully.'
nfev: 218922
nit: 625
success: True
x: array([ 14.33053717, 12.20158703, 9.97549255, 8.1718939 ,
6.36190483, 4.58724228, 2.61478521, 0.66474105,
-1.46625252, -3.60065138, -6.27127806, -9.65032254,
-14.06390123, -18.287351 , -23.44011235, -28.39033479,
-32.35825176, -43.38390248, 6.13059504, 2.01974223,
0.64631137, 0.67555403, 2.44873606])
今回は STREAM > CHAOS > VOLTAGE > FREEZE > AIR となっており,CHAOS のほうが大きく影響しているのがわかります。
こちらも実際に評価を行ってみましょう。
def pred2(x, sgr_arr):
sc = score(x, sgr_arr)
return np.argmax(sc, axis=1) + 1
train_pred2_arr = pred2(result.x, train_sgr_arr)
test_pred2_arr = pred2(result.x, test_sgr_arr)
train_cm2 = ConfusionMatrix(train_level_arr, train_pred2_arr)
test_cm2 = ConfusionMatrix(test_level_arr, test_pred2_arr)
print('====================')
print('Train Score')
print(' Accuracy: {}'.format(train_cm2.Overall_ACC))
print(' Fmeasure: {}'.format(train_cm2.F1_Macro))
print('====================')
print('Test Score')
print(' Accuracy: {}'.format(test_cm2.Overall_ACC))
print(' Fmeasure: {}'.format(test_cm2.F1_Macro))
print('====================')
====================
Train Score
Accuracy: 0.4960235640648012
Fmeasure: 0.48246495009640167
====================
Test Score
Accuracy: 0.5454545454545454
Fmeasure: 0.5121482282311358
====================
今度は 54.5% の正解率となりました。先程よりはかなり改善されていますが,まだまだです。
plt.figure(figsize=(10, 10), dpi=200)
sns.heatmap(pd.DataFrame(test_cm2.table), annot=True, square=True, cmap='Blues')
plt.show()
こちらはおおむね直線上に出ていますが,低難易度と高難易度がやはりイマイチです。
結論
結論としては,“それっぽい値を得ることはできるが,実用できる域には至らない” といったところでしょうか。今回の試みは StepMania の譜面に難易度をつける参考にできればなと思って行いましたが,やはり最終的にはプレイして調整する必要がありそうですね。
それともうひとつ,ロジスティック回帰では当然各定数項が大小関係の制約を受けるわけですが,今回のコードでは乱数によっては制約通りにいかず,異常な値で収束判定となる可能性があります3。SciPy の最適化関数では不等式での制約を与えることができるのですが,どうも私のほうでは未知の形の制約が渡されたみたいなエラーが出てうまくいきませんでした。探索時間も無駄になるのでこちら解決できた方がいれば教授願いたいです。