最近，カラオケ採点（LiveDAM の精密採点 DX）にハマっている． カラオケ採点では，「童謡を歌うと高得点がとれる」という話をよく聞くし，実際自分で歌ってみても他の曲より高い点が出ることが多い． そこで，カラオケで高得点がとれる曲というのはどんな曲なのか，逆に難しい曲とは何なのか，ということを調べてみようと思った． 最終的には，ある人がある曲を歌った時の点が予測できたら面白いと思う．

精密採点 DX のデータについて

精密採点 DX の採点方法は，公式の公開情報とスコアラーの試行錯誤によって，ある程度知られている． 精密採点 DX における最終的な得点（総合点）は，基本5項目である音程・安定性・表現力・リズム・ビブラート&ロングトーン（以下VL）の評価，および隠しボーナスの補正から決まるとされる． これらの項目別評価はそれぞれ 100 点満点で行われており，歌うたびにサーバーに記録されていて，あとから DAM の公式サイトで調べることができる．

これらのデータを集約したサービスとして，精密集計 DX (http://clubdam.info/) というサイトがある． 今回は，このサイトから引っ張ってきた 40 人分の歌唱履歴 44,777 件をデータとして使う． データにどんな項目があるかは，精密集計 DX のサンプル http://clubdam.info/user/m/song を参照のこと．

計算方法

モデルを立てる

今回のモデルでは，基本5項目の点数を曲の難易度とユーザーの能力によって説明し，これによって点数を予測する．

ユーザー ，曲 に対して， ユーザー が曲 を歌ったときの項目別得点を 次元ベクトル で表す． 今回は基本5項目だけ使うので で， とする．

に対して，次のようなモデルを考える．

(1)  

と はパラメータ（あるいは潜在変数）で， は 行列， は 次元ベクトルである． は曲ごとの特性， はユーザー固有の能力を表していると考えればよい． ユーザー能力は項目別得点以上の情報を持ちえないので， を想定している．

解釈のために制約を入れる

(1) は次元削減も含んだ一般的なモデルで，実際には に適当な制約を入れる必要がある． 今回は制約条件として とし， が対角行列 であるとする． この制約はかなり厳しいものではあるが，次のような利点がある．

• と が同じ空間になり， の各要素が基本項目に対応するようになるので，ユーザー能力の解釈が簡単になる． たとえば， の1つ目の要素はユーザー が持つ（曲に依らない）音程センスだと解釈できる．
• が対角行列なので，各項目が独立した線形モデルとして計算できる． または のどちらかを固定してしまえば線形単回帰になり，簡単に計算できる．

以下，ユーザー能力は の5次元ベクトルとして扱う．

切片項をどうするか

(1) のモデルでは切片項を入れていない．切片項を入れようとすると， と の対称性を考えて次のモデルを使うことになる．

(2)  

切片項付きのモデルを使わない理由は， の解釈が難しくなるというのが大きい． 切片項を入れない場合は と の原点がそろっているので，「ユーザー本来の能力はこれくらいだが選曲によって○倍になった」という解釈ができるが， 切片項を入れた場合はこれができなくなる． もちろん，切片項を入れた場合はユーザー能力が「選曲に依存する部分」と「しない部分」に分かれて，それはそれで面白いが……

切片項を入れない理由としてはもう一つ，自由度が著しく上がってしまうというのもある． 今回のデータでは「1曲に対して歌っている人が1人しかいない」ようなケースがあり，曲ごとのパラメータを増やしてしまうと推定が不安定になる． なので，曲ごとの情報はすべて に押し込んでしまって，切片項を付けないほうがよさそうだと判断した．

同じユーザーが同じ曲を歌う

一般に，同じユーザーが同じ曲を何回も歌うのはごくありふれたことである． 当然歌うたびに違う得点が出るわけだが，当然この差はホワイトノイズではなく，バイアスが含まれている． たとえば回数を重ねれば，だんだん得点は上がっていくことが予想できる．

(1) のモデルでは，ユーザー・曲が同じ歌唱履歴がデータに含まれていても，推定自体は問題なく実行できる． しかし，バイアスは全く考慮していないので，このようなデータは推定結果を悪くする原因となりうる（実際になった）． そこで今回は，同じユーザー・曲を含む歌唱履歴は，総合点の一番高いものだけを使うことにした．

もちろん他にベターな解決方法はいろいろあると思うので，ここは要検討．

最適化

と を求めるために，次の最適化問題を解く．

(3)  

ここで， 番目の歌唱履歴に対して， は項目別点数， は歌ったユーザー， は歌った曲を表す． 歌唱履歴の数 はユーザー数 ，曲数 とはまったく関係ないことに注意．

前述のとおり， と の片方を固定すれば線形単回帰で計算できるので，交互最適化によってこの問題を解く． 各パラメータはすべて独立に推定することができ，たとえばユーザー の音程能力 は次のように推定できる．

(4)  

ここで はクロネッカーのデルタで，2つの引数が等しいとき1となる関数である． つまり，この場合はユーザー が含まれる歌唱履歴だけ取ってきて， が説明変数であると思って単回帰を実行すればよい．

を推定するときには を説明変数と思って，同様に更新を行う．

基本5項目の評価予測

というわけでいろいろ説明したが，ここからは実際にデータを使って分析した結果を書いていく． まずは予測の結果から．

予測誤差を計算するため，データの 5% をテストデータ，残りを訓練データとして，MSE で精度を比較した． 比較対象のベースラインとして，単にユーザーごとの曲間平均を取ったものを使う． たとえば，ユーザー が曲 を歌った時の音程の予測値は で与える．

それぞれの MSE は次のようになった．

	Proposed	Baseline
MSE	66.593	74.928

多少はベースラインに勝っているとはいえ，まあなかなか難しそう． ちなみに，音程についての真値対予測値のプロットがこちら．

やはり誤差分散がかなり大きく予測は難しそうというのが正直なところ． ただ，あまりバイアスは出ていないようなので，大域的には線形モデルでも十分対応できそうではある．

上の図についてあえていいところを挙げておくと，ベースラインは平均を使っているので 100 点に近づくほど下方修正がかかっているが， 提案手法ではきちんと 100 点の近くまで予測できている． もちろん提案手法は 100 点を超える予測値を出す可能性もあるので，何とも言えないけど．

訓練データの選別

実は上の分析では，訓練データのうち4回以上歌われている曲だけを使って分析を行っている． というのも，すべてのデータを使うと，かえって MSE が悪くなってしまったからである．

一応，しきい値を「全部使う」から「6 回以上」まで動かして MSE を比較してみると，次のようになる．

	全部	2-	3-	4-	5-	6-
MSE	85.943	71.748	70.230	66.593	72.045	73.160
データ数	9040	5206	3694	2704	2036	1511

データを全部使う場合，「1回しか歌われていない曲」が 4,000 曲近くあるので， 推定精度は極めて悪い．しきい値を増やしていくと精度は上がっていくが， データ数は減っていくのでどこかで頭打ちになる． 今回の場合は， 4 回以上歌われた曲を使った場合がベストであるという結果になった．

曲の難易度

前述のとおり，推定された の値の大小を見ることで，曲の難易度を測ることができる． が大きければ「個人の能力以上の点が出る曲」つまり易曲，小さければ「個人の能力の割に点が伸びない曲」つまり難曲ということになる．

音程について， の大きい曲を並べてみるとこんな感じ．

	曲名	アーティスト
1.126853	キラキラぼし	童謡
1.103351	君が代	国歌
1.084746	故郷	唱歌・抒情歌
1.080799	お正月	童謡
1.079042	仰げば尊し(生音)	唱歌・抒情歌
1.072712	Won(*3*)Chu KissMe!	SAKURA*TRICK
1.072103	ウミ	唱歌・抒情歌
1.069104	君は僕に似ている	See-Saw
1.066168	仰げば尊し	唱歌・抒情歌
1.061675	夏の約束	堀江由衣

見事に上位を童謡が占める結果に． トップのキラキラぼしは が 1.12 もあるので，たとえば他の曲ではだいたい音程 80 しか取れない人でも，この曲を歌うと音程 90 がとれることになる．スゴイ！

逆に の小さい順に並べてみると……

	曲名	アーティスト
0.8492424	エアーマンが倒せない	Team.ねこかん[猫]
0.8538813	STAND PROUD	橋本仁
0.8550817	M	浜崎あゆみ
0.8662832	キセキ	GReeeeN
0.8713905	ジョジョ～その血の運命～(生音)	富永TOMMY弘明
0.8926162	ultra soul(生音)	B’z
0.9021115	決意の朝に(生音)	Aqua Timez
0.9025316	アンハッピーリフレイン	wowaka feat.初音ミク
0.9046287	炉心融解	iroha feat.鏡音リン
0.9076202	君じゃなきゃダメみたい	オーイシマサヨシ

シャウト系(?)が入ってきてなるほどという感じ．エアーマンはテンポ速い・音域広い・跳躍多いで納得のトップ． 浜崎は何が難しいのかわからない．今度歌ってみよう．

改善点

結果を並べていて気付いたこと． 今回は曲の難易度とユーザー能力の積で得点を予測したが，単に和でもいいのでは？という気がした． つまり，こんなモデルでもいいかもしれない．

(5)  

さすがにこれだと簡単すぎるかもしれないが，現状線形モデルでそこそこ説明できているので，こういうのもアリかもしれない． (5) はランダム効果モデルの式と似ているので，そちらのノウハウが使えるかも．

別の正攻法的アプローチとしては，今回入れていた「 が対角」という条件を緩めることもできる． このタイプの方向性としては，たとえば PCA で使うような制約を入れたりして計算するなどがある． が対角でなくなると各項目間の相関も取れるようになるので，また違った結果が出てくるだろう． 少なくとも， の次元は よりは落とせそうである．

逆に，高次の項を導入しても，それほど結果は変わらないのではないかと思う． より複雑なモデルにするとしたら，たとえば 95 ～ 100 点のあたりを部分的にノンパラにして，端のひずみに対応するなどがある． とくにスコアラーにとっては 100 点近くの精度が非常に重要なので，この方向性もいいかもしれない．

予測はこのデータからは厳しそうではあるが，今回は無視していた歌唱日時などのデータを使えば十分改善できそうではある．

いろいろ試してみて，また面白い結果が出れば報告します．

ちなみに最初に紹介した集計サイトでも，新機能として似たようなことをやろうとしているようだ．期待．

おまけ

こんな論文を見つけた．ディープなにおいがする．

今回使用したスクリプト

## grid / ggplot —- library(ggplot2) library(grid) make.grid <- function(row, col) { grid.newpage() l <- grid.layout(row, col) v <- viewport(layout=l) pushViewport(v) } print.at <- function(o, i, j) { print(o, vp=viewport(layout.pos.row=i, layout.pos.col=j)) } end.grid <- function() { popViewport() } ## データ読み込み ---- dat <- read.csv("user.csv", fileEncoding="utf-8") elemname <- c("chartInterval", "chartStability", "chartExpressiveness", "chartVibratoLongtone", "chartRhythm") p <- length(elemname) lv.user <- levels(dat$user) lv.reqno <- levels(dat$requestNo) n.user <- length(lv.user) n.reqno <- length(lv.reqno) # (曲, ユーザー) ペアで総合点が最も高いものだけ使う filt.idx <- tapply(1 : nrow(dat), factor(paste0(dat$requestNo, dat$user)), function(idx) { sub.dat <- dat[idx, ] idx[which.max(dat[idx, "totalPoint"])] }) dat <- dat[filt.idx, ] # pitch を数値に変換 replace <- function(reptb, x) { for (i in 1 : nrow(reptb)) x <- gsub(reptb[i, 1], reptb[i, 2], x); x } dat$highPitch <- gsub("♭", "b", dat$highPitch) dat$lowPitch <- gsub("♭", "b", dat$lowPitch) reptb.pitch <- data.frame( c("lowAb", "lowA", "lowBb", "lowB", "lowC", "lowDb", "lowD", "lowEb", "lowE", "～lowF", "lowGb", "lowG", "m1Ab", "m1A", "m1Bb", "m1B", "m1C", "m1Db", "m1D", "m1Eb", "m1E", "m1F", "m1Gb", "m1G", "m2Ab", "m2A", "m2Bb", "m2B", "m2C", "m2Db", "m2D", "m2Eb", "m2E", "m2F", "m2Gb", "m2G", "hihiAb", "hihiA", "hihiBb", "hihiB～", "hihiC", "hihiDb", "hihiD", "hihiEb", "hihiE", "hihiF", "hihiGb", "hihiG", "hiAb", "hiA", "hiBb", "hiB", "hiC", "hiDb", "hiD", "hiEb", "hiE", "hiF", "hiGb", "hiG", "Ab", "A", "Bb", "B", "C", "Db", "D", "Eb", "E", "F", "Gb", "G" ), c(32 : 43, 44 : 55, 56 : 67, 80 : 91, 68 : 79, 68 : 79) ) dat$highPitch <- as.numeric(replace(reptb.pitch, dat$highPitch)) dat$lowPitch <- as.numeric(replace(reptb.pitch, dat$lowPitch)) # 曲テーブルを作成 songs <- data.frame( reqno=lv.reqno, artist=factor(tapply(as.character(dat$artist), dat$requestNo, function(x) x[1])), contents=factor(tapply(as.character(dat$contents), dat$requestNo, function(x) x[1])), highPitch=tapply(dat$highPitch, dat$requestNo, function(x) x[1]), lowPitch=tapply(dat$lowPitch, dat$requestNo, function(x) x[1]) ) songs$diffPitch <- songs$highPitch - songs$lowPitch + 1 ## 訓練・テストデータを作る ---- # 全体の 5% をテストに回す set.seed(0) index.test <- sample(1 : nrow(dat), nrow(dat) * 0.05) dat.test <- dat[index.test, ] dat.train <- dat[-index.test, ] # 訓練データにない曲を使うテストデータをはじく（要検討） dat.test <- dat.test[table(dat.train$requestNo)[dat.test$requestNo] > 0, ]

# Temporary variables y <- as.matrix(dat.train[, elemname]) user.test <- as.numeric(dat.test$user) reqno.test <- as.numeric(dat.test$requestNo) y.test <- as.matrix(dat.test[, elemname]) # データ数のチェック c(train=nrow(dat.train), test=nrow(dat.test)) ## ベースライン: ユーザー内平均を予測値とする ---- mse <- function(true.y, pred.y) { sum((pred.y - true.y)^2) / length(true.y) } pred.y.by.mean <- apply(y, 2, function(y) tapply(y, dat.train$user, mean))[user.test, ] mse(y.test, pred.y.by.mean) ## diaglm 実装 ---- library(parallel) diaglm <- function(dat, threshold=4, weighted=FALSE, verbose=TRUE, cl=NULL) { if (!is.null(cl)) sapply <- function(...) parSapply(cl, ...) # 訓練には threshold 回以上出現する曲だけ使う dat <- dat[table(dat$requestNo)[dat$requestNo] >= threshold, ] reqno <- as.numeric(dat$requestNo) user <- as.numeric(dat$user) y <- as.matrix(dat[, elemname]) # 初期値の設定 a <- matrix(1, n.reqno, p) x <- apply(y, 2, function(y) tapply(y, dat$user, mean)) if (weighted) { w <- apply(y, 2, function(y) table(y)[as.character(y)]) } else { w <- matrix(1, nrow(dat), 5) } # 交互最適化 iter <- 0 err <- Inf while (err > 1e-07) { iter <- iter + 1 a0 <- a x0 <- x j <- which(table(dat$requestNo) > 0) a[j, ] <- t(sapply(j, function(j, y, xx, w, reqno, user) { sub.y <- y[reqno == j, , drop=FALSE] sub.x <- xx[user[reqno == j], , drop=FALSE] sub.w <- w[reqno == j, , drop=FALSE] sub.wx <- sqrt(sub.w) * sub.x diag((t(sub.wx) %*% sub.y) / (t(sub.wx) %*% sub.x)) }, y=y, xx=x, w=w, reqno=reqno, user=user)) x <- t(sapply(1 : n.user, function(i, y, a, w, reqno, user) { sub.y <- y[user == i, , drop=FALSE] sub.a <- a[reqno[user == i], , drop=FALSE] sub.w <- w[user == i, , drop=FALSE] sub.wa <- sqrt(sub.w) * sub.a diag((t(sub.wa) %*% sub.y) / (t(sub.wa) %*% sub.a)) }, y=y, a=a, w=w, reqno=reqno, user=user)) err <- sum((a - a0)^2) / (n.reqno * p) + sum(((x - x0) * 0.01)^2) / (n.user * p) if (verbose) { cat(paste0("#", iter, ": ", round(err, digits=8), "\n")) } } resid <- a[reqno, ] * x[user, ] - y r2 <- sapply(1 : p, function(k) { (var(y[, k]) - sum(resid[, k]^2) / nrow(dat)) / var(y[, k]) }) colnames(a) <- elemname colnames(x) <- elemname list(a=a, x=x, resid=resid, r2=r2) } ## 推定 ---- cl <- makeCluster(4) # diaglm diaglm.std <- diaglm(dat.train, threshold=4, weighted=FALSE, cl=cl) pred.y.by.diaglm.std <- diaglm.std$a[reqno.test, ] * diaglm.std$x[user.test, ] c(mse=mse(y.test, pred.y.by.diaglm.std), r2=diaglm.std$r2) # weighted diaglm diaglm.w <- diaglm(dat.train, threshold=0, weighted=TRUE, cl=cl) pred.y.by.diaglm.w <- diaglm.w$a[reqno.test, ] * diaglm.w$x[user.test, ] c(mse=mse(y.test, pred.y.by.diaglm.w), r2=diaglm.w$r2) # threshold summary.diaglm.t <- t(sapply(1 : 6, function(threshold) { diaglm.t <- diaglm(dat.train, threshold=threshold) pred <- diaglm.t$a[reqno.test, ] * diaglm.t$x[user.test, ] c(mse=mse(y.test, pred), r2=diaglm.t$r2) })) summary.diaglm.t stopCluster(cl) ## Interval のみで mean, diaglm を比較 int.df <- data.frame( true = y.test[, "chartInterval"], pred.mean = pred.y.by.mean[, "chartInterval"], pred.diaglm = pred.y.by.diaglm.std[, "chartInterval"] ) g1 <- ggplot(data=int.df, aes(x=true, y=pred.diaglm)) + geom_point() + geom_abline(slope=1) + xlim(50, 100) + ylim(50, 100) + xlab("True int") + ylab("Predicted int (proposed)") g2 <- ggplot(data=int.df, aes(x=true, y=pred.mean)) + geom_point() + geom_abline(slope=1) + xlim(50, 100) + ylim(50, 100) + xlab("True int") + ylab("Predicted int (baseline)") svg("interval.svg", width=8, height=4) make.grid(1, 2) print.at(g1, 1, 1) print.at(g2, 1, 2) end.grid() dev.off() ## 残差の大きいサンプルを見る ---- int.df$resid.diaglm <- int.df$pred.diaglm - int.df$true int.df$resid.mean <- int.df$pred.mean - int.df$true int.df <- int.df[order(abs(int.df$resid.diaglm), decreasing=TRUE), ] int.df[1:20, c("true", "resid.diaglm", "resid.mean")] ## 音程の取りづらい曲は？ ---- cl <- makeCluster(4) diaglm.full <- diaglm(dat, threshold=4, weighted=FALSE, cl=cl) stopCluster(cl) df <- data.frame( a=diaglm.full$a[, 1], songs[, c("diffPitch", "contents", "artist")] ) df[order(df$a), ][1 : 10, ] df[order(df$a, decreasing=TRUE), ][1 : 10, ] [/code]