HUNTER x HUNTERの連載期間を推定
みなさん、ハンターハンターは好きですよね。
緻密な心理戦、登場人物の成長と戦略の変化など、冨樫先生の才能をじかに感じられる喜びがあります。
しかし、リアルタイム勢にとって最もスリリングなのは、「いつまで連載するのか」という点です。
さて、ご存じの通り、海外有志がHxHの連載を記録しています。
このデータをもとに、ベイズモデリングでハンターハンターの連載期間の推定をします。
参考文献は
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)
と入力します。
プロットすると、このようになります。縦軸が割合、横軸が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 |
ちゃんと収束したようです。
続いて、生存曲線。
カプランマイヤーとよく似たモデルができているようです。
実際に必要なのはshape 1.144、 scale 4.872です。
これをワイブル分布に用いると、来週や再来週HxHが掲載されるかが予測できるはずです。
眠いので、おやすみなさい。
beast regard
アヒル本4章問題のこり
前回に続き、松浦健太郎先生著の「StanとRでベイズ統計モデリング (Wonderful R)」（以下アヒル本）の練習問題をやります。
4章（５）が残っています。
Y1は、平均0、標準偏差5の正規分布に従う標本数30のベクトル、
Y2は、平均1、標準偏差4の正規分布に従う標本数30のベクトルのt検定をせよだそうです。
データ設定
以下コード
set.seed(123)
N1<- 30
N2<- 20
Y1 <-rnorm(n=N1, mean = 0, sd =5)
Y2 <- rnorm (n=N2, mean = 1, sd = 4)
以上
以上をexercise4-data.Rの名前で保存します。
stanコードを書きます。
全く同じでは、流石につまらないので、ベクトル化してみます。
以下コード
data {
int N1;
int N2;
vector[N1] Y1;
vector[N2] Y2;
}
parameters {
real mu1;
real mu2;
real<lower = 0> sigma1;
real<lower = 0> sigma2;
}
model {
Y1 ~ normal(mu1, sigma1);
Y2 ~ normal(mu2, sigma2);
}
以上をex5v.stanと名付けて保存します。
データセクションの
real Y1[N1];
が実数列の書き方、
ベクトル化すると
vector[N1] Y1;
になります。処理が長い場合にベクトル化のほうがやや高速化するとか。
今回は単純なのであまり変わりません。
デメリットは可読性が下がり、バグが入りやすくなるそうです。
なので最初にforでモデルを書いて動くことを確認したのちに、
ベクトル化してchainを増やしたり、サンプリングを増やしたりするそうです。
つづいてRのコード
library(rstan)
source('~/Rstan/chap4/exercise4-data.R')
data <- list(N1 = N1, N2 = N2, Y1 = Y1, Y2 =Y2)
fit <- stan(file='~/Rstan/chap4/ex5v.stan', data = data, seed = 1234)
save.image('result-ex5v.RData')
ms5 <- rstan::extract(fit)
prob5 <- mean(ms5$mu1 < ms5$mu2)
prob5
以上をex5v.Rと名付けて保存しました
prob5は 0.9335でした
| mean | se_mean | sd | X2.5. | X25. | X50. | X75. | X97.5. | n_eff | Rhat | |
| mu1 | -0.24 | 0.02 | 0.98 | -2.21 | -0.86 | -0.23 | 0.42 | 1.67 | 3297.23 | 1.00 |
| mu2 | 1.59 | 0.01 | 0.81 | -0.01 | 1.07 | 1.61 | 2.11 | 3.18 | 3662.62 | 1.00 |
| sigma1 | 5.15 | 0.01 | 0.71 | 3.93 | 4.65 | 5.08 | 5.57 | 6.75 | 3217.31 | 1.00 |
| sigma2 | 3.62 | 0.01 | 0.63 | 2.62 | 3.17 | 3.55 | 3.98 | 5.07 | 3279.02 | 1.00 |
| lp__ | -95.37 | 0.04 | 1.51 | -99.20 | -96.11 | -95.02 | -94.25 | -93.52 | 1765.64 | 1.00 |
結果も収束していますが、問題（４）と同じくmu2が大きいです。
収束をしらべます。
以下Rコード
library(ggmcmc)
load('result-ex5v.RData')
write.table(data.frame(summary(fit)$summary),
file='~/Rstan/chap4/fit-summary5v.txt', sep='\t', quote=FALSE, col.names=NA)
p <- ggs_traceplot(ggs(fit, inc_warmup=TRUE, stan_include_auxiliar=TRUE))
p <- p + theme_bw(base_size=18)
p <- p + labs(color='Chain')
ggsave(p, file='20181026XYZ.png', dpi=300, w=7, h=10)
以上
収束結果を図示します。
いい感じに自己相関なくサンプリングできている気がします。
best regard
Rstan単回帰練習
初心者です。
伝統的統計は記述統計の発展であり、（手続きを理解できている場合は）データの直感的理解を促進します。
一方、機械学習は的確な予測を与えますが、カーネルや多層化など直感を超える処理を含む印象があります。
さて、ベイズモデルはパラメーターが分布で得られることから、モデリングの「理解」と「予測」を両立することが特徴的だと考えています。
さて、久保拓弥「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・ＭＣＭＣ」と豊田秀樹「基礎からのベイズ統計」を読んだので
松浦健太郎「StanとRでベイズ統計モデリング」（通称アヒル本）を学んでいます。
ゴールドスタンダードの順番です。
書評は多くの方々がされており、自身の実力も不足していることから割愛します。
右も左もわからないので、GitHubの写経を行い、一行ずつ理解していくつもりです。
アヒル本の第4章から行います。
GitHubのアドレス
（１）データ設定
以下コード
set.seed(123)
N1<- 30
N2<- 20
Y1 <-rnorm(n=N1, mean = 0, sd =5)
Y2 <- rnorm (n=N2, mean = 1, sd = 4)
以上
以上をexercise4-data.Rの名前で保存します。
Y1は、平均0、標準偏差5の正規分布に従う標本数30のベクトル、
Y2は、平均1、標準偏差4の正規分布に従う標本数30のベクトルと設定できました。
以下作図コード
source('~/Rstan/chap4/excercise4-data.R')
library(ggplot2)
d1 <- data.frame(group = 1, Y = Y1)
d2 <- data.frame(group = 2, Y = Y2)
d <- rbind(d1, d2)
d$group <- as.factor(d$group)
p <- ggplot(data=d, aes(x = group, y = Y, group = group, col = group))
p <- p + geom_boxplot(outlier.size = 0)
p <- p + geom_point(position=position_jitter(w = 0.4, h = 0), size = 2)
ggsave('file = fig-ex1.png', plot = p, dpi =300, w = 4, h = 3)
以上
初心者なので今まではhist関数で大まかにしか作図してませんでしたが、ggplot2パッケージを進められています。
一つ目のpをプロット
二つ目のpをプロット
三つ目のpをプロット
pというオブジェクトに少しずつ情報を追加していくコンセプトのパッケージのようです。
（２）stanでStudentのt検定
まず、コンソールにlibrary(rstan)を打ち込みます。
このパッケージを読み込まないと、Rstudioがstanの入力を補助してくれませんでした。
以下stanコード
data {
int N1;
int N2;
real Y1[N1];
real Y2[N2];
}
parameters {
real mu1;
real mu2;
real<lower = 0> sigma;
}
model {
for (i in 1:N1)
Y1[i] ~ normal(mu1, sigma);
for (i in 1:N2)
Y2[i] ~ normal(mu2, sigma);
}
以上
上記をex3.stanで保存する。
Rではなくstanなので、コードの最後の行が改行だけであることが必要でした。
また、初心者なので;を書き忘れてrstanに怒られます。
HatenaBlogにコピペするとtabが反映されないので、読みにくいです。
以下R側で実行するコード
library(rstan)
source('~/Rstan/chap4/exercise4-data.R')
data <- list(N1 = N1, N2 = N2, Y1 = Y1, Y2 =Y2)
fit <- stan(file='~/Rstan/chap4/ex3.stan', data = data, seed = 1234)
save.image('result-ex3.RData')
以上
Rhatは収束しているように思えました。
Inference for Stan model: ex3.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu1 -0.23 0.01 0.83 -1.89 -0.79 -0.22 0.34 1.35 3610 1
mu2 1.61 0.02 1.00 -0.33 0.91 1.61 2.28 3.59 3391 1
sigma 4.49 0.01 0.47 3.68 4.16 4.44 4.77 5.53 3446 1
lp__ -97.76 0.03 1.30 -101.21 -98.35 -97.40 -96.83 -96.33 1666 1
どう考えてもExcelに張り付けて表にして張るべきです。
次回に活かします。
（４）mu1 < mu2の確率を求めよ
以下Rコード
library(rstan)
load('result-ex3.RData')
ms <- rstan::extract(fit)
prob <- mean(ms$mu1 < ms$mu2)
以上
GitHub通りではなく、extract関数は多様なパッケージで実装されているんので、パッケージを指定しています。
(豊田本(たしか)の方針に従っています。)
とりあえず、単回帰の練習問題はできました。
(5)はparametersの部分で
real<lower=0> sigama;
を
real<lower=0> sigama1;
real<lower=0> sigama2;
に変え、modelも変更するとできるため略します。
こちらは一応ggmcmcのパッケージでTrace Plotを作図しました。
写経すると、「できない」and「おろそかにしている」部分が明示されるので、身が引き締まります。
best regard
Rstan install
2018 10/21時点
Rバージョン3.5.1
windows10
Rtools35.exeを使用
参考文献
Rtoolsインストール方法
https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started-(Japanese)
Rtoolsのpath編集
https://sudori.info/stat/stat_rcpp_01.html#002win
GCCのPATHの場所の参考
https://notchained.hatenablog.com/entry/2016/02/27/211818
手順1
Rをインストール
Rstudioをインストール
Rtoolsをインストール
ここまでは調べて言うとおりにすれば容易でした。
RtoolsのPath編集の参考ブログのとおり
SystemPropertiesAdvanced
からシステム環境設定を開き、GCCのパスを書きかえを行います。
現在、GCCの保存先が新しくなっているので、
GCCのPATHのPATHの参考のブログをもとにC:\Rtools\mingw_64\binにGCCがあることを確認し、PATHを設定します。
おそらく、これでRstanのインストールの準備が完了します。
>install.paclages("rstam")
>library(rstan)
でRの環境下でstanを使えるようになります。
追記
使えていません。
.shlib_internal(commandArgs(TRUE)) でエラー:
C++14 standard requested but CXX14 is not defined
PATHを間違っているのかもしれません。
