HUNTER x HUNTERの連載期間を推定 - 逆にありかもしれない

みなさん、ハンターハンターは好きですよね。

緻密な心理戦、登場人物の成長と戦略の変化など、冨樫先生の才能をじかに感じられる喜びがあります。

しかし、リアルタイム勢にとって最もスリリングなのは、「いつまで連載するのか」という点です。

さて、ご存じの通り、海外有志がHxHの連載を記録しています。

hiatus-hiatus.github.io

このデータをもとに、ベイズモデリングでハンターハンターの連載期間の推定をします。

参考文献は

ajhjhaf.hatenablog.com

Bayesian Inference With Stan ~062~ | |

の二つのサイトです。

AR過程で推定する戦略もあるようですが、ワイブル分布のほうがstanが単純だったので、ワイブル分布を採用します。

ある微小時間のイベント回数がワイブル分布に従っているというモデルだと思います。

f(t) ~ weibull(t, shape, scale)で、

f(t)の累積密度関数をF(t)とすると、

F(t)は時点tまでの死亡率、1-F(t)が生存関数になると思います。

さて、データをまとめると

week	arrest
12	1
12	1
8	1
1	1
4	1
2	1
11	1
1	1
1	1
3	1
10	1
3	1
2	1
3	1
9	1
5	1
2	1
1	1
5	1
2	1
5	1
2	1
3	1
1	1
3	1
3	1
3	1
4	1
1	1
13	1
4	1
2	1
4	1
3	1
3	1
2	1
3	1
4	1
3	1
3	1
3	1
3	1
3	1
2	1
1	1
1	1
1	1
2	1
1	1
1	1
2	1
2	1
1	1
2	1
3	1
3	1
2	1
3	1
3	1
3	1
3	1
3	1
2	1
2	1
1	1
4	1
10	1
10	1
10	1
21	1
32	1
2	1
5	1
3	1
1	1
11	1
10	1
10	1
3	0

となります。

arrest = 1　なら連載終了を観測、arrest = 0なら観測打ち切りです。

これをHxHdata1.csvとして保存しました。

カプラン・マイヤー曲線の作成

まずベイズではなく、通常の生存時間曲線を描きます。

HxHdata1.csvにstatus列を追加して、status = 1 なら連載中、status = 0なら休載中というデータをつくり、休載期間も合わせてデータし、HxHsf.csvとして保存しました。

survivalのパッケージをインストール、展開し、

Hsf <- survfit(Surv(week, arrest) ~ status, data=HxHsf)

と入力します。

f:id:oooaki47ooo:20181028231232p:plain

プロットすると、このようになります。縦軸が割合、横軸がweekです。

破線が連載期間、実線が休載期間の生存時間曲線です。

二度、長期の連載期間が存在するために、生存期間曲線がいびつです。

これを見る限り、皆さんの実感どおり12週間以上連載する確率は低いように思われます。

続いて、stanでモデリングを行います。

比例ハザードモデルにしないので、stanコードは比較的単純です。

以下

data {
int N;
int week[N];
real arrest[N];
}

parameters {
real shape;
real scale;
}

model{
for(n in 1:N){
if (arrest[n] == 0){
target += weibull_lccdf(week[n]| shape, scale);
}else{
target += weibull_lpdf(week[n]| shape, scale);
}
}
}

generated quantities{
real pred_Y[N];

for(n in 1:N){
pred_Y[n] = (1-weibull_cdf(week[n], shape, scale));
}
}

以上

参考文献のコードはほぼ同じでしが、generated quantitiesのセクションが異なっています。

私の目的は生存時間曲線の作成であるため、「Rで学ぶベイズ統計学」さんの方法に従いました。

加えて、「Easy to Type」さんはpystanのユーザーのようである点もネックでした。

ワイブル分布の詳細はアヒル本や上記ブログを参考にしていただければ助かります。

指数分布やレイリー分布の一般化だと思います。

続いて、Rのコード。

library(rstan)
HxHdata1 <- read_delim("~/Rstan/HxH/HxHdata1.csv",
";", escape_double = FALSE, trim_ws = TRUE)
N <- length(HxHdata1$week)
d <- HxHdata1

data <- list(week = d$week, arrest = d$arrest, N = N)
fit <- stan(file='~/Rstan/HxH/HxHmodel1.stan',
data = data, seed = 1234)

save.image('~/Rstan/HxH/HxH_result1.RData')

write.table(data.frame(summary(fit)$summary),
file='~/Rstan/HxH/model1_fit_summary.txt', sep='\t', quote=FALSE, col.names=NA)

library("tidyverse")

load('~/Rstan/HxH/HxH_result1.RData')

PRED1 <-rstan::extract(fit)$pred_Y %>%
data.frame() %>%
apply(., 2, quantile, prob = c(0.025, 0.10, 0.5, 0.90, 0.975)) %>%
t() %>%
`colnames<-`(c("cred2.5", "cred10", "cred50", "cred90", "cred97.5")) %>%
cbind(d)
PRED2 <- do.call(rbind, PRED1)

G1 <- ggplot(data.frame(PRED2))
G2 <- G1 + geom_line(data = PRED1, aes(week, y = cred50), colour = "#4169e1")
G3 <- G2 + geom_ribbon(data = PRED1, aes(week, ymin = cred2.5, ymax = cred97.5), alpha = 0.1, fill = "#4169e1")
G4 <- G3 + geom_ribbon(data = PRED1, aes(week, ymin = cred10, ymax = cred90), alpha = 0.3, fill = "#4169e1")
G5 <- G4 + theme_classic()
G6 <- G5 + theme(axis.text = element_text(size = 15), axis.title = element_text(size = 20))

以上

初心者なので、パイプ関数をうまく使えないため、参考のコードから変更できていません。

少し余裕ができたら勉強します。

結果です。

	mean	se_mean	sd	X2.5.	X25.	X50.	X75.	X97.5.	n_eff	Rhat
shape	1.1439864914	0.0018472169	0.0906771403	0.9668414276	1.083567748	1.1425098503	1.2064211912	1.3174975306	2409.682732168	0.9994322193
scale	4.8718696817	0.0116456112	0.5322446165	3.9225171134	4.5060828684	4.8579300028	5.2085245521	5.9620653922	2088.8053980174	1.0002681537
lp__	-195.0351128272	0.0251603167	1.0325026309	-197.8198902233	-195.4283494359	-194.7405790845	-194.2976895707	-194.0150298631	1684.0311806587	1.0024096459