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)

- MasK

Answer 1

"Value" 部分 文档 描述了 MCMCregressChange 返回的内容，说明模型的对数边际似然存储在属性 logmarglike 中。因此，可以像下面这样访问它：

attr(model1, "logmarglike")

这些属性值在运行代码时也会被报告：

print(BayesFactor(model0, model1, model2))

关于模型的区别，MCMCresidualBreakAnalysis是MCMCregressChange的特殊情况，即X为单变量时。实际上，MCMCregressChange的代码会检查X中的列数，如果为1，则重新格式化输入参数以调用MCMCresidualBreakAnalysis。由于后者没有特定的额外参数，因此了解MCMCregressChange更加通用，而且使用它就足够了。 MCMCresidualBreakAnalysis描述中的一条注释进一步强化了这一点：

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

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

- merv

Answer 2

除了MCMCpack之外，我认为一些专门设计用于变点检测的贝叶斯模型可能会很有用。在R中，有三个可能的包，分别是、和。在拟合模型和处理数据类型方面，和更加灵活。是一种特定于同时贝叶斯时间序列分解（类似于）和变点检测（类似于）的方法；还报告后验对数边际似然，可用于比较关于变点的备选假设（更准确地说，是关于变点数量的备选先验）。

以下是使用和对您的样本数据进行快速结果展示的一些内容。

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的图表：

以下是只指定了一个变点的model1的图表。请注意，在先验中，我们只指定变点的数量为1，Rbeast仍会找出最可能的位置；它估计其随时间发生的概率（即绿色Pr(tcp)曲线），并确定最可能的位置（即垂直虚线）。order_t曲线描述了需要充分拟合趋势所需的多项式曲线的平均阶数；这里它是一条零线，因为我们将其固定为零（即torder.minmax=c(0,0)）。

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

- zhaokg

---

忘了说我是Rbeast的作者。显然，我竭尽全力地满足自己分享软件包并尝试提供潜在帮助的自负心态。 - zhaokg