关于贝叶斯结构时间序列模型的问题

机器算法验证 时间序列 贝叶斯 马尔可夫链蒙特卡罗 bsts
2022-03-12 01:03:17

我正在研究在 R 中使用bsts 包的贝叶斯结构时间序列模型的结果的稳定性。以下代码为三种不同数量的 MCMC 绘制（100、1000、 10000）。对于这三种情况，我重复估计 10 次并将 R 方存储在数据框中：

library(bsts)
data(AirPassengers)
y <- log(AirPassengers)
df=data.frame(niter100=NA,niter1000=NA,niter10000=NA)
dfcolumn=0
for (j in c(100,1000,10000)) {
dfcolumn=dfcolumn+1
  for (i in 1:10) {
  ss <- AddLocalLinearTrend(list(), y)
  model_benchmark <- bsts(y,state.specification = ss,niter = j)
  summary=summary(model_benchmark)
  df[i,dfcolumn]=summary$rsquare
  }
}

结果如下：

> df
        niter100 niter1000  niter10000
    1  0.9058217 0.8959122  0.8352333
    2  0.9058217 0.6254595  0.7148984
    3  0.9058217 0.8956490  0.8840317
    4  0.9058217 0.9071929  0.8682971
    5  0.9050454 0.9076566  0.9017717
    6  0.9050454 0.8904109  0.8416038
    7  0.9050454 0.9073501  0.8674943
    8  0.9050454 0.9059262  0.8563360
    9  0.9050454 0.9070879  0.8585177
    10 0.9050454 0.6612644  0.8920700

我的期望是，随着我们增加 MCMC 抽签次数，结果应该会变得更加准确和稳定。然而，上面的测试向我表明，随着 MCMC 绘制次数的增加，结果往往会变得不太稳定，即对于 100 次迭代，R 方的变化很小，而对于 1000 次迭代，它在 0.63 和 0.91 之间变化。这可能是什么原因？有什么策略可以在实际应用中处理这个问题吗？

1个回答

我可以看到您的示例存在几个潜在问题：

  1. 你没有指定一个seed所以bsts会使用系统时钟和连续蒙特卡罗运行之间的串行相关性会弄乱你的统计数据
  2. 您选择的指标rsquare可能不是您认为的那样（请参阅帮助summary.bsts）
  3. 您的模型不太适合数据，因此可能需要大量样本才能收敛

将您的代码修改为地址 1 和 2...

library(bsts)
data(AirPassengers)
y <- log(AirPassengers)

res <- list()

for (j in c(100,1000,10000)) {
  res.inner <- list()
  for (i in 1:10) {
    ss <- AddLocalLinearTrend(list(), y)
    # ss <- AddSeasonal(ss, y, nseasons = 12)
    seed <- floor(i*j*(as.numeric(Sys.time()) %% pi))
    model_benchmark <- bsts(y, state.specification = ss, niter = j, seed=seed)
    x <- summary(model_benchmark)
    res.inner[[i]] <- c(x$residual.sd, x$prediction.sd, x$rsquare, x$relative.gof)
  }
  df.inner <- Reduce(rbind, res.inner)
  colnames(df.inner) <- c("residual.sd", "prediction.sd", "rsquare", "relative.gof")
  res[[j]] <- df.inner
  print(res[[j]])
}

X <- Reduce(rbind, lapply(res, function(x) {if (length(x) > 0) apply(x,2,sd)}))
row.names(X) <- c("100", "1000", "10000")

X

...产生与您相似的结果：

       residual.sd prediction.sd     rsquare relative.gof
100   0.0005273839  2.215716e-05 0.000728682 0.0005939483
1000  0.0163049235  5.001772e-04 0.026663130 0.0133087130
10000 0.0244355253  6.447072e-04 0.042745282 0.0170885684

现在，如果我们通过取消注释这一行来添加季节性术语

# ss <- AddSeasonal(ss, y, nseasons = 12)

我们得到：

       residual.sd prediction.sd      rsquare relative.gof
100   0.0022619056  2.987453e-04 0.0007276338 0.0031640530
1000  0.0018138912  1.890575e-04 0.0005448601 0.0020016078
10000 0.0007862781  2.854752e-05 0.0002245090 0.0002997004

所以看起来罪魁祸首是第 3 位——局部线性趋势并不适合高度季节性的 AirPassenger 数据集。

其它你可能感兴趣的问题
上一篇非线性内生性 下一篇使用 Epanechnikov 内核的 np 包内核密度估计