UPDATE: 2022-12-19 20:42:42
ブートストラップ信頼区間(パーセンタイル法)を書くことがあったので、その備忘録。今回の例では、ワイブル分布 (Weibull distribution) を例にする。この分布は、サバイバル分析などで出てくるように、人の死亡確率をモデリングする際に役立つ。人の死亡確率というと角が立ちそうだが、ECサイトの会員の離脱率なんかをイメージすると良いかもしれない。つまり、観測している期間において死亡確率は一定というよりも、死亡確率は時間と共に高くなり、変化するという仮定を起きたい時に便利。この反対もしかり。その事象が発生するまでの時間を確率変数と考えると、その確率変数が従う分布はワイブル分布となる。
ブートストラップサンプリングとは、得られたサンプルデータから、リサンプリングを行うこと。そこから、リサンプルデータから統計量を求め、推測を行う方法のことをブートストラップ〇〇とか言う。ブートストラップ信頼区間(パーセンタイル法)ではリサンプリングされた値を使ってパーセンタイルを計算すると、その値が取りうる信頼区間が得られる。これを各値で計算していくことで得られる。
{fitdistrplus}
パッケージを使っている。{fitdistrplus}
パッケージの詳細はこちら。
library(fitdistrplus)
library(tidyverse)
set.seed(123)
<- 100
n <- rweibull(n = n, shape = 2, scale = 10)
df_prob <- seq(0.1, 30, len = n)
x
<- 100
loops <- map_dfc(1:loops, function(i) {
booted_df <- sample(df_prob, size = length(df_prob), replace = TRUE)
xi <- suppressWarnings(fitdist(xi, distr = "weibull",method = "mle"))
MLE_est dweibull(x, shape = MLE_est$estimate["shape"], scale = MLE_est$estimate["scale"])
})
<- booted_df %>%
dat1 bind_cols(x %>% as_tibble() %>% rename(x = value)) %>%
gather(key = boot_ind, val = y, -x)
<- booted_df %>%
dat2 array_tree(., 1) %>%
map_dfc(., quantile, c(0.025, 0.5, 0.975)) %>%
t() %>% as_tibble() %>%
bind_cols(x %>% as_tibble()) %>%
setNames(c("2.5%", "50%", "97.5%", "x")) %>%
gather(key = quantile, val = y, -x)
ggplot() +
geom_line(data = dat1, aes(x, y, group = boot_ind), col = "#006E4F", alpha = 0.2, size = 0.25) +
geom_line(data = dat2, aes(x, y, group = quantile), col = "#006E4F", size = 0.5, linetype = "longdash") +
theme_classic() + xlab("x") + ylab("Probability density") +
ggtitle("Weibull Probability Distribution")