種田梨沙にはまっている。最近ではご注文はうさぎですかのリゼのバイトする姿やお着替えシーンでブヒったり、極黒のブリュンヒルデの寧子のおっぱいのほくろがどうなってんのか小一時間観察したりと楽しい日々である。
そんな種田梨沙だが、ゆゆ式の縁あたりから、なんとなく「百合アニメ」に出演が多いような気がするのが最近気になっているところである。というのも、きんいろモザイクの綾とか、そもそも最初の主役級キャラである神世界よりの早季がﾎﾞﾉﾎﾞｫ…してたとか、割れたりんごでファンになったわとかいろいろあるけど、種田梨沙が出ると百合アニメになりがちじゃないかと思っている。
で、百合アニメと種田梨沙の出演歴から、本当に偏りがあるかを検定して確かめようと思ったわけだが、いまから種田梨沙に張り付いて、数あるアニメのオーディション時にランダム割付しよう!!…というRCTは無理で、今までの出演歴をデータベースから漁ってデータを取るしかない。しかし、こういうデータを漁ると、例えば日常系アニメが多いまんがタイムきらら系列のアニメならば百合アニメになる確率は高そうだし、特定の音響監督に偏りがちという現象が既に報告されているように、単純にアニメ出演歴と百合アニメの関係をとると、バイアスの影響をもろに受ける(音響監督についてはキャスティング情報のbag-of-声優モデルを用いた音響監督推定問題として既に取り組まれている)。ということで、いままで勉強してきたPropensity scoreを使って、出演情報の裏に隠れている、監督や制作会社などの変数を調整して、種田梨沙が百合アニメに出やすいのかどうか考える。
結論としては偏って出演していそうだけども、統計学的に有意な結果とまではならなかった。
.lainというアニメ・声優データベースからデータをパースして、種田梨沙の出演が急増した2012年以降のTVアニメデータを対象とする。
百合アニメかどうかの判定は、自分と師匠のふたりが独立に評価し、どちらか一方が百合アニメと判定した場合に百合アニメとした。結果、403アニメ中92アニメが百合アニメということにした。
パラメータの数が多くなりすぎてもアレなので、アニメの出演声優は先頭からクレジットされている声優5人まで、403作品に7回以上出演している女性声優、アニメ制作関係者を対象にして、女性声優48人、スタッフ50人までにした。
PSを行うときに、(アウトカム達成症例数) / (独立変数)8 だとバイアスを生じやすくなるらしいので、まあこんなもんか。
補正をしない状態でフィッシャーの正確確率検定を行うと、
| 百合 | 非百合 | |
| 種田(+) | 9 | 16 |
| 種田(-) | 83 | 295 |
OR=2 くらいだがp=0.1367 で百合アニメに出やすそうとは言えないっぽい。
データが用意できたので、アニメ制作スタッフから種田梨沙が起用される確率(PS)を計算し、マッチングを行うと、種田梨沙が出ている百合アニメ(と判定したもの)が25アニメしかなかったので、マッチング後に21ペアが残った。
| 百合 | |||
| 種田(+) | 種田(-) | ||
| 種田(+) | 1 | 3 | |
| 非百合 | 種田(-) | 7 | 10 |
McNemarの正確確率検定でOR=2.3, p=0.3438 とこれまた統計学的有意差は見られなかった。
というわけで、ORが2と2倍くらい百合アニメになりそうな予感はして、いままでの感覚と大きく外れてはいなさそうなのだけれども、統計学的には有意というわけではなかったようだ。
データの収集
# Python import os import re import time import progressbar wd = "/taneda/" files = os.listdir(wd) staff = re.compile("<title>.* - スタッフ - \.lain</title>") a = re.compile('">.+?</a>') title = " - アニメデータベース" cast = "キャスト" media = "<dt>メディア</dt>" staffs = 'href="/mediadb/staff/view/' hosobi = "<dt>放送日</dt>" onair = re.compile("\d+年\d+月") cvre = re.compile('<a href="/voicedb/profile/.*">.*</a>') widgets = ["progress:", Percentage(), Bar()] pbar = ProgressBar(maxval=len(files), widgets=widgets).start() w0 = open("cvkyoen.txt", "w") w0.write("\t".join(["title", "media", "cv", "staff", "year", "month", "day"]) + "\n") for f in range(len(files)): #pbar.update(pbar.currval + 1) g = open(wd + files[f], "rU") cvlist = [] res2 = "" staff_res = "" onair_res = "" for i in range(1000): tmp = g.readline() if title in tmp: # 作品タイトル res1 = tmp.strip().split(title)[0].split("<title>")[-1] elif media in tmp: # 放送メディア res2 = g.readline().strip()[4:-5] elif len(cvre.findall(tmp)) > 0: # 声優 cvlist += [cvre.findall(tmp)[0].split("</a>")[0].split(">")[-1]] elif staffs in tmp: # 制作スタッフ staff_res = map(lambda x: x[2:-4], a.findall(tmp)) elif hosobi in tmp: # 放送開始日 for j in range(3): tmp = g.readline() if len(onair.findall(tmp)) > 0: onair_res = onair.findall(tmp)[0] if res2 == "TVアニメ" and len(onair_res) > 0: w0.write(res1 + "\t" + res2 + "\t" + ",".join(cvlist) + "\t" + ",".join(staff_res)+ "\t" + "\t".join(map(lambda x: re.split("\D", onair_res)[x], [0, 3, 6])) + "\n") w0.close()
library(arules) library(rbounds) library(exact2x2) dat <- read.delim("cvkyoen1.txt") # ニコニコ動画から 女性声優一覧 を取っておく cvlist <- as.character(read.csv("cv_women.txt", header=FALSE)$V1) # 先頭何人の声優を考えるかと、出演が何回以上かの制限 head_n <- 5 cv_n <- 7 staff_n <- 7 cv <- lapply(strsplit(as.character(dat$cv), ","), unique) taneda <- mapply(function(x) "種田梨沙" %in% x, cv) + 0 cv <- lapply(cv, head, head_n) cv <- as(as(cv, "transactions"), "matrix") cv <- cv[, colSums(cv) >= cv_n] cv <- cv[, !is.na(match(colnames(cv), cvlist))] fisher.exact(table(taneda, yuri)) staff <- lapply(strsplit(as.character(dat$staff), ","), unique) staff <- as(as(staff, "transactions"), "matrix") staff <- staff[, colSums(staff) >= staff_n] yuri <- (dat$yuri | dat$yuri2) + 0 dat1 <- as.data.frame(cbind(yuri, cv, staff)) dat1 <- dat1[, !duplicated(colnames(dat1))] # 出演回数の頻度 par(mfrow=c(1, 2)) barplot(table(colSums(cv)), main="女性声優の出演アニメ数") barplot(table(colSums(staff)), main="アニメ制作関係者の出演数") # クール別の百合アニメの本数の状況 season <- cut(dat$month, c(1, 3, 6, 9, 12), include.lowest=TRUE) res <- array(0, c(3, 4, 2)) dimnames(res) <- list(rev(unique(dat$year)), c("冬", "春", "夏", "秋"), c("no", "yes")) for(i in seq(3)){ idx <- dat$year == 2011 + i tmp <- tapply(yuri[idx], season[idx], table) res[i, , ] <- matrix(unlist(tmp), nc=2, byrow=TRUE) } tmp <- t(cbind(c(t(res[,,1])), c(t(res[,,2])))) par(mar=c(4, 4.5, 2, 2), xpd=TRUE) b <- barplot(tmp, las=1, col=grey(c(0.9, 0.7)), ylab="アニメ本数", cex.lab=1.6, cex.axis=1.2) p <- par()$usr text(b, 0, c("冬", "春", "夏", "秋"), adj=c(NA, 2), cex=1.2) x <- tapply(b, rep(1:3, each=4), range) for(i in seq(x)){ segments(x[[i]][1], -6.7, x[[i]][2], lwd=2)} text(sapply(x, mean), 0, paste(2012:2014, "年"), adj=c(NA, 4.5), cex=1.2) up <- 2 text(b, up, tmp[2,]) text(b, tmp[1,]+up, tmp[1,]) legend("topright", legend=c("非百合", "百合"), cex=1.5, pch=22, pt.bg=grey(c(0.9, 0.5)), bty="n") # 種田梨沙が起用される確率 g1 <- glm(taneda ~ staff, family = binomial(logit)) D <- taneda # 治療する Y <- yuri # 結果 X <- g1$fitted.values # PS dens <- mapply(function(x) density(X[D==x]), 0:1, SIMPLIFY=FALSE) m.obj <- Match(Y = Y, Tr = D, X = X, M = 1, replace=FALSE, estimand = "ATT", caliper=1) m <- m.obj$MatchLoopC D1 <- D[m[, 1:2]] Y1 <- Y[m[, 1:2]] X1 <- X[m[, 1:2]] # PS dens1 <- mapply(function(x) density(X1[D1==x]), 0:1, SIMPLIFY=FALSE) xl <- c(-1, 1)*max(dens[[1]]$y, dens[[2]]$y, dens1[[1]]$y, dens1[[2]]$y) yl <- c(0, 1) par(mar=c(5, 4.5, 2, 2)) plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="Propensity score (PS)", cex.lab=1.6, cex.axis=1.6) abline(v=0, lty=2) x <- length(X1)/length(X) points(dens[[2]]$y, dens[[2]]$x, type="l", col="red", lwd=2, lty=3) points(-dens[[1]]$y, dens[[1]]$x, type="l", col="blue", lwd=2, lty=3) points(dens1[[2]]$y*x, dens1[[2]]$x, type="l", col="red", lwd=2) points(-dens1[[1]]$y*x, dens1[[1]]$x, type="l", col="blue", lwd=2) text(mean(c(0, par()$usr[1])), par()$usr[3], "種田梨沙を起用しない", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="blue") text(mean(c(0, par()$usr[2])), par()$usr[3], "種田梨沙を起用する", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="red") legend("topleft", legend=c("Raw", "PS match"), lty=c(3, 1), bty="n", cex=1.5, lwd=3) tab <- table(tail(Y1, length(Y1)/2), head(Y1, length(Y1)/2)) mcnemar.exact(tab) binarysens(m.obj, Gamma=2, GammaInc=0.1) # sensitivity analysis
