声優の主役力をRstanでやる

MikuHatsune2015-05-30

声優の主役力の推定をモデル化してRstanでやってみる話の導入。
声優統計第五号に、主役力 : キャストの表記順に着目したプレイヤーレーティング、という論文があるが、これはTrue skill というモデル(これとかこれとか)を用いている。1対1のペアを作るが、同一アニメでのキャストは上から重要人物が挙げられているだろうという仮説に基づく。この記事で言えば、.lainによると

順位 アニメ役 声優
1 大宮忍 西明日香
2 アリス・カータレット 田中真奈美
3 小路綾 種田梨沙
4 猪熊陽子 内山夕実
5 九条カレン 東山奈央
6 大宮勇 田村ゆかり
7 烏丸さくら(烏丸先生) 佐藤聡美
8 松原穂乃花 諏訪彩花
9 忍のお母さん 高橋美佳子
10 以下省略 数十人いる

となっている。先の論文では
西明日香は下位のキャストにすべて勝っている
田中真奈美は下位のキャストにすべて勝っている
・・・
と繰り返している。これをモデルとしている。
 
今回は、1クールあたりに何本のアニメに出演しているかと、上位の役をどれくらい取っているかをそれとなくモデル化する。
とりあえず、上位の役をどれくらい取っているかだけを簡単にモデル化してみる。
 
