如何从Fable包预测中提取模拟数据？答案

【问题标题】：How to extract simulation data from Fable package forecast?如何从Fable包预测中提取模拟数据？
【发布时间】：2023-01-26 17:22:01
【问题描述】：

在运行底部发布的 R 代码时，我使用寓言包中的 forecast() 函数并运行 1000 个模拟样本路径，以 80% 和 95% 的置信度得出接下来 10 个周期的预测，如下所示：

生成的寓言对象在 R Studio 控制台中如下所示：

我想从上面的 Fable 对象访问模拟路径，这样我就可以绘制预测的分布，例如在第 20 期，如下面示例中的概念所示。任何想法如何做到这一点？

代码：

library(feasts)
library(fable)
library(fabletools)
library(ggplot2)
library(tsibble)

tmp <- data.frame(
  Month = c(1,2,3,4,5,6,7,8,9,10),
  StateX = c(1527,1297,933,832,701,488,424,353,302,280)
  ) %>%
  as_tsibble(index = Month)

fit <-  tmp %>% model(NAIVE(StateX))

fc <- fit %>% forecast(h = 10, bootstrap = TRUE, times = 1000)

autoplot(fc, tmp) +
  labs(title="Transitions to Dead State X", y="Units" )

【问题讨论】：

您可以使用分布上的 parameters() 函数从分布中获取参数（在本例中为样本分布中的样本）。试试parameters(fc$StateX)。
我试过 parameters() 但我收到错误消息“参数错误（fc$StateX）：找不到函数“参数”。parameters() 是包的一部分吗？

标签： r forecasting fable-r

【解决方案1】：

请参阅此链接上的帖子，它更完全地解决了这个问题：https://community.rstudio.com/t/how-to-derive-a-forecast-density-for-a-specified-point-in-time-with-the-r-fable-package-forecast-function/158329。包 ggdist 提供了必要的功能。适配OP代码示例，见下方解决代码：

library(dplyr)
library(fable)
library(ggdist)
library(ggplot2)
library(tsibble)

tmp <- data.frame(
  Month = c(1,2,3,4,5,6,7,8,9,10),
  StateX = c(1527,1297,933,832,701,488,424,353,302,280)
  ) %>%
  as_tsibble(index = Month)

# Fit model
fit <- tmp |>
  model(naive = NAIVE(StateX))

# Generate forecasts
fc <- bind_rows(
  forecast(fit, h = 10) |> mutate(.model = "Theoretical"),
  forecast(fit, h = 10, bootstrap = TRUE, times = 1000) |> mutate(.model = "Bootstrapped")
)

# Density + histogram plot (alpha below sets opacity)
ggplot() +
  stat_dist_slab(aes(dist = StateX, y = 0, fill = "Theoretical"),
                 data = fc |> filter(.model == "Theoretical"),
                 colour = NA, alpha = 0.3
  ) +
  geom_histogram(aes(x = x, y = after_stat(ndensity), fill = "Bootstrapped"),
                 data = tibble(x = fc |>
                                 filter(.model == "Bootstrapped") |>
                                 pull(StateX) |>
                                 unlist()),
                 colour = NA, bins = 50, alpha = 0.3
  ) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  labs(x = "Units reaching dead state X (at h=10)", y = "Forecast density") +
  theme(legend.position = "none")

【讨论】：