はじめに

この記事は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推定値にしています。

最強のM-1漫才師は誰だ

最強のM-1漫才師は誰だ

はじめに

結論；ブラマヨが一番すごいんじゃないか

手続き