最強のM-1漫才師は誰だ
はじめに
この記事はStan Advent Calendar2017，12月11日のエントリー記事です。
Stanというかベイズ統計を勉強するには，犬4匹本もいいんですが，コワい本もまたいいですね。
コワい本の中に，７人の科学者というお題がありまして。７人がある実験でスコアをとったのだけど，どうも最初の二人は未熟者で変なスコアを出している。残りの５人の熟練者から考えるとスコアは10.00点ぐらいになるはずなんだけど・・・という話があります。これをモデリングするにあたって，真の値は一つなんだけど，評定者ごとに「測定誤差」を持っている，というモデルを立ててます。すなわち，科学者$$k$$によって得られた計測値を$$X_{k}$$とすると，$$X_{k} \sim N(\mu,\sigma_k)$$というわけです。
これを応用して，「評定者のくせを考慮したスコア」というのを算出することを考えました。 データはそう，M-1グランプリです。
データの出典はこちらです。 そして今年の分はビデオを見ながら手入力しました。
結論；ブラマヨが一番すごいんじゃないか
この結論に至る手続きを，以下解説いたします。
手続き
データを取り込みます。審査員は毎回変わりますから，欠損値が多いデータセットです。 元のファイルも置いておきますね。UTF-8のcsvファイルをこちらからダウンロードしてください。
## 年代 演者 上沼恵美子 松本人志 ## Min. : 1 笑い飯 : 9 Min. :81.00 Min. :70.00 ## 1st Qu.: 4 麒麟 : 5 1st Qu.:89.00 1st Qu.:85.00 ## Median : 7 ハライチ : 4 Median :90.50 Median :89.00 ## Mean : 7 フットボールアワー: 4 Mean :91.17 Mean :88.12 ## 3rd Qu.:10 千鳥 : 4 3rd Qu.:95.00 3rd Qu.:93.00 ## Max. :13 POISON GIRL BAND : 3 Max. :98.00 Max. :97.00 ## (Other) :90 NA's :73 NA's :46 ## 博多大吉 春風亭小朝 中川礼二 渡辺正行 ## Min. :50.00 Min. :65.00 Min. :87.00 Min. :75.00 ## 1st Qu.:70.00 1st Qu.:80.00 1st Qu.:89.00 1st Qu.:87.00 ## Median :84.00 Median :88.00 Median :90.00 Median :89.00 ## Mean :79.13 Mean :84.28 Mean :90.71 Mean :88.16 ## 3rd Qu.:91.00 3rd Qu.:90.00 3rd Qu.:92.25 3rd Qu.:90.50 ## Max. :97.00 Max. :95.00 Max. :95.00 Max. :95.00 ## NA's :72 NA's :90 NA's :91 NA's :64 ## オール巨人 増田英彦 岩尾望 吉田敬 ## Min. :79.00 Min. :85 Min. :85.00 Min. :83.00 ## 1st Qu.:86.00 1st Qu.:87 1st Qu.:89.00 1st Qu.:85.00 ## Median :88.50 Median :89 Median :90.00 Median :89.00 ## Mean :88.46 Mean :89 Mean :90.67 Mean :87.56 ## 3rd Qu.:91.00 3rd Qu.:91 3rd Qu.:92.00 3rd Qu.:90.00 ## Max. :96.00 Max. :93 Max. :96.00 Max. :93.00 ## NA's :73 NA's :110 NA's :110 NA's :110 ## 徳井義実 富澤たけし 石田明 佐藤哲夫 ## Min. :88.00 Min. :89.00 Min. :83 Min. :88.00 ## 1st Qu.:89.00 1st Qu.:91.00 1st Qu.:87 1st Qu.:89.00 ## Median :89.00 Median :92.00 Median :88 Median :90.00 ## Mean :90.56 Mean :91.89 Mean :89 Mean :90.56 ## 3rd Qu.:91.00 3rd Qu.:93.00 3rd Qu.:92 3rd Qu.:92.00 ## Max. :96.00 Max. :94.00 Max. :94 Max. :93.00 ## NA's :110 NA's :110 NA's :110 NA's :110 ## 哲夫 島田紳助 南原清隆 大竹まこと ## Min. :87.00 Min. : 50.00 Min. :76.00 Min. :75.00 ## 1st Qu.:88.00 1st Qu.: 80.25 1st Qu.:84.75 1st Qu.:82.00 ## Median :90.00 Median : 86.00 Median :88.00 Median :85.00 ## Mean :89.78 Mean : 84.77 Mean :87.78 Mean :85.18 ## 3rd Qu.:91.00 3rd Qu.: 90.00 3rd Qu.:90.25 3rd Qu.:89.00 ## Max. :93.00 Max. :100.00 Max. :98.00 Max. :97.00 ## NA's :110 NA's :37 NA's :83 NA's :47 ## 宮迫博之 中田カウス 東国原英夫 ラサール石井 ## Min. :88 Min. :79.00 Min. :85.00 Min. :68.00 ## 1st Qu.:90 1st Qu.:86.00 1st Qu.:86.00 1st Qu.:82.00 ## Median :91 Median :90.00 Median :88.00 Median :86.00 ## Mean :92 Mean :89.98 Mean :87.67 Mean :86.84 ## 3rd Qu.:93 3rd Qu.:95.00 3rd Qu.:89.00 3rd Qu.:92.00 ## Max. :98 Max. :98.00 Max. :92.00 Max. :96.00 ## NA's :110 NA's :38 NA's :110 NA's :64 ## 島田洋七 西川きよし 立川談志 鴻上尚史 ## Min. :75.00 Min. :75.00 Min. :50 Min. :73.00 ## 1st Qu.:84.00 1st Qu.:80.00 1st Qu.:70 1st Qu.:82.25 ## Median :90.00 Median :86.00 Median :70 Median :83.50 ## Mean :88.09 Mean :85.68 Mean :70 Mean :81.80 ## 3rd Qu.:92.00 3rd Qu.:90.00 3rd Qu.:70 3rd Qu.:84.00 ## Max. :98.00 Max. :97.00 Max. :80 Max. :85.00 ## NA's :74 NA's :100 NA's :110 NA's :109 ## 青島幸男 ## Min. :75.00 ## 1st Qu.:76.25 ## Median :80.00 ## Mean :81.50 ## 3rd Qu.:85.00 ## Max. :90.00 ## NA's :109
データは次のようにして，縦長に整形しました。
m1 %>% tidyr::gather(審査員,val,-年代,-演者) %>% na.omit ->m1.long
m1.long$審査員 <- factor(m1.long$審査員)
m1.long$演者 <- factor(m1.long$演者)
r.name <- levels(m1.long$審査員)
p.name <- levels(m1.long$演者)
datastan <- list(L=nrow(m1.long),
N=max(as.numeric(m1.long$演者)),
M=max(as.numeric(m1.long$審査員)),
idX=as.numeric(m1.long$演者),
idY=as.numeric(m1.long$審査員),
X=m1.long$val)
モデルは先ほどと同様で，漫才師(演者)$$k$$の真のお笑いの実力$$\theta_k$$を持っていたとして，審査員(評定者)$$j$$がつけた得点$$X_{jk}$$は$$X_{jk}\sim N(\theta_k,\sigma_j)$$としています。
data{
int<lower=1> L; //data Length
int<lower=1> N; //number of players
int<lower=1> M; //number of rators
int idX[L]; //player ID index
int idY[L]; //rator ID index
real X[L]; // scores
}
parameters{
real<lower=0> theta[N];
real<lower=0> sig[M];
}
model{
for(l in 1:L){
X[l] ~ normal(theta[idX[l]],sig[idY[l]]);
}
theta ~ normal(50,100);
sig ~ cauchy(0,5);
}
generated quantities{
real PredScore[N,M];
for(n in 1:N){
for(m in 1:M){
PredScore[n,m] = normal_rng(theta[n],sig[m]);
}
}
}
これで推定した結果をグラフにします。まずは各演者の「漫才力$$\theta_k$$」のプロット。 点がEAP推定値，横幅は50％と95％のHDIです。
かわいそうなDonDokoDon。っていうか皆さん覚えてます？ぐっさんと平畠のコンビなんですが。
ただ，このモデルの中に含まれていることは，「評定者のスコアは100点満点で絶対的な感覚を持って採点している」ということ。だんだんスコアがインフレしているかも・・・というのは仮定していません。初期のメンバーはスコアが低くなりがちですが，ひょっとしたら審査員もみんな慣れてきて，スコアが上がってきたとかってことがあるかもです。
次に審査員の$$\sigma_j$$の結果。これは歪んでいるのでMED推定値にしています。 
立川談志師匠はまあいいとして，意外と大吉先生は評価の幅が大きい。一番小さいのは哲夫（笑い飯）ですが。審査員歴が一番長い松本人志は5点弱のブレのようですね。
Who are king of kings ?
事後分布から
さて，これらのスコアを元に，もし歴代のM-1選手たちが一堂に会し，実力一発勝負をしたらどうなるでしょうか? MCMCサンプルはバーンイン期間を除いて10000点とりましたので，$$\theta_k$$の同時確率空間における順序で一位を取った回数をグラフにしました。それがこちら。
# 順位 colnames(player) <- p.name player %>% dplyr::mutate(Winner=as.factor(p.name[apply(player,1,which.max)])) %>% ggplot(aes(x=Winner,fill=Winner)) + stat_count()+ coord_flip()+ theme(legend.position = 'none')
強いのはブラマヨ。ね。
事後予測分布から
違う角度から検討します。実力の通りではなく，評定者との関係というのもあるでしょうから，事後予測分布を使った予測をして見ます。
全てのプレイヤーが，全ての評定者に評価されたと考え，事後予測分布の平均スコアをもとに予測される得点順に並べたのがこちら！
pred <- rstan::extract(fit,pars="PredScore")$PredScore
predAverage <- data.frame(PPDscore=round(apply(pred,2,mean)))
predAverage$player <- factor(p.name)
predAverage %>% ggplot(aes(x=reorder(x=player,PPDscore,identity),y=PPDscore,fill=player))+
geom_bar(stat="identity") + coord_flip() +
theme(legend.position = 'none') +
ylab("Median Scores from PPD") + xlab("演者")
1位はブラックマヨネーズ！
2位はパンクブーブーとミキが同率！
3位はアンタッチャブル，オードリー，サンドウィッチマン，かまいたち，とろサーモン，トレンディエンジェルが並ぶ！
という結果になりますた。 ちなみに，2017年の結果を入れなければ，ミキは入ってこないし（初挑戦だから当然），ジャルジャルは3位に入れてたんですよw
もちろん，順番の問題もあるでしょう。また複数回ファイナリストになった人の実力はひとつ（成長や変化がない）というモデルだという限界はあります。
なにより，これはあくまでも確率です。実現値は常に一つ。可能性は色々あるんだけど，結果は結果。そこに生まれる物語や感動について，ベイジアンは言葉を持たず，ただ感動に浸り，漫才という芸術に感じ入るしかないものだと思います。
Enjoy Bayes and OWARAI!
追伸 博多大吉先生の評価についてのコメント，ラジオクラウドで聴けます。期間限定かも。是非〜 ＞ 博多大吉のM-1審査員をやって・・・ https://radiocloud.jp/archive/tama954/?content_id=24438
2件のコメント
コメントは受け付けていません。