こんなツイートを見かけた。
808 のサークルにアンケート調査をして、新刊を落としたサークルと落とさなかったサークルで、制作にとりかかった時期の分布を出している。
ここで、落とさなかったサークルは平均64日前から、落としたサークルは平均48日前から制作に取り組んでいる、と考察しているが、ヒストグラムの分布を見るにニ峰性で、要約統計量として平均を出すのも、中央値や最頻値を出すのもよろしくない。
落としたサークルの平均48日が、ヒストグラムのみぎから3番目の40-59日のビンに入るので、まあ、悪くないとして、落としたことのないサークルの平均64日が、最も頻度の少ない60-79日のビンに収まってしまうことに違和感がバリバリである。
ここで、数値から離れて、コミケで同人誌を制作するというモチベーションを考えてみると、60-70日前というのは、夏冬ともにコミケ当選の時期と被っており、当選前から制作しようと考えるか、当選後に制作しようと考えるか、という仮定をおいてみる。
すると、このヒストグラムは、新刊を落としたことのある/ない、に加えて、コミケ当選（60-70日程度前）の前/後で制作を開始する、が混合した分布と考えられる。ビンの高さは左右非対称のように見えるので、正規分布ではなくガンマ分布を想定した。
ガンマ分布は形パラメータをふたつとり、
推定した結果、
[[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
さて、コミケ当選前後で制作を開始するという分類を入れたのでこの別で見てみる。
以下の行列は制作開始前の平均日数である。行はコミケ当選後に制作を始めるかコミケ当選前から制作を始めるサークルかで、列は新刊を落としたことがあるか、ないか、である。コミケ当選後では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)
