読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRでベイズ統計モデリングします. たまに書評.

kivantiumさんのブログ記事「アニメキャラのバストサイズとPixiv R-18タグ率の関係」の追加解析

Stan

@kivantiumさんの以下のツイートが面白そうすぎて追加解析してみました。特に2つ目のツイートが重要で、これがないと階層ベイズでやってみようという気は起りませんでした。

kivantiumさんはブログ記事およびデータもすぐにアップして下さっています(アニメキャラのバストサイズとPixiv R-18タグ率の関係)。今回はこちらのデータをそのまま使わせていただきます。

解析する前にまずは散布図行列を見ます。描き方はこの記事参照。しかし、行列の下三角部分が物足りなくて面倒でしたが色々カスタマイズしました。

20150515_ggpairs1.png

クリックすると大きくなります。色は濃い虹色みたいな色はアニメの違い、薄い赤・青はヒロインでないか・ヒロインかを表します。

まずは対角線のヒストグラムから見ていきます。各アニメでキャラがどれくらいいるかとか、バストの分布は一様分布に近いけど貧乳は少なくて巨乳は多そうだなーとか、ratioの分布は右にtail引いてる感じだなーとか読み取ります。

次に右上を見ます。アニメを区別しないでのSpearmanの順位相関係数(✕100)なのですが、区別しないとどれも強い相関なさそうかな。強いて言えば、ヒロインはバストが大きくなさそうとか、バストとratioは正の相関あるかもぐらいです。

次に左下の図を見ます。特に重要なのは行列成分で言えば4行1列目の図。アニメごとにかなりratioが違う感じです。ratioの平均が大きいアニメと平均が小さいアニメを比べると、平均が大きいアニメの方が同一アニメ内のratioのバラつきが大きいような印象を受けます。次に4行3列目の図を見ます。同一アニメ内だとバストとratioは相関ありそうですが、その依存の傾きはアニメに依存してそうです。このデータの数からは先ほどのバラつきはどうもバスト依存性によるものと考えておくのが無難そうです。

ここまで来て、モデルを構築します。今回は作品数がかなり多いので二項分布を正規分布で近似してratioを直接モデリングしてもそこまで問題ないと思いますが、計算してみたらやはり予測分布のベイズ信頼区間の下限がゼロ以下になってしまうキャラが何人かいましたので、おとなしく二項分布回帰を使います。モデル式は以下になります。

20150515_image01.png

20150515_image02.png

20150515_image03.png

20150515_image04.png

ここでnはキャラの添え字で1~Nになります。inv_logitというのは1/(1+exp(-x))のことです。関数の形はwikipediaを参照して下さい。この関数は(-∞,∞)を[0,1]に変換します。20150515_image05.pngはR18になる確率を表します。このようにモデル化することで、inv_logitの中身が例えば1だとすると、オッズと呼ばれる量(=R18になる確率/R18にならない確率)がexp(1)≒2.7になると解釈できます。このあたりは久保緑本が詳しいので知りたい人はぜひそちらを読んで下さい。

20150515_image06.pngはキャラ共通の切片項、20150515_image07.pngはアニメの影響です。20150515_image08.pngはそのキャラがどのアニメに属しているかのアニメの添え字が返ります。20150515_image09.pngはキャラ共通のバストの影響(傾き)、20150515_image11.pngはアニメに依存するバストの影響(傾き)です。20150515_image10.pngはバストのデータが入ります。「-80」しておくことで20150515_image07.pngがバストを80に揃えた時のアニメの影響となって比較しやすいです。逆に「-80」がないとバストがゼロの時のアニメの影響となるのでそのままでは比較しづらくなります(他の値を代入していけば比較できる)。20150515_image12.pngはヒロインの影響(傾き)で20150515_image13.pngにはヒロイン情報のデータが入ります。

20150515_image14.pngにはR18作品数、20150515_image15.pngには総作品数が入ります。20150515_image15.pngと先ほどの20150515_image05.pngから二項分布による確率的な変動を込みにして最終的なR18作品数になると考えます。

各アニメのキャラが少ないので、アニメの影響やアニメに依存するバストの影響を制約なしで推定しようとすると失敗します。そこで階層ベイズの出番です。アニメの影響は平均ゼロ・標準偏差20150515_image16.pngに従うという緩い制約を入れます。20150515_image16.pngは無情報事前分布としてデータから推定します。同様にアニメに依存するバストの影響も平均ゼロ・標準偏差20150515_image17.pngに従うとします。

これらを推定するには、Stanと{rstan}パッケージを使うとラクで速いです。Stanコードは以下になります。

data {
   int<lower=0> N;
   int<lower=1> A;
   int<lower=1,upper=A> Anime[N];
   int<lower=0,upper=1> isHeroine[N];
   real<lower=0> B[N];
   int<lower=1> M[N];
   int<lower=1> Y[N];
}

parameters {
   real b0;
   real bB;
   real bH;
   real rA[A];
   real rB[A];
   real<lower=0> s_rA;
   real<lower=0> s_rB;
}

transformed parameters {
   real q[N];
   for (n in 1:N)
      q[n] <- inv_logit(b0 + rA[Anime[n]] + (bB + rB[Anime[n]])*(B[n] - 80) + bH*isHeroine[n]);
}

model {
   for (a in 1:A){
      rA[a] ~ normal(0, s_rA);
      rB[a] ~ normal(0, s_rB);
   }

   for (n in 1:N)
      Y[n] ~ binomial(M[n], q[n]);
}

generated quantities {
   real y_pred[N];
   for (n in 1:N)
      y_pred[n] <- binomial_rng(M[n], q[n]);
}

