廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[RStan]ディリクレ・多項モデルで μ's とAqours の人気の差を調べる

背景の整理

前に書いたエントリが引用されて嬉しかったのでもう一ネタ書きます.

モデル

キャラクター i の得票数 y_i はパラメータ p_i の多項分布に従うとします。

ここで p_i は、以下のディリクレ分布に従うとします。

\displaystyle p_i  \sim \mathrm{Dirichlet}(\exp(\beta_0 + \beta_1 x_i))

で、x_i はキャラクター i が μ's に所属しているとき 1、さもなくば 0 の値を取る変数とします。

係数 \beta_1 はグループが人気に与える効果の強さを表すと解釈できます。

Stan で書くとこうなりました。

data{
  int<lower=1> M;
  int<lower=1> V[M];
  matrix[M,2] X;
}
parameters{
  simplex[M] theta;
  vector[2] beta;
}
model{
  theta ~ dirichlet(exp(X*beta));
  V ~ multinomial(theta);
}

解析対象

こんなデータです。

f:id:abrahamcow:20170624124125p:plain

横軸が得票数。

結果

Stan で MCMC を2000回して収束しました。

サンプリングは一瞬で終わります。

トレースプロットはこんな感じです。

f:id:abrahamcow:20190912164553p:plain

当てはまりを見るために最尤推定値(赤)と提案モデルによる推定値(青)をプロットしてみます。

f:id:abrahamcow:20190912165110p:plain

次にグループ効果 \exp(\beta_1) の事後分布を見ます。

f:id:abrahamcow:20190912164837p:plain

青い実線が事後期待値、点線が95%信用区間です。

μ's は Aqours にくらべ 2倍以上人気であるといえそうです。

付録

R のコードを貼り付けます。

library(rstan)
library(cowplot)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
pop <-"name vote group
矢澤にこ 384 muse
西木野真姫 384 muse
南ことり 370 muse
東條希 336 muse
綾瀬絵里 283 muse
園田海未 263 muse
小泉花陽 251 muse
星空凛 226 muse
津島善子 190 aqours
黒澤ルビィ 171 aqours
国木田花丸 155 aqours
渡辺曜 141 aqours
高坂穂乃果 128 muse
松浦果南 110 aqours
黒澤ダイヤ 103 aqours
桜内梨子 93 aqours
小原鞠莉 80 aqours
高海千歌 74 aqours"
pop_dat <-read.table(text = pop,header = TRUE)
ggplot(pop_dat,aes(x=reorder(name,vote),y=vote,fill=group))+
  geom_col()+
  coord_flip()+
  xlab("")+
  theme_cowplot(font_family = "HiraKakuProN-W3")
pop_dat4stan <- list(M=nrow(pop_dat),V=pop_dat$vote,X=model.matrix(~group,data=pop_dat))
lovelive <-stan_model("~/Documents/lovelive.stan")
smp_pop <-sampling(lovelive,pop_dat4stan)
traceplot(smp_pop,pars=c("beta","theta"))
all(summary(smp_pop)$summary[,"Rhat"]<1.1)
ex_pop <-extract(smp_pop)
pop_dat$est <- get_posterior_mean(smp_pop,pars="theta")[,"mean-all chains"]
pop_dat$lower <- apply(ex_pop$theta,2,function(x)quantile(x,probs = 0.025))
pop_dat$upper <- apply(ex_pop$theta,2,function(x)quantile(x,probs = 0.975))

ggplot(pop_dat,aes(x=reorder(name,vote)))+
  geom_linerange(aes(ymin=lower,ymax=upper,colour="pred"))+
  geom_point(aes(y=est,colour="pred"),size=3,alpha=0.7)+
  geom_point(aes(y=vote/sum(vote),colour="obs"),size=3,alpha=0.7)+
  coord_flip()+
  xlab("")+
  theme_cowplot(font_family = "HiraKakuProN-W3")
hist(exp(ex_pop$beta[,2]),breaks = "scott",xlab = "",main = "")
abline(v=quantile(exp(ex_pop$beta[,2]),probs = 0.025),col="royalblue",lwd=2,lty=2)
abline(v=quantile(exp(ex_pop$beta[,2]),probs = 0.975),col="royalblue",lwd=2,lty=2)
abline(v=mean(exp(ex_pop$beta[,2])),col="royalblue",lwd=2)