そもそもの解析のモチベーションとして、大好きなのに主役を演じる機会があまりないような声優がいるような気がしている。というわけで今回の解析に登場する声優を軽く紹介すると
種田梨沙:別に主役が少ないというわけではなく、好きだから挙げた。うれシード。
小松未可子凪のあすからの美海やニセコイのつぐみといった、うーんまあ二番手なのはしょうがないよねなキャラ、かと思ったらモーレツ宇宙海賊のマリカとかKのネコとか、一番の主役キャラおったな、という感じ。でも基本的に二番手な印象。
内山夕実きんいろモザイクの陽子はサイコーである。Aチャンネルのナギも好きだったが、「内山夕実さん主役!!これ見なきゃ!!(使命感」となった記憶がない…
石原夏織:夏待ちの青である。変猫の小豆梓である。これだけでもう十分な滑り台感があるが、よくよく思い出してみると淫乱シュークリームでもあった。
喜多村英梨:まどまぎの青により元祖滑り台?? 主役キャラはたくさんあるのだがそれと同じくらいに二番手キャラも多い??
 
という5人をリストアップした。
rstan の練習を兼ねるので、ごく簡単なモデルにする。
上記の表で、種田梨沙を例にとると、きんいろモザイク(1期)には総勢44人(男性10人、女性34人。重複は除外した)の声優が出演する。この内、種田梨沙は女性キャラで3番めに挙げられている。種田梨沙は3番目に挙げられているので、女性キャラでは3番目、逆順では32番目ということにする。
種田梨沙が出演しているアニメは29本、出演順位は逆順でこんな感じ。

N_rank <- c(8, 7, 23, 4, 4, 6, 6, 4, 4, 5, 32, 3, 5, 6, 9, 3, 9, 8, 5, 14, 3, 2, 14, 18, 13, 7, 1, 4, 13)
standata <- list(N_anime=length(N_rank), N_rank=N_rank)

モデルとしては、キャスト順位はポアソン分布によるとして、
N_{rank}\sim poisson(p_{main})
でサンプリングされるとする。
このモデルはこちら。

# castmodel.stan
data{
	int<lower = 1> N_anime;
	int<lower = 1> N_rank[N_anime];
}
parameters{
	real p_main;
}
model{
	for(n in 1:N_anime)
		N_rank[n] ~ poisson(p_main);		
}

種田梨沙の収束の様子はこんな感じ。

 
2000年以降にデビューし、計20本以上のアニメに出演した声優について計算すると、結果は下の通り。
主役力は
石原夏織 > 喜多村英梨 > 内山夕実 > 種田梨沙 > 小松未可子
となる。
ただし、このモデルでは、総出演声優が多いときに上位を取ったアニメの影響がでかすぎると思われる。

cvs <- names(debut)[first$debut_year>=2000 & sapply(debut, sum)>=20]
res <- matrix(0, length(cvs), 3000)
rownames(res) <- cvs
for(tmp_cv in cvs){
	standata <- list(N_anime=sum(dat1$cv==tmp_cv), N_rank=dat1$rev[dat1$cv==tmp_cv])
	fit <- sampling(stanmodel, data=standata, chains=3, warmup=200, iter=1200, seed=123)
	res[tmp_cv,] <- extract(fit)[[1]]
}
m <- t(res)
m <- m[, order(colMeans(m))]
cvs <- c("種田梨沙", "小松未可子", "内山夕実", "石原夏織", "喜多村英梨")
idx <- match(cvs, colnames(m))
cols <- replace(rep(0, ncol(m)), idx, 2)
par(mar=c(5.5, 4, 2, 2))
b0 <- boxplot(m, cex=0.3, pch=16, xaxt="n", col=cols, ylab="Expectation of main cast")
axis(1, at=match(cvs, colnames(m)), labels=cvs, las=2)
  種田梨沙 小松未可子   内山夕実   石原夏織 喜多村英梨 
  8.327675   8.107743   9.810561  11.556734  10.330768 


 
というわけでモデルを修正する。
ある声優がN役キャストがいるアニメで、k番目の役を取ったすると、その人は\frac{k}{N}の主役力があるとする。こうすると[0,\hspace{3}1]にスケーリングされる。
ここで、いい感じに[0,\hspace{3}1]でサンプリングできそうな分布を探すと、ベータ分布が良さそうである。
0以上の実数\alpha, \betaによって決まり、期待値は\frac{\alpha}{\alpha+\beta}である。これをモデル化しよう。
上記の N_rank に対する、女性声優の出演人数は

n_female <- c(25, 19, 25, 4, 19, 8, 8, 9, 4, 8, 34, 11, 5, 16, 27, 7, 16, 22, 7, 15, 27, 2, 14, 45, 24, 8, 15, 26, 22)

モデルとしては、主役力はベータ分布によるとして、
N_{rank}\sim beta(\alpha,\hspace{5}\beta)
でサンプリングされるとする。
このモデルはこちら。

# castmodel2.stan
data{
	int<lower = 1> N_anime;
	real<lower=0, upper=1> N_rank[N_anime];
}
parameters{
	real<lower=0> alpha;
	real<lower=0> beta;
}
model{
	for(n in 1:N_anime)
		N_rank[n] ~ beta(alpha, beta);		
}

で、これをやってみるのだが

Error : 

Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

error occurred during calling the sampler; sampling not done

!?!?!?!?!?!?
とエラーになる。
stan-referenceの Initialization によると、初期値は[-2, 2]の一様分布からサンプリングされるけど、制限がない場合にはうんぬんかんぬん書いてある。
初期値の対応としてレプリカ交換法があるらしいが、とりあえずrstanを動かすために\beta=2に決め打ちしたら動いた。

# castmodel3.stan
data{
	int<lower = 1> N_anime;
	real<lower=0, upper=1> N_rank[N_anime];
	real<lower=0> beta;
}
parameters{
	real<lower=0> alpha;
}
model{
	for(n in 1:N_anime)
		N_rank[n] ~ beta(alpha, beta);		
}
standata <- list(N_anime=length(N_rank), N_rank=N_rank/n_female, beta=2)
fit <- sampling(stanmodel, data=standata, chains=3, warmup=200, iter=1200, seed=123)


 
このモデルでの主役力はこんな感じ。これは昇順にしてある。
石原夏織はこのモデルでも上位にきたが、これはおそらくマギでのマギ役が女性声優で1番に挙げられているし、思い返せばジャージ部とかあったなと。
喜多村英梨はやはりメインヒロイン役が多いんじゃないかという印象。ただし、2番手のキャラがたくさんあってもこのモデルだと主役力は1に近くなりうる。
小松未可子種田梨沙は…ようわからん。種田梨沙はモブキャラ出演も多いからこうなっている?
内山夕実はこの中では際立って主役力が落ちてしまった。内山夕実メインヒロインのラブコメに期待したい。

  内山夕実   種田梨沙 小松未可子 喜多村英梨   石原夏織 
 0.4834522  0.5441191  0.5737306  0.6281150  0.6490718 
cvs <- names(debut)[first$debut_year>=2000 & sapply(debut, sum)>=20]
res <- matrix(0, length(cvs), 3000)
rownames(res) <- cvs
beta <- 2
for(tmp_cv in cvs){
	s <- ifelse(dat1$sex[dat1$cv==tmp_cv][1]=="女性", "n_female", "n_male")
	standata <- list(N_anime=sum(dat1$cv==tmp_cv), N_rank=dat1$rev[dat1$cv==tmp_cv]/dat1[dat1$cv==tmp_cv, s], beta=beta)
	fit <- sampling(stanmodel, data=standata, chains=3, warmup=200, iter=1200, seed=123)
	res[tmp_cv,] <- extract(fit)[[1]]
}

m <- t(res/(res+beta))
m <- m[, order(colMeans(m))]]
cvs <- c("種田梨沙", "小松未可子", "内山夕実", "石原夏織", "喜多村英梨")
idx <- match(cvs, colnames(m))
cols <- replace(rep(0, ncol(m)), idx, 2)
par(mar=c(5.5, 4, 2, 2))
b0 <- boxplot(m, cex=0.3, pch=16, xaxt="n", col=cols, ylab="Expectation of main cast")
axis(1, at=match(cvs, colnames(m)), labels=cvs, las=2)