R语言 从ACF相关图中提取置信区间

bzzcjhmw  于 2023-02-10  发布在  其他
关注(0)|答案(3)|浏览(313)

在R中,我们可以运行时间序列的ACF相关图,置信区间带将以浅蓝色绘制。但是当我拉ACF对象的结构时,我找不到这些值。有人知道如何提取置信区间带的值吗?
例如:

List of 6
 $ acf   : num [1:27, 1, 1] 1 0.06453 -0.06354 0.00213 -0.01324 ...
 $ type  : chr "correlation"
 $ n.used: int 501
 $ lag   : num [1:27, 1, 1] 0 1 2 3 4 5 6 7 8 9 ...
 $ series: chr "tser[i:(i + 500)]"
 $ snames: NULL
 - attr(*, "class")= chr "acf"

t5fffqht

t5fffqht1#

我已经看过这个函数了,但我找不到一个简单的方法来提取置信区间。区域是在plot.acf函数中计算的。要查看这个函数,请使用

getS3method("plot", "acf")

在这个函数中,有一个变量clim,这就是你要找的变量。最简单的方法是把plot.acf复制到myplot.acf,但是返回clim的值。

tzdcorbm

tzdcorbm2#

我知道这个问题已经很老了,但是如果有人想要置信区间的值,那就是置信水平的z值除以观测值的平方。在plot.acf函数中,计算如下:

clim0 <- if (with.ci) 
    qnorm((1 + ci)/2)/sqrt(x$n.used)

其中with.ci是指示用户是否想要绘制置信区间的逻辑值,并且ci是期望的置信水平(例如.95、.9等)
编辑:这是置信区间,如果你假设滞后值是白色噪声,如果不是这种情况,有一个修正,你可以应用

clim <- clim0 * sqrt(cumsum(c(1, 2 * x$acf[-1, i, j]^2)))

您可以阅读更多有关here的信息

50pmv0ei

50pmv0ei3#

好的,我们有一个序列X,我们使用内置的stats::acf函数来计算自相关函数的值。

X     <- c(seq(20,10,-1),seq(1,20))
X_ACF <- acf(X) # by default the same as `acf(X, ci.type="white")`

您将得到一个图,其中置信区间为常数值acf(X, ci.type="white")(对于默认白噪声零假设)或非常数值acf(X, ci.type="ma")(对于移动平均假设)。* 有关差异的信息,请参阅plot.acf的文档。*

然而,与直觉相反的是,acf()返回的对象中并没有包含这些图中的置信区间数据,但是,你仍然可以很容易地自己得到这些数据,为了直接回答你的问题,这里有一个函数可以从一个"acf"对象中得到这些置信区间(受@csgillespie建议的启发):

get_clim <- function(x, ci=0.95, ci.type="white"){
  #' Gets confidence limit data from acf object `x`
  if (!ci.type %in% c("white", "ma")) stop('`ci.type` must be "white" or "ma"')
  if (class(x) != "acf") stop('pass in object of class "acf"')
  clim0 <- qnorm((1 + ci)/2) / sqrt(x$n.used)
  if (ci.type == "ma") {
    clim <- clim0 * sqrt(cumsum(c(1, 2 * x$acf[-1]^2))) 
    return(clim[-length(clim)])
  } else {
    return(clim0)
  }
}

像这样使用它
一个一个二个一个一个一个三个一个一个一个一个一个一个四个一个一个一个一个一个五个一个
现在,为了证明这是有效的,并且因为它可能是有用的,这里有一个函数,它使ggplot2图与上面的默认基本R图相对应。

library(ggplot2)
theme_set(theme_minimal())
ggplot_acf <- function(
    x,
    ci=0.95, ci.type="white", ci.col = "blue"){
  #' Replicates plot.acf() but using ggplot by default instead of base R plot
  #' `x` must be an object of class "acf" such as that outputted by `acf()`
  #' `ci.type` must be "white" or "ma"
  if (!ci.type %in% c("white", "ma")) stop('`ci.type` must be "white" or "ma"')
  if (class(x) != "acf") stop('pass in object of class "acf"')
  with.ci <- ci > 0 && x$type != "covariance"
  with.ci.ma <- with.ci && ci.type == "ma" && x$type == "correlation"
  if(with.ci.ma && x$lag[1L, 1L, 1L] != 0L) {
    warning("can use ci.type=\"ma\" only if first lag is 0")
    with.ci.ma <- FALSE
  }
  clim <- get_clim(x, ci=ci, ci.type=ci.type)
  df <- data.frame(lag = x$lag, acf=x$acf)
  p <- ggplot(df, aes(x=lag)) +
    geom_linerange(aes(ymax=acf, ymin=0)) +
    labs(y="ACF", x="Lag")
  if (with.ci) {
    if (ci.type == "white") {
      p <- p + 
        geom_hline(yintercept = 0-clim, lty = 2, col = ci.col) +
        geom_hline(yintercept = 0+clim, lty = 2, col = ci.col)
    } else if (with.ci.ma && ci.type == "ma") { # ci.type="ma" not allowed for pacf
      dfclim <- df[-1,]
      dfclim$clim <- clim
      p <- p + 
        geom_line(data = dfclim, aes(y = 0-clim), lty = 2, col = ci.col) +
        geom_line(data = dfclim, aes(y = 0+clim), lty = 2, col = ci.col)
    }
  }
  return(p)
}

为了检查这是否正常工作,我们将生成的ggplot对象绘制在plot.acf生成的相应基本R图旁边。

library(patchwork)
p11 <- ggplot_acf(X_ACF, ci.type="white") + labs(subtitle="ggplot version")
p12 <- wrap_elements(panel=~plot(X_ACF, ci.type="white"))  + labs(subtitle="base R version")

old_par <- par(mar = c(0,0,0,0), bg = NA)
(p11+p12) 
par(old_par)

p21 <- ggplot_acf(X_ACF, ci.type="ma") + labs(subtitle="ggplot version")
p22 <- wrap_elements(panel=~plot(X_ACF, ci.type="ma")) + labs(subtitle="base R version")

old_par <- par(mar = c(0,0,0,0), bg = NA)
(p21+p22) 
par(old_par)

相关问题