DanceDanceRevolution¹ は KONAMI が開発している音楽ゲームのひとつです。DanceDanceRevolution には譜面ごとに難易度が振られていて²,その譜面をプレイするのがどれほど難しいかを示しています。

それとは別に,グルーブレーダーという仕組みがあって,その譜面の傾向を示します。各要素は以下の 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 の譜面に難易度をつける参考にできればなと思って行いましたが,やはり最終的にはプレイして調整する必要がありそうですね。

それともうひとつ,ロジスティック回帰では当然各定数項が大小関係の制約を受けるわけですが,今回のコードでは乱数によっては制約通りにいかず,異常な値で収束判定となる可能性があります³。SciPy の最適化関数では不等式での制約を与えることができるのですが,どうも私のほうでは未知の形の制約が渡されたみたいなエラーが出てうまくいきませんでした。探索時間も無駄になるのでこちら解決できた方がいれば教授願いたいです。

長くて省略したくなるが Qiita で省略するとメモリについて調べたい人の邪魔になりそうな気がしたので省略は行わないこととする。 ↩
かつて足のアイコンで示されていたので “足 16” などと表記される。 ↩
実際範囲を調整する段階で一度そのような現象に見舞われた。 ↩