前回「精密採点 DX のデータから曲の難易度を推定する」の記事 では，ユーザーと曲から基本5項目の項目別点数を推定した． 今回は，この項目別点数が実際に総合点にどのように影響しているのかを，ノンパラ回帰を使って推定する．

データについて

データは前回と同じく，精密採点 DX における 40 人分 44,777 件の歌唱履歴である．このうち，総合点 100 のサンプルは何らかの打ち切りを受けている可能性が高いので，この条件にマッチする 1,875 件は除外した．ちなみにこれらのサンプルは，入れても入れなくても結果的にほとんど影響はなかった．

ノンパラ回帰による分析

今回は，音程・安定性・表現力・リズム・VL の基本5項目の点数を説明変数とし，サンプルベースのノンパラメトリック回帰によって分析をした．通常の線形回帰でなくノンパラを用いたのは，次のような理由による．

• 真のモデルは精密採点 DX の開発者によって人工的に作られたモデルである．したがって，各項目の総合点への影響は全体としてはおそらく線形ではないが，ごく簡単な関数の組み合わせによって表せると推測される．
• 項目別点数は 0 ～ 100 の値を取る離散変数なので，5項目分をすべて合わせてもパラメータの個数は高々 500 しかない．サンプル数は 4 万以上あるので，ノンパラでも十分に推定できる．

パラメータの推定は，歌唱履歴 に対する を総合点， を項目別点数として，最小二乗法によって行う．

(1)  

ここで は基本5項目のインデックス集合， はクロネッカーのデルタ． がパラメータである．データから原点が 0 であることが分かっているので，切片項は含めないことにした．

この推定は，各項目点をダミー変数で置きなおせば，通常の線形回帰のソルバーで解くことができる．

残差についての考察

モデルが立てられたのでさっそく推定した係数を見てみたいところではあるが，その前に残差について興味深い結果がある．次のグラフは，残差のヒストグラムである．

ヒストグラムを見ると，平均がおよそ 0 の大きな山のほかに，平均 -7 ほどにも山があるのがわかる．大きな山は通常の残差の分布である．一方で，小さい山の方も左右に対称であり，しかもバラつきの幅は大きな山とほぼ同じである．このことは，「総合点を一定値下げる何らかの要因（ペナルティ）がある」ことを示唆している．

これについて，データに含まれている自分の歌唱履歴のうち残差が極端に小さいものを調べてみた．その結果，残差が極端にずれているのはすべて「途中で演奏停止した」場合で，自分がいままで演奏停止した曲が過不足なく含まれていた．したがって，このペナルティは，演奏停止によるペナルティだと思って間違いないだろう．

残差については，もうひとつ面白いプロットがある．次の図は，残差をユーザーごとに並べて3色で色分けしたものである．

これを見ると，「通常の残差」の中でも，ユーザーごとに多少の差が表れている．たとえば，一番右側の人は残差の平均が -0.5 ほどであるのに対し，右から2番目の人は 1.5 ほどもあり，分布も明らかに異なっている．

これは，スコアラーの間でよく知られている「裏加点」による影響だと思われる．裏加点は声の倍音の強さなどによると言われているので，ユーザーごとに裏加点の乗りやすさが異なるのは十分にあり得ることである．もちろん同じユーザーの中でもバラつきがあり，部屋やマイクによっても変わってくると思われるが，残念ながらこれらの情報はデータに含まれていないので，推定に加えることはできない．幸い，ユーザー情報を付加すると，全体の残差はほぼ正規分布に収まったので，推定結果への影響は小さそうではある．

なお，残差で -5 ～ -10 あたりに出ている点が，先ほどの演奏停止ペナルティ付きの点である．ユーザーによってここの密度がかなり違っており，演奏停止をよく押す人とほとんど押さない人の差が出ているのがわかる．

ノンパラ回帰による分析（改）

