Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-06-12

MikuHatsune2018-06-12

なぜ人は新刊を落とすのか

こんなツイートを見かけた。

C93 で買いそびれていた。

 

808 のサークルにアンケート調査をして、新刊を落としたサークルと落とさなかったサークルで、制作にとりかかった時期の分布を出している。

ここで、落とさなかったサークルは平均64日前から、落としたサークルは平均48日前から制作に取り組んでいる、と考察しているが、ヒストグラムの分布を見るにニ峰性で、要約統計量として平均を出すのも、中央値や最頻値を出すのもよろしくない。

落としたサークルの平均48日が、ヒストグラムのみぎから3番目の40-59日のビンに入るので、まあ、悪くないとして、落としたことのないサークルの平均64日が、最も頻度の少ない60-79日のビンに収まってしまうことに違和感がバリバリである。

f:id:MikuHatsune:20180612162623p:image

 

ここで、数値から離れて、コミケ同人誌を制作するというモチベーションを考えてみると、60-70日前というのは、夏冬ともにコミケ当選の時期と被っており、当選前から制作しようと考えるか、当選後に制作しようと考えるか、という仮定をおいてみる。

すると、このヒストグラムは、新刊を落としたことのある/ない、に加えて、コミケ当選(60-70日程度前)の前/後で制作を開始する、が混合した分布と考えられる。ビンの高さは左右非対称のように見えるので、正規分布ではなくガンマ分布を想定した。

 

ガンマ分布は形パラメータをふたつとり、G(s_1,s_2) と書くと、コミケ当選後に制作を始めるサークルの割合をpとすると、コミケ当選前から制作を始めるサークルの割合は1-p となり、混合分布はpG_1+(1-p)G_2 となる。これをヒストグラムの累積分布から、落としたことのあるサークル、落としたことのないサークル、のパラメータ各2つずつ、とp ゴリ押しで推定する。

 

推定した結果、pコミケ当選後から制作するサークルのガンマ分布パラメータふたつ、コミケ当選前から制作するサークルのガンマ分布パラメータふたつ、の順で並べると以下の通りである。

[[1]]
[1]  0.8531444  4.0315755  0.1158495 95.3276615  0.9989746

[[2]]
[1]  0.7308913  7.7311609  0.1909049 98.0969186  1.0017727

混合ガンマ分布での平均値は以下の通りで、コミケ当選前後で制作を開始するサークルを別にしなければ、12日程度の差がある。

[1] 43.69920 55.95106

さて、コミケ当選前後で制作を開始するという分類を入れたのでこの別で見てみる。

f:id:MikuHatsune:20180612163620p:image

以下の行列は制作開始前の平均日数である。行はコミケ当選後に制作を始めるかコミケ当選前から制作を始めるサークルかで、列は新刊を落としたことがあるか、ないか、である。コミケ当選後では6日程度の差があるが、コミケ当選前では2日である。総製作期間に対する割合を考慮するとあまり差がないと考えられる。

         [,1]     [,2]
[1,] 34.80011 40.49744
[2,] 95.42551 97.92333

声優統計を書いていたときはコミケ当選前からネタ探しはして、実際にとりかかるのはコミケ当選が発表されてからすぐだったのでだいたい60日前くらいだろうか。8回書かせてもらって落としたことはない(ドヤァア

txt <- "
19 71 23
39 169 105
59 117 119
79 8 11
99 57 63
100 22 43
"

d <- as.matrix(read.table(text=txt, row=1))

day <- c(0, as.numeric(rownames(d)))

cols <- c("skyblue", "hotpink")
led <- sprintf("新刊を落とした経験が%s", c("ある", "ない"))
yl <- c(0, 180)
ab <- seq(0, 180, by=20)
barplot(t(d), beside=TRUE, col=NA, border=NA, las=1, ylim=yl, xaxt="n", yaxt="n", xlab="日数")
abline(h=ab, col=grey(0.8))
barplot(t(d), beside=TRUE, col=cols, las=1, ylim=yl, add=TRUE, yaxt="n")
axis(2, at=ab, labels=ab, las=1)
legend("topright", legend=led, col=cols, pch=15, bg="white")

# cumulative density
dens <- sweep(apply(d, 2, cumsum), 2, colSums(d), "/")

# 最適化でゴリ押しする
y <- head(dens[,1], -1)
resid <- function(par){
  yhat <- par[1]*pgamma(day[2:(length(day)-1)], par[2], par[3]) + (1-par[1])*pgamma(day[2:(length(day)-1)], par[4], par[5])
  sum((y-yhat)^2)
}

niter <- 10000
res <- vector("list", niter)
fun <- function() optim(c(runif(1), rpois(1, 40), runif(1, 0, 10), rpois(1, 90), runif(1, 0, 10)), resid, method="BFGS")
for(i in seq(niter)){
  res[[i]] <- try(fun(), silent=TRUE)
}

res <- res[sapply(res, length) > 1] # 推定がうまくいかなかった場合を除く
res[[which.min(sapply(res, "[[", "value"))]]
ps <- res[[which.min(sapply(res, "[[", "value"))]]$par
# ps には新刊を落としたことがある場合と落としたことがない場合をlist でいれる

X <- 0:150
es <- mapply(function(z) z[1]*dgamma(X, z[2], z[3]) + (1-z[1])*dgamma(X, z[4], z[5]), ps)
matplot(es, type="l", col=cols, lty=1, lwd=4, xaxt="n", xlab="日数", ylab="Density")
abline(v=seq(0, 100, by=20), lty=3, col=1)
legend("topright", legend=led, col=cols, pch=15, bg="white")
axis(1, at=seq(0, max(X), by=10), seq(0, max(X), by=10))

# コミケ当選前後をわけない
colSums(es * X)
# コミケ当選前後をわける
mapply(function(z) c(z[2]/z[3], z[4]/z[5]), ps)

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/MikuHatsune/20180612/1528790421