1.2. bcp贝叶斯检验

InitCluster(20)
lst <- dt_dlply(df_mon, .(site),
  ~ cpt_bcp(.$value, .$date),
  .progress = "text", .parallel = TRUE
) %>% rm_empty()
# lst = map(lst, ~set_names(.x, c("date", "cpt", "prob")))
length(lst)

# map_dbl(lst, nrow) %>% table() %>% print()
# lst_homo = lst[sites_homo] %>% rm_empty() #%>% length()

info = map(lst, function(d) {
  d[which.max(prob), ]
  # data.table(TP = d$date[1], 
  #   mean_a = d$mean[1], mean_b = d$mean[2], 
  #   var_a = d$variance[1], var_b = d$variance[2]
  # )
}) %>% melt_list("site")
# bcp检测不出来
write_fig({
  # par(mfrow = c(1, 3))
  # hist(info[, var_b - var_a])
  # hist(info[, mean_b - mean_a])
  hist(info[, year(date)], main = "突变年份")
  # hist(info[n_miss > 0 & site %in% siteHomoInfo$site, ]$perc_miss*100)
}, 'Figure2_TP_bcp.pdf', 10, 5, show = F)
# cpt检验
# 均值下降,方差变大,且2005附近时间;
code("Figure2_TP_bcp.pdf")


kongdd/RHtestsHelper documentation built on April 18, 2023, 1:57 a.m.