先ほどのモデル式とほとんど同じです。37-41行目は推定するだけなら必要ないのですが、予測分布を求めたいので追加しました。計算時間はSurface Pro 3core i5)で1chainあたりおよそ4秒でした。

キックするRコードは以下のとおりです。

library(rstan)

d <- read.csv('input/data.csv', header=TRUE, fileEncoding='UTF-8')
d$anime <- factor(d$anime, levels=unique(d$anime))

stanmodel <- stan_model(file='model.stan')
standata <- list(N=nrow(d), A=length(unique(d$anime)), Anime=as.numeric(d$anime), isHeroine=d$isHeroine, B=d$B, M=d$total, Y=d$R18)
fit <- sampling(stanmodel, data=standata, chains=3, warmup=200, iter=1200, seed=123)
  • 6行目はコンパイルだけをしています。
  • 8行目はMCMCサンプリングだけをしています。

結果は以下になります。 まずはパラメータの事後平均などから。

meanse_meansd2.5%25%50%75%97.5%n_effRhat
b0-2.321 0.013 0.297 -2.906 -2.496 -2.334 -2.154 -1.690 528 1.01
bB0.049 0.001 0.022 0.005 0.036 0.049 0.063 0.092 466 1.00
bH-0.092 0.001 0.030 -0.153 -0.112 -0.091 -0.072 -0.035 2551 1.00
rA[1]0.533 0.013 0.298 -0.103 0.365 0.543 0.711 1.120 528 1.01
rA[2]-0.665 0.013 0.297 -1.302 -0.830 -0.652 -0.487 -0.082 527 1.01
rA[3]0.392 0.013 0.298 -0.240 0.223 0.404 0.567 0.985 533 1.01
rA[4]-0.125 0.013 0.299 -0.772 -0.293 -0.112 0.055 0.468 529 1.01
rA[5]0.753 0.013 0.324 0.090 0.560 0.760 0.940 1.420 610 1.01
rA[6]-0.860 0.013 0.329 -1.550 -1.057 -0.849 -0.657 -0.238 688 1.01
rA[7]-0.732 0.013 0.303 -1.375 -0.910 -0.719 -0.551 -0.131 523 1.01
rA[8]0.669 0.017 0.525 -0.280 0.314 0.651 0.999 1.716 942 1.00
rA[9]0.121 0.014 0.510 -0.904 -0.196 0.116 0.444 1.094 1373 1.00
rB[1]0.017 0.001 0.022 -0.026 0.002 0.016 0.030 0.062 487 1.00
rB[2]-0.005 0.001 0.022 -0.048 -0.019 -0.005 0.009 0.040 467 1.00
rB[3]0.079 0.001 0.024 0.034 0.064 0.078 0.094 0.127 526 1.00
rB[4]-0.019 0.001 0.022 -0.062 -0.033 -0.019 -0.006 0.026 477 1.00
rB[5]-0.015 0.001 0.026 -0.065 -0.031 -0.015 0.001 0.035 472 1.00
rB[6]0.031 0.001 0.032 -0.028 0.009 0.030 0.052 0.094 1008 1.00
rB[7]-0.064 0.001 0.023 -0.109 -0.078 -0.064 -0.049 -0.018 492 1.00
rB[8]-0.028 0.001 0.056 -0.153 -0.060 -0.025 0.007 0.073 1703 1.00
rB[9]0.011 0.001 0.049 -0.082 -0.020 0.009 0.040 0.111 1718 1.00
s_rA0.810 0.008 0.278 0.448 0.618 0.755 0.931 1.482 1090 1.00
s_rB0.058 0.001 0.022 0.030 0.042 0.053 0.067 0.113 1107 1.00
.................................
lp__-43842.205 0.127 3.563 -43850.350 -43844.259 -43841.812 -43839.713 -43836.209 790 1.00

bBは共通のバストの影響(1cmあたりの傾き)ですが、ゼロより大きい確率は0.984でしたのでたしかに影響はあると言えそうです。しかしその影響度合いは5cmあたりで考えてみても、0.049*5=0.245でアニメの影響の標準偏差であるs_rAの値(0.810)よりだいぶ小さそうです。ヒロインの影響もゼロより小さい確率は0.9993でしたのでこちらも影響があると言えそうです。ヒロインはR18になりにくいと言えるということです。しかしその影響度合いはやはり小さそうです。

アニメ別にrA, rBを見てみましょう。

20150515_violinplot_rA-2.png

色はアニメの違い、点はMCMCサンプルの中央値、線はMCMCサンプルから求めた80%ベイズ信頼区間、分布はMCMCサンプルの分布です。散布図行列で見たようなアニメの影響が出ています。ただしキャラが1人しかいないとやはり推定の幅は広くなります。

20150515_violinplot_rB-2.png

こちらはアニメに依存するバストの影響(傾き)です。点や線はさきほどの図と同じです。『ソードアートオンライン』や『冴えない彼女の育て方』がバストの影響が強そうです。

最後に、実際のR18作品数とこのモデルから予測されたR18作品数の80%予測区間を比べました(クリックすると大きくなります)。

20150515_visual_predict_check-2.png

点と線は前の図と同様です。『ラブライブ』の一部のキャラ以外はかなりうまく推定できているようです。しかし、予測よりもR18作品数が少ない『小泉花陽』さんと『星空凛』さん、予測よりもR18作品数が多い『東條希』さんと『矢澤にこ』さんあたりが目立ちます。このモデルでは表現しきれていない部分があるのでしょう。僕は現在のアニメの知識があまりないので、このあたりの考察はドメイン知識を豊富に持つ識者にゆだねます。どんなデータを追加で取ればよさそうかを考えるきっかけになるかもしれません。kivantiumさん、データ公開や考察など大変ありがとうございました!