使用 Reversible-Jump MCMC 贝叶斯方法在 R 中更改点包答案

【问题标题】：Change point package in R using Reversible-Jump MCMC Bayesian approach使用 Reversible-Jump MCMC 贝叶斯方法在 R 中更改点包
【发布时间】：2022-01-16 00:13:07
【问题描述】：

我正在寻找一种方法（或至少是一个 R 包）来使用 Reversible-jump MCMC 方法执行贝叶斯变点分析。

我将应用它来检测台风时间序列中的变化点。

这是我的参考论文：https://journals.ametsoc.org/doi/pdf/10.1175/JCLI-D-13-00744.1

我还想为每个变化点绘制后验概率质量函数。

这里是示例数据：

structure(list(V1 = c(7L, 6L, 4L, 4L, 4L, 2L, 5L, 4L, 4L, 4L, 
6L, 7L, 8L, 6L, 10L, 7L, 9L, 5L, 1L, 4L, 5L, 5L, 2L, 2L, 5L, 
1L, 2L, 4L, 0L, 3L, 6L, 3L, 6L, 1L, 5L, 3L, 4L, 0L, 2L, 4L)), class = 
"data.frame", row.names = c(NA, 
-40L))

我找到了这个 R-package，但它没有应用于 changepoint 分析： https://cran.r-project.org/web/packages/rjmcmc/rjmcmc.pdf

谁能指出我正确的包或至少帮助我如何在 R 中做到这一点？如有任何帮助，我将不胜感激。

【问题讨论】：

我可能误解了您的问题，但看起来您只需要使用hist 创建一个直方图，只要您已经拥有来自 mcmc 的数据。否则，this package 可能会有所帮助。
@Trusky--我检查了那个，但没有可逆跳转 MCMC。

标签： r bayesian

【解决方案1】：

以防万一这个老问题仍然需要额外的想法，R 中的大多数贝叶斯变点或断点检测包都是通过 Gibbs 采样而不是 Reversible-Jump MCMC 采样实现的。一个例外是我实现的一个包Rbeast (https://cran.r-project.org/web/packages/Rbeast/index.html)，它使用混合可逆跳跃MCMC 采样器来估计变化点概率。需要注意的是，Rbeast 是为类高斯数据而非泊松数据（即参考论文中的计数数据）制定的。在任何情况下，您仍然可以在计数数据上对其进行测试，例如，通过对数转换。在 SouthWood 的书Ecological Methods（第五版）中，第 475 页给出了一个示例，在应用Rbeast 之前，通过平方根对计数数据进行转换，使其更像高斯。

以下是使用您提供的示例数据的一些示例代码和快速结果。请注意，Rbeast 同时进行时间序列分解（如果存在周期性分量）和变化点检测，这使其不同于仅分解时间序列或仅检测中断的 stl 或 changepoint 函数。因为您的样本数据没有季节性/周期性成分，所以在下面的示例中，season='none' 是在 beast 函数中指定的：

Y= c(7L, 6L, 4L, 4L, 4L, 2L, 5L, 4L, 4L, 4L, 6L, 7L, 8L, 6L, 10L, 7L, 9L, 5L, 1L, 
     4L, 5L, 5L, 2L, 2L, 5L, 1L , 2L, 4L, 0L, 3L, 6L, 3L, 6L, 1L, 5L,3L, 4L, 0L, 2L, 4L)
        
library(Rbeast)

# Rbeast is also a tool for decomposing a time series into seasonal and trend components, 
# but here your data is trend-only, so season='none' is used.
out=beast(Y,season='none')
plot(out)

垂直虚线指出最可能的变化点位置。绿色的 Pr(tcp) 图描绘了在单个时间点发生变化点的概率，大致对应于您参考论文的图 d 和 e，但没有区分变化点的顺序（例如，第一个变化点，第二个变化点，... )。

您要查找的后验概率质量函数存储在输出out$trend$ncpPr 中。这是它的条形图。

ncpPr=out$trend$ncpPr[,1]
ncp  = (1:length(ncpPr)) -1
barplot( ncpPr ~ ncp, data=data.frame(ncp, ncpPr), xlab='num of changepoints', ylab='posterior prob')

用于此类分析的另一个出色软件包是 bfast，它不是基于 MCMC 的。

【讨论】：