残差についての考察を踏まえて，次のように推定方法を変更する．

• ユーザーをダミー変数として説明変数に加える．
• ペナルティを受けていると思われるサンプルを除外する．

ユーザーを説明変数に加えるので，目的関数は次のようになる．

(2)  

ユーザーごとの裏加点の平均が 0 になってほしいので，推定の際には の制約を入れている．

演奏停止をしたかどうかのフラグはデータに含まれていないので，この判定には外れ値検出の方法を用いる．今回は，ユーザーごとの残差の 40% 点・中央値・60% 点を使って正規近似を計算し，99.99%点より外側のサンプルを除外した．全サンプル数が 4 万程度なので，もし全体が正規分布に従うなら，せいぜい数個のサンプルだけが除外される計算になる．

次の図では，外れ値と判定されたサンプルを赤で示している．

推定された係数

それでは，いよいよ推定結果の紹介に移る．次の図は，推定された を項目別にプロットしたもので，各項目の を足し合わせると（裏加点を考慮しない）総合点の予測値となる．たとえばすべての項目で 80 を取ったとすると， の値はそれぞれ 72.4, 1.9, 4.2, 3.5, 2.1 なので，総合点はおよそ 84.1 だと計算できる．

音程以外の4項目については，80 と 90 にはっきりとした境目があり，分割されたそれぞれの区間では線形となっているのがわかる．これは想定していた通りの結果で，80 と 90 を境に関数形を変えているのはほぼ間違いないとと思われる．

この図から，総合点の 7 ～ 8 割は音程によって決まるが，音程以外の項目で 80 以上が取れるようになれれば，総合点が飛躍的に伸びることがわかる．実際，音程以外の4項目で 100 を取ったときの の合計がおよそ 22.3 なので，音程が 72 あれば 90 点，80 もあれば 95 点まで到達可能となる．総合 90 点以下で音程に自信のない人は，ビブラートや表現力を練習すると得点アップにつながるだろう．

次に，各項目 80 以上の部分を拡大してプロットしたのが下の図である．この図では各項目での最高値が 0 となるようにスケールを合わせている．

先にコメントしておくと，今回の分析の結果では，すべての項目で最高点を取ったとしても（音程 95 ・それ以外 100）， の合計は 99.886 となり，100 点にわずかに届かない．もっともこの差は誤差の範囲内であるし，実際には裏加点があるので 100 点を超えることは十分に可能である．一応図中には 100 点のボーダーとなるラインを書き込んでおいた（0 の上の黒い横線）．

音程は 95 でピークを取るまで，他より急な傾きを保っている．たとえば音程 80 での最高到達点は 95 点だが，音程 90 まで上げれば最高到達点は 99 点まで伸びる．一方で，他の項目は 90 を超えると伸び率が鈍っている．したがって，95 点前後の人が 98 ～ 99 点を目指すには，90 以下の苦手項目をひとまず 90 まで乗せてから，音程を 95 まで上げていくのが効率的だといえるだろう．

一方，99 点より上を目指す場合は，苦手な項目が一つでもあると厳しくなる．リズム以外でひとつでも 90 を取ってしまうと，およそ -1 点の減点となるので，総合で 99 点にすら乗らないことになる．100 点を取るには，まず音程を 95% にし，比較的稼ぎやすいリズムと表現をしっかり最高値近く (98 以上) まで上げたうえで，安定性と VL を噛み合わせるしかない．

以上の結果は，これまでのスコアラーの知見とほぼ一致しているので，推定はおおよそ正しくできたとみてよいだろう．

裏加点について

最後に，ユーザーごとの裏加点の乗り方を見てみたい．次の図は，ユーザーごとの裏加点パラメータに残差を加えたボックスプロットである．

