[PYTHON] [DanceDanceRevolution] Is it possible to predict the difficulty level (foot) from the value of the groove radar?

Groove radar numbers can be calculated and calculated exactly from the score itself. The formula has not been published, but it has been revealed with considerable accuracy by volunteer players.

On the other hand, the difficulty value is artificially determined by the production side. Therefore, the difficulty level may be reviewed at the timing of version upgrade.

Then, is it possible to estimate the difficulty level from the numerical value of the groove radar? let's do it.

[^ DDR]: It's long and I want to omit it, but if I omit it with Qiita, I feel that it will be an obstacle for people who want to check memory, so I will not omit it. [^ level]: Since it was once indicated by a foot icon, it is written as “foot 16”.

environment

• Python3 + JupyterLab
  • Matplotlib
  • NumPy
  • Pandas
  • PyCM
  • SciPy
  • Seaborn

Preparation

Import etc.

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

Read the data of each musical score. As data, the data of the old song and the new song of DanceDanceRevolution A20 from the BEMANI wiki at that time was set to CSV. Place it at here. This time, the training data that uses the old song for fitting and the new song will be the evaluation data. Now, let's read it and make it a 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)

3395 rows × 9 columns

385 rows × 9 columns

Furthermore, we will standardize the numerical values of each groove radar. Make sure that the average of the training data is 0 and the standard deviation is 1, and the same operation is performed for the evaluation data.

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])

3395 rows × 5 columns

385 rows × 5 columns

Then, extract the 2nd floor tensor that shows the groove radar of each musical score and the 1st floor tensor that shows the difficulty level of each musical score.

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

Multiple regression analysis by least squares method

Multiple regression analysis is based on the following concept.

There is an explanatory variable group $ x_n $ and an objective variable $ y . In this case, the explanatory variables are the values of the groove radar. The objective variable is the difficulty level. At this time, $ y' = k_0 + k_1x_1 + k_2x_2 + \cdots + k_nx_n $$ Considering the coefficient groups $ k_n $ and $ y'$, the square error of $ m $ data $ e ^ 2: = \ sum_ {i = 1} ^ m \ left (y'_i --y_i \ right) ) $ K_n $ will be searched so that ^ 2 $ is the smallest. This time, we will find such an optimization problem by the differential evolution method using SciPy.

First, define the function you want to minimize. It's $ 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()

Give this to SciPy's differential_evolution function. As for the search range, I give it a range that seems to be sufficient while trying various things.

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 ])

Looking at this result, it seems that STREAM has the most influence, followed by VOLTAGE, CHAOS, FREEZE, AIR.

Now, let's evaluate using the parameters actually obtained.

First, define a function for prediction. Returns a tensor of the predicted difficulty that gives this function a parameter and a tensor of groove radar.

def pred1(x, sgr_arr):
    est = estimate(x, sgr_arr).clip(1., 19.)
    return np.round(est).astype(int)

A confusion matrix object is created by giving the return value of this function and the actual difficulty to PyCM's ConfusionMatrix. Access the property of this and find the correct answer rate and macro F value.

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
====================

The correct answer rate was 31.4%. This is a much lower result than I expected. Let's heatmap the confusion matrix with Seaborn.

plt.figure(figsize=(10, 10), dpi=200)
sns.heatmap(pd.DataFrame(test_cm1.table), annot=True, square=True, cmap='Blues')
plt.show()

The horizontal axis is the actual difficulty level, and the vertical axis is the predicted difficulty level. Looking at this, it seems that songs with a low difficulty level are high, those with a certain level of difficulty are evaluated low, and those with a higher difficulty level are overestimated. The heat map is not a straight line but bends inward to draw a bow shape.

Maximum likelihood estimation method (order) Logistic regression analysis

Then try logistic regression analysis. Logistic regression analysis is suitable for analysis that takes some probability as the objective variable. Consider the following formula. $ y' = \frac{1}{1 + \exp\left[-\left(k_0 + k_1x_1 + k_2x_2 + \cdots + k_nx_n\right)\right]} $ This $ y'$ takes a value from 0 to 1. When the objective variable $ y $ is determined by 0 or 1 in $ m $ data, the likelihood $ l = \ prod_ {i = 1} ^ m y_iy_i'+ (1 --y_i) (1 --y_i' ) Consider maximizing $. It is difficult to understand if you write it in a mathematical formula, but $ y'$ is the probability of being positive, and if $ y $ is 1 that is positive, then $ y'$ as it is, and if $ y $ is 0, that is, negative. Probability, that is, use the value obtained by subtracting $ y'$ from 1. If you multiply $ y'$, which is the interval of $ (0, 1) $, the value will become too small and difficult to handle on a computer, so log-likelihood $ \ log l = \ sum_ {i = 1} ^ m \ log \ left [y_iy_i'+ (1 --y_i) (1 --y_i') \ right] It is normal to maximize $.

This time we're not dealing with negative or positive, but where in the ordered class they belong. In such a case, assume multiple logistic curves that differ only in the constant term $ k_0 $, and consider them as “probability of difficulty level 2 or higher”, “probability of difficulty level 3 or higher”, ..., respectively. The "probability of difficulty level 2" is the "probability of difficulty level 2 or higher" minus the "probability of difficulty level 3 or higher", so the likelihood can be calculated from this.

First, convert the difficulty level to a one-hot format 2nd floor tensor for calculation.

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.

Then give the function to be minimized. Since it is a minimization, we define the above log-likelihood with a minus.

def upperscore(x, sgr_arr):
    x_const = np.append(np.append(oo, x[:18].copy()), -oo) #Insert infinity at both ends for probability 1 above 1 and probability 0 greater than 20
    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) #Shift and pull,Remove the end to get the probability of each difficulty

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

Perform a search. Please note that it will take considerably longer than before.

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])

This time, STREAM> CHAOS> VOLTAGE> FREEZE> AIR, and you can see that CHAOS has a greater influence.

Let's actually evaluate this as well.

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
====================

This time, the correct answer rate was 54.5%. It's a lot better than before, but it's still a long way off.

plt.figure(figsize=(10, 10), dpi=200)
sns.heatmap(pd.DataFrame(test_cm2.table), annot=True, square=True, cmap='Blues')
plt.show()

This is generally on a straight line, but the low and high difficulty levels are still not good.

Conclusion

The bottom line is, "You can get that kind of value, but it's not practical." I made this attempt with the hope that it could be used as a reference for adding difficulty to the score of StepMania, but in the end it seems that it will be necessary to play and adjust.

Another thing is that in logistic regression, each constant term is naturally constrained by the magnitude relation, but in this code, depending on the random number, the constraint may not be met, and an abnormal value may result in a convergence test [^ bound]. ]. SciPy's optimization function can give constraints with inequalities, but it didn't work for me because I got an error like a constraint of unknown form was passed. Search time is also wasted, so if anyone can solve it, I would like to ask a professor.

[^ bound]: I encountered such a phenomenon once at the stage of adjusting the actual range.