使用 Reversible-Jump MCMC Bayesian 方法在 R 中更改点包

Question

我正在寻找一种方法（或至少是一个 R 包）来使用可逆跳跃 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 包，但它不适用于变点分析： https://cran.r-project.org/web/packages/rjmcmc/rjmcmc.pdf

任何人都可以在正确的包中指出我或者至少帮助我如何在 R 中执行此操作吗？我将不胜感激任何帮助。

Answer 1

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

之前通过平方根转换计数数据以使其更像高斯分布

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

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 的。

使用 Reversible-Jump MCMC Bayesian 方法在 R 中更改点包

Change point package in R using Reversible-Jump MCMC Bayesian approach

r

bayesian