こうしてみると，人によって裏加点の乗り方はかなり違っているのがわかる．高い人は平均で 1.7 もの裏加点が乗っているが，1.7 も加点されるということは音程が 88 でも総合 100 点に到達可能という計算になる．逆に，一番右の人などはヒゲの上端すら 0 に届くか届かないかなので，総合 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() } plot.to.file <- FALSE ## データの読み込み ---- dat <- read.csv("user.csv", fileEncoding="utf-8") dat$date <- as.POSIXct(dat$date, tz="JST") dat <- dat[order(dat$user, dat$date), ] enam <- c("chartInterval", "chartStability", "chartExpressiveness", "chartRhythm", "chartVibratoLongtone") ## とりあえずlm ---- lm.std <- lm(totalPoint ~ chartInterval + chartStability + chartExpressiveness + chartVibratoLongtone + chartRhythm, data=dat) summary(lm.std) ## ノンパラlm ---- mask.nonp <- rep(TRUE, nrow(dat)) # 総合点100点を除外 (打ち切り対策) mask.nonp <- mask.nonp & dat$totalPoint < 100 ## 各項目30以上だけ使う (singular 対策) #mask.nonp <- mask.nonp & apply(dat[, enam], 1, function(row) all(row[enam] >= 30))

# Singular 対策 mask.nonp <- mask.nonp & dat$chartRhythm != 0 mask.nonp <- mask.nonp & dat$chartRhythm != 2 mask.nonp <- mask.nonp & dat$chartExpressiveness != 13 # 推定 do.lm <- function(mask.nonp, use.user=FALSE) { enam <- if (use.user) c(enam, "user") else enam dat.lm <- dat[mask.nonp, c("totalPoint", enam)] for (en in enam) dat.lm[[en]] <- factor(dat.lm[[en]]) design.lm <- data.frame(model.matrix(~ . - 1, dat.lm)) if ("chartInterval0" %in% colnames(design.lm)) { design.lm <- design.lm[, -which(colnames(design.lm) == "chartInterval0")] } if (use.user) { ucol <- paste0("user", levels(dat$user)[-1]) design.lm[, ucol] <- design.lm[, ucol] - as.numeric(dat$user[mask.nonp] == levels(dat$user)[1]) } lm(totalPoint ~ . - 1, design.lm) } lm.nonp <- do.lm(mask.nonp) summary(lm.nonp) # 係数のプロット coef.nonp <- coef(lm.nonp) print.coef <- function(coef.nonp, lower=40) { value <- c(-1, lower : 100) df <- data.frame( elements = factor(rep(enam, each=length(value)), levels=enam), value = rep(value, 5), beta = as.vector(sapply(enam, function(en) { coef.nonp[paste0(en, value)] })) ) f <- function(en) { ggplot(df[df$elements == en | df$value == -1, ], aes(x=value, y=beta)) + geom_point(aes(color=elements)) + scale_x_continuous(en, limits=c(lower, 100), breaks=seq(lower, 100, 10)) + theme(legend.position="none") } make.grid(3, 2) print.at(f("chartInterval"), 1, 1) print.at(f("chartStability"), 1, 2) print.at(f("chartExpressiveness"), 2, 1) print.at(f("chartRhythm"), 2, 2) print.at(f("chartVibratoLongtone"), 3, 1) end.grid() } #print.coef(coef.nonp) # 残差の分布を見る resid.nonp <- resid(lm.nonp) print.resid.point <- function(resid.nonp, mask.nonp) { df <- data.frame( index = 1 : length(resid.nonp), resid = resid.nonp, user = factor(as.numeric(dat$user[mask.nonp]) %% 3, levels=0:2) ) make.grid(1, 1) print.at( ggplot(df, aes(x=index, y=resid)) + geom_point(aes(color=user), size=1) + theme(legend.position="none"), 1, 1) end.grid() } print.resid.hist <- function(resid.nonp) { df <- data.frame(resid = resid.nonp) make.grid(1, 1) print.at( ggplot(df, aes(x=resid)) + geom_histogram(aes(fill=0)) + theme(legend.position="none"), 1, 1) end.grid() } if (plot.to.file) svg("total-resid-index.svg", width=6, height=4) print.resid.point(resid.nonp, mask.nonp) if (plot.to.file) dev.off() if (plot.to.file) svg("total-hist-resid.svg", width=6, height=4) print.resid.hist(resid.nonp) if (plot.to.file) dev.off() ## 外れ値を除去 ---- # ユーザーごとで 99.99% 点より外側のサンプルを除外 quant.resid <- tapply(resid.nonp, dat$user[mask.nonp], function(r) quantile(r, probs=c(0.4, 0.5, 0.6))) mask.ol <- mapply( function(resid, user) { quant <- quant.resid[[user]] resid >= quant[“50%”] – (quant[“60%”] – quant[“40%”]) * (qnorm(0.9999) / (2*qnorm(0.6))) }, resid.nonp, dat$user[mask.nonp] ) mask.nonp2 <- mask.nonp mask.nonp2[mask.nonp] <- mask.ol # 判定された外れ値を見る olresid.df <- data.frame( index = 1 : length(resid.nonp), resid = resid.nonp, col = ifelse(mask.nonp2[mask.nonp], "darkgray", "red") ) if (plot.to.file) svg("total-resid-index-ol.svg", width=6, height=4) make.grid(1, 1) print.at( ggplot(olresid.df, aes(x=index, y=resid)) + geom_point(color=olresid.df$col, size=1) + theme(legend.position="none"), 1, 1) end.grid() if (plot.to.file) dev.off() # 推定 lm.nonp2 <- do.lm(mask.nonp2, use.user=TRUE) summary(lm.nonp2) # 係数のプロット coef.nonp2 <- coef(lm.nonp2) if (plot.to.file) svg("total-coef-elem.svg", width=8, height=12) print.coef(coef.nonp2, lower=40) if (plot.to.file) dev.off() # 各項目 80 以上での拡大図 df <- data.frame( elements = factor(rep(enam, each=21), levels=enam), value = rep(80 : 100, 5), beta = as.vector(sapply(enam, function(en) { ecoef <- coef.nonp2[paste0(en, 80 : 100)] ecoef - max(ecoef, na.rm=TRUE) })) ) border <- 100 - sum(coef.nonp2[paste0(enam, 100)[-1]]) - coef.nonp2["chartInterval95"] border if (plot.to.file) svg("total-coef-elem-o80.svg", width=6, height=4) make.grid(1, 1) print.at( ggplot(df, aes(x=value, y=beta, group=elements)) + geom_hline(yintercept=border, color=1) + geom_point(aes(color=elements)) + geom_line(aes(color=elements)) + xlim(c(80, 100)), 1, 1) end.grid() if (plot.to.file) dev.off() # 残差の分布を見る print.resid.hist(resid(lm.nonp2)) # 裏加点の box-plot ura.per.user <- coef.nonp2[paste0("user", levels(dat$user)[-1])] ura.per.user <- c(-sum(ura.per.user), ura.per.user) df <- data.frame( ura = resid(lm.nonp2) + ura.per.user[as.numeric(dat$user[mask.nonp2])], user = dat$user[mask.nonp2] ) df <- df[order(ura.per.user[as.numeric(dat$user[mask.nonp2])], decreasing=TRUE), ] df$user <- factor(as.numeric(df$user), levels=unique(as.numeric(df$user))) df$ucol <- factor(as.numeric(df$user) %% 3, levels=0:2) if (plot.to.file) svg("total-resid-user.svg", width=6, height=4) make.grid(1, 1) print.at( ggplot(df, aes(user, ura)) + geom_boxplot(aes(color=ucol), outlier.size=1) + theme(legend.position="none", axis.ticks=element_blank(), axis.text.x=element_blank()), 1, 1) end.grid() if (plot.to.file) dev.off() [/code]