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

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

調査の結果、アニメキャラのバストサイズとPixivでR-18タグが付く割合の相関係数は0.286・相関検定の結果はp=0.14であったため、バストサイズとR-18タグがつく割合には相関が無いという結論になりました pic.twitter.com/16bt8tyk4p

— 川奈 清 (@kivantium) 2015, 5月 12

同一作品の間では確かに相関があるっぽかったけど、作品間による差の方が大きかったため相関が消えてしまった

— 川奈 清 (@kivantium) 2015, 5月 12

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

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

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

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

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

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

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

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

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

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

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

これらを推定するには、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 3（core 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サンプリングだけをしています。

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

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
b0	-2.321	0.013	0.297	-2.906	-2.496	-2.334	-2.154	-1.690	528	1.01
bB	0.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_rA	0.810	0.008	0.278	0.448	0.618	0.755	0.931	1.482	1090	1.00
s_rB	0.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を見てみましょう。

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

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

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

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