R中的MCMC变点模型

Question

我想运行MCMC线性高斯多变点模型来检测连续值的时间序列向量的变点。

在这样做的时候，我正在考虑使用MCMCregressChange函数，但我在这里有几个问题：

(1)如何获得这些模型的对数边际似然？

(2) MCMCregressChange函数和MCMCresidualBreakAnalysis函数有什么区别？

R脚本如下所示。如果你能帮我解决这个问题，我将非常高兴。

library(MCMCpack)
set.seed(1234)
n <- 100
x1 <- runif(n, min = 0, max  = 1)
x2 <- runif(n, min = 1, max  = 2)
X <- c(x1,x2)

B0 <- 0.1
sigma.mu=sd(X)
sigma.var=var(X)

model0 <- MCMCregressChange(X ~ 1, m=0, b0=mean(X),   mcmc=100, burnin=100, verbose = 1000,
                                   sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
model1 <- MCMCregressChange(X ~ 1, m=1, b0=mean(X),   mcmc=100, burnin=100, verbose = 1000,
                                   sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
model2 <- MCMCregressChange(X ~ 1, m=2, b0=mean(X),   mcmc=100, burnin=100, verbose = 1000,
                                   sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")

print(BayesFactor(model0, model1, model2))

plotState(model0)
plotChangepoint(model0)

plotState(model1)
plotChangepoint(model1)

plotState(model2)
plotChangepoint(model2)

Answer 1

the documentation的"Value“小节描述了MCMCregressChange返回的内容，指出模型的对数边际似然存储在属性logmarglike中。因此，可以像这样访问它

attr(model1, "logmarglike")

在代码中运行行时，也会报告以下属性值：

print(BayesFactor(model0, model1, model2))

至于模型之间的差异，MCMCresidualBreakAnalysis是MCMCregressChange的一个特例，即当X是单变量时。实际上，MCMCregressChange的代码检查X中的列数是否为1，如果是，则将输入参数重新格式化为对MCMCresidualBreakAnalysis的调用。由于也没有专门针对后者的额外参数，因此了解MCMCregressChange更加通用，所有人都应该使用。

在MCMCresidualBreakAnalysis描述中有一条说明加强了这一点：

“此代码主要是为testpanelSubjectBreak__的内部使用而编写的。”

也就是说，虽然它是一个导出函数，但它主要是一个源于特定用例的便利函数。

Answer 2

除了MCMCpack之外，我认为一些专门为变化点检测设计的贝叶斯模型可能会很有用。在R中，有三个可能的包：bcp、mcp和Rbeast。bcp和mcp在模型拟合和处理的数据类型方面更加通用。Rbeast是一种专门用于同时进行贝叶斯时间序列分解(类似于stl)和变点检测(类似于changepoint)的方法；Rbeast还报告可用于比较变点替代假设的后验对数边际似然(更准确地说，用于比较变点数字的替代先验)。

下面是对使用bcp和Rbeast的示例数据的一些快速结果。

set.seed(1234)
n  = 100
x1 = runif(n, min = 0, max  = 1)
x2 = runif(n, min = 1, max  = 2)
X  = c(x1,x2)

library(bcp)
fit=bcp(X)
plot(X)

​

平均而言，bcp精确定位大约2个变化点；最佳位置由概率曲线中的峰值指示。它也可以从sum(fit$posterior.prob,na.rm = TRUE)获得。我认为bcp在这里拟合了分段常数模型；上面绘制的平均曲线是许多MCMC采样的分段常数模型的平均值，这给出了检测到的变化点周围的不规则性。

作为一种时间序列分解模型，Rbeast以“Y=季节/周期(如果存在)+趋势+误差”的形式对时间序列进行拟合。季节分量和趋势分量分别被建模为分段调和曲线和分段线性(多项式)曲线。假设样本数据中没有周期成分，则在以下代码中使用season='none'来拟合仅趋势模型。此外，作为变化点模型，Rbeast允许用户指定可能的变化点数量的范围；如果允许的最小和最大变化点数量相同，并且Rbeast会将变化点的数量固定为常量；例如，tcp.minmax=c(0,0)指定趋势没有变化点。对于每个分段趋势/分段，可以将多项式阶数设置为允许的最小和最大阶数范围。下面，我们将最小和最大阶数固定为零，以便将每个线段拟合为一条恒定线(即torder.minmax=c(0,0))。

library(Rbeast)

model0 = beast(X, season='none', tcp.minmax = c(0,0), torder.minmax = c(0,0) ) # no changepoint
model1 = beast(X, season='none', tcp.minmax = c(1,1), torder.minmax = c(0,0) ) # 1 changepoint
model2 = beast(X, season='none', tcp.minmax = c(2,2), torder.minmax = c(0,0) ) # 2 changepoints

plot(model0)
plot(model1)
plot(model2)

# These are the posterior log marginal likelihoods; the numbers will vary slightly
# across runs due to the MCMC nature.
model0$marg_lik    #   -460.6778
model1$marg_lik    #   -313.9160 (the most likely)
model2$marg_lik)   #   -315.8801 

下面是假设没有变化点的model0图：

​

下面是仅指定了1个变更点的model1的曲线图。请注意，在前面的代码中，我们只是将变化点的数量指定为1，Rbeast仍然会找到最可能的位置；它会估计它随着时间的推移出现的概率(即，绿色Pr(tcp)曲线)加上识别最可能的位置(即，垂直虚线)。order_t曲线描绘了充分拟合趋势所需的多项式曲线的平均阶数；这里它是一条零线，因为我们将其固定为零(即torder.minmax=c(0,0))。

​

​

下面是假设有两个变化点的model2图。估计的变化点概率与bcp结果大致相同。在实践中，运行Rbeast最合理的方法不是通过将变化点的数量固定为已知常量来指定强先验，而是指定一个较大的范围，并让模型确定数量和位置。

​

​