KOC2017の分析
はじめに
キングオブコント，2017が先日開催されましたね。いや，今年も面白かった。
にゃんこスター，すごかったなあ，一本目。あんなんずるいやん，と言いながらずっとゲラゲラ笑ってました。お笑いの面白いところは，ルールがあるようでないとか，ある程度方向性が決まったらそれを崩そうという勢力がでてきて，常に例外・規格外・問題外みたいなのが生まれてくる土壌があるところだと思っています。審査員長の松本人志がいう，「ある程度いやらしさが必要」というのはそうところだと思うの。そんなの卑怯じゃない，そんなのマナー違反じゃない，と思いながら笑わずに得られないというかね。あかんと思いながら笑ってしまっている自分がいるというかね。
でもまあ，2本同じ系列だったのは残念なところ。かまいたち・山内に言われなくても「一回冷静に」なっちゃったもんね。違う切り口で笑わせてもらいたいと思うわけですよ。流石にあのやり方，二回連続では笑わない。そういう意味では，かつての2700はよかった。右肘左肘の次にやった「キリンスマッシュ・キリンレシーブ」のネタってみなさん覚えてます？あれすごいことですよ。色物すぎて賞は取れなかったけどねえ。
今回この記事を書こうと思って，審査員のスコアをググって見たんですが，キーワード「キングオブコント 採点」で調べると，やれ松本人志の基準がぶれているだの，アンガールズに甘かっただのといった意見があるわけですが，個人的には我々が共感できないレベルでの評価軸を持っててこその審査員だと思うので，私はそういう反論に組みするものではないです。ただ見て楽しい人と，舞台に立って実際にやる人，やってる経験がある人から考える上手い下手ってあると思うもんね。もちろん好みもね。
ということで，審査員の好みとか，にゃんこスターの特殊性みたいなのがデータから観れるかどうか，分析してみます。
まとめ
だらだら読むと大変なので，先に結論を。
- 相関係数をベイズ推定するときはサンプルサイズをたっぷりね
- ウィッシャートブンプを使ったらその辺すこしはましだわね
- MDSとMCMCの組み合わせで幻想的な花火が描けますよ
ここから本題。
library(rstan)## Loading required package: ggplot2## Loading required package: StanHeaders## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
dat <- matrix(
c(87, 89, 88, 87, 83,
93, 92, 93, 89, 85,
90, 93, 95, 93, 93,
89, 90, 92, 92, 89,
85, 86, 85, 85, 80,
90, 93, 87, 90, 95,
90, 89, 97, 93, 97,
90, 82, 83, 87, 90,
86, 85, 82, 84, 82,
85, 85, 82, 86, 84,
93, 93, 87, 89, 90,
90, 92, 90, 93, 93,
95, 97, 92, 91, 92,
94, 95, 95, 97, 96,
90, 92, 92, 94, 94),ncol=5,byrow=T)まずは単純な審査者の評定の類似性。まあこれは相関係数でいいと思うのです。 ということで，しれっとstanで推定してみます。
stancode1 <-"
data{
int<lower=0> N; //number of raters
int<lower=0> M; //number of actors
vector[N] Y[M];
}
parameters{
vector[N] mu;
vector<lower=0>[N] sig;
corr_matrix[N] rho;
}
transformed parameters{
cov_matrix[N] V;
V = quad_form_diag(rho,sig);
}
model{
for(i in 1:M)
Y[i] ~ multi_normal(mu,V);
sig ~ student_t(4,0,5);
rho ~ lkj_corr(1);
}
"
# stanによる相関係数の推定
model1 <- stan_model(model_code=stancode1,model_name="cor estimation")
datastan <- list(N=ncol(dat),M=nrow(dat),Y=as.matrix(dat))
fit <- sampling(model1,datastan)transformed parametersブロックのところでquad_form_diagという謎関数を使っていますが，これはStanマニュアルを読んでいただければ「あー，行列とベクトルを組み合わせた計算するのね」とお分りいただけると思います。
あとは素直にコード化しているので，読んでいただければわかるかなと。 で，相関係数はちょっと歪みがあるかもしれないので，get_posterior_meanで取り出すのではなくmedian(MAP)推定値が欲しい。 ということで，MCMCサンプルから中央値を取り出して行列にするコードをR側で準備。
# サンプルの取り出し
rhos <- rstan::extract(fit,pars="rho")
# 受け皿の準備
tmp <- matrix(nrow=ncol(dat),ncol=ncol(dat))
colnames(tmp) <- colnames(dat)
rownames(tmp) <- colnames(dat)
for(i in 1:(ncol(dat)-1)){
for(j in (i+1):ncol(dat)){
tmp[i,j] <- median(rhos$rho[,i,j])
tmp[j,i] <- tmp[i,j]
}
}
diag(tmp) <- 1
tmp## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.5725501 0.3563585 0.3797303 0.4199263
## [2,] 0.5725501 1.0000000 0.4650221 0.4943685 0.3092411
## [3,] 0.3563585 0.4650221 1.0000000 0.6881482 0.3889331
## [4,] 0.3797303 0.4943685 0.6881482 1.0000000 0.7111259
## [5,] 0.4199263 0.3092411 0.3889331 0.7111259 1.0000000これでできましたよっと。 念のために確認。普通の記述統計的推定値と近い値が得られていますかね？
cor(dat)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.7647386 0.6058435 0.6378767 0.6382824
## [2,] 0.7647386 1.0000000 0.6930632 0.7072440 0.5726141
## [3,] 0.6058435 0.6930632 1.0000000 0.8482605 0.6519889
## [4,] 0.6378767 0.7072440 0.8482605 1.0000000 0.8585755
## [5,] 0.6382824 0.5726141 0.6519889 0.8585755 1.0000000む，ちょっと違うな。なんでやろ。
ということでちょっとシミュレーションとかして確認してたんですが，サンプルサイズが小さい（今回15回分しかない）のが問題だったかも。irisデータとか使うと，まあまあな数字になるんですけど。それにしてもずれてるなあ。
とツイッターで呟いたら，どこからともなく(@simizu706$$)「偏差積和行列渡して，ウィッシャート分布使いなはれ。あと標準偏差の推定やめなはれ」というアドバイスが。ちょっとやってみましょう。(2017.12.04.コードのコメント修正)
stancode2 <-'
data{
int<lower=0>N; //number of raters
int<lower =0>M; //number of actors
matrix[N,N] spd; //偏差積和行列
}
transformed data{
vector[N] sig;
for(n in 1:N) sig[n] = (spd[n,n]/(M-1))^0.5;
}
parameters{
corr_matrix[N] rho;
}
transformed parameters{
cov_matrix[N] V;
V = quad_form_diag(rho,sig);
}
model{
rho ~ lkj_corr(1);
spd ~ wishart(M-1,V); //M-1は自由度
}
'
# stanによる相関係数の推定
model2 <- stan_model(model_code=stancode2,model_name="cor with wishart")
tmp <- matrix(as.numeric(scale(dat)),ncol=ncol(dat))
spd <- t(tmp) %*% tmp
datastan <- list(N=ncol(dat),M=nrow(dat),spd=spd)
fit2 <- sampling(model2,datastan)# サンプルの取り出し
rhos <- rstan::extract(fit2,pars="rho")
# 受け皿の準備
tmp <- matrix(nrow=ncol(dat),ncol=ncol(dat))
colnames(tmp) <- colnames(dat)
rownames(tmp) <- colnames(dat)
for(i in 1:(ncol(tmp)-1)){
for(j in (i+1):ncol(tmp)){
tmp[i,j] <- median(rhos$rho[,i,j])
tmp[j,i] <- tmp[i,j]
}
}
diag(tmp) <- 1
tmp## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.6026052 0.4030287 0.4139016 0.4544265
## [2,] 0.6026052 1.0000000 0.5073684 0.5239749 0.3529765
## [3,] 0.4030287 0.5073684 1.0000000 0.7338763 0.4616989
## [4,] 0.4139016 0.5239749 0.7338763 1.0000000 0.7528307
## [5,] 0.4544265 0.3529765 0.4616989 0.7528307 1.0000000ふぬぅ，まあいいでしょう。 これをみると，バナナマンの二人はやっぱり評価が似てるんですよね。松本さんと大竹さんも評価が似ている。
似てる，似てないからマップを作るのは，そうだね，MDSだね。
result.mds <- cmdscale((1-tmp),2)
plot(result.mds,type="n")
text(result.mds,colnames(tmp))
こうしてみると，左右はっきり分かれます。左側がバナナマン的な評価次元。設楽さんはわりとロジックの通ったシュールな世界とかを作るので，やや論理的な方向，右側は勢いを評価する次元というところかな。
上下で考えると，上はシュールで下はインパクトみたいな感じもしますが。
ちなみにこれはMAP推定値をプロットしているんであって，実際はたくさんのMCMCサンプルがある=座標も分布しているわけです。 ちょっとこの分布をプロットしてみましょう。
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.4-4configs <- data.frame()
for(i in 1:dim(rhos$rho)[1]){
# 各サンプルをMDS
tmp<- cmdscale(1-rhos$rho[i,,])
# プロットを整えるためにプロクラステス回転
ret.tmp <- data.frame(vegan::procrustes(result.mds,tmp)$Yrot)
ret.tmp$iter <- i
ret.tmp$rator <- factor(1:5,label=c("設楽","日村","三村","大竹","松本"))
configs <- rbind(configs,ret.tmp)
}
library(ggplot2)
# マカーの呪文
old = theme_set(theme_gray(base_family="HiraKakuProN-W3"))
# 散布図
g <- ggplot(configs,aes(x=X1,y=X2,color=rator))+geom_point(alpha=0.2)
g 
ぼんやり雲のように重ねてみました。こうしてみると，三村さんのポジションがちょっと広く散らばっているようですね。他の人と相対的にみると，ちょっと立ち位置の確信度が低い，という言い方ができるかもしれません。
相関行列に対してMDSをする，というのはない話ではないのですが，MCMCサンプルについてそういうことをやるっていうのは珍しいかな，と思ってブログにしてみました。こういう遊び方もありだと思うんです。
おまけ
最後におまけです。 行と列を入れ替えて，すなわち演者の方についてMDSをやってみましょう。公平を期すために？1回目のネタだけを対象にやってみたいと思います。自由度の関係でMCMCもできないから，記述統計量で。
# 行列逆転
dat <- t(dat[1:10,])
colnames(dat) <- c("わらふじなるお","ジャングルポケット","かまいたち","アンガールズ","パーパー","さらば青春の光",
"にゃんこスター","アキナ","GAG少年楽団","ゾフィー")
result.mds2 <- as.data.frame(cmdscale(1-cor(dat),2))
g <- ggplot(result.mds2,aes(x=V1,y=V2,label=colnames(dat))) + geom_point()+geom_text(family="HiraKakuProN-W3")
g
こうしてみると，やっぱりにゃんこスターは独特のところに位置づくんだな。うん。
おまけのおまけ
@simizu706$$ さんから，相関行列のMCMCサンプルから行列で推定値を取り出す方法，教えてもらいました。感謝。
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionfit2 %>% rstan::extract() %$% rho %>% apply(c(2,3),median)##
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.6026052 0.4030287 0.4139016 0.4544265
## [2,] 0.6026052 1.0000000 0.5073684 0.5239749 0.3529765
## [3,] 0.4030287 0.5073684 1.0000000 0.7338763 0.4616989
## [4,] 0.4139016 0.5239749 0.7338763 1.0000000 0.7528307
## [5,] 0.4544265 0.3529765 0.4616989 0.7528307 1.0000000