R 中 arima() 函数的计算复杂度是多少?
What is the computational complexity of arima() function in R?
如果我的时间序列有 n
个成员并且我想拟合 ARIMA(1,0,1)
模型,那么大 O
表示法的复杂度是多少?
在下面的例子中我需要知道第二行代码的复杂度:
series <- arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=1000)
result <- arima(series, order=c(1,0,1))
谢谢。
它是 O(n)
复杂度。一个故事吗?见下文。
正如我在评论中所说,我们可以通过回归模型来衡量它。作为玩具演示,请考虑以下数据收集和回归过程。
我们首先定义一个函数来测量ARMA(1,1)
模型(或ARIMA(1,0,1)
)的模型拟合时间。 (请注意,我在这里使用基本的 system.time()
。您可以考虑使用包 microbenchmark
中的 microbenchmark()
。但在下文中,我将使用相当大的 n
, 以降低时间测量的灵敏度。)
t_arma <- function(N) {
series <- arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=N)
as.numeric((system.time(arima(series, order=c(1,0,1)))))[3]
}
我们需要收集,比如 100 条数据。我们尝试 100 个越来越大 n
,并测量模型拟合时间 t
。
k <- 100; t <- numeric(k)
n <- seq(20000, by = 1000, length = k) ## start from n = 20000 (big enough)
for (i in 1:k) t[i] <- t_arma(n[i])
现在,如果我们假设复杂度是:a * (n ^ b)
,或O(n ^ b)
,我们可以通过回归模型估计a
,b
:
log(t) ~ log(a) + b * log(n)
我们对斜坡估算特别感兴趣:b
。
所以让我们调用lm()
fit <- lm(log(t) ~ log(n))
#Coefficients:
#(Intercept) log(n)
# -9.2185 0.8646
我们还可以画出log(n)
v.s的散点图。 log(t)
,以及拟合线:
plot(log(n), log(t), main = "scatter plot")
lines(log(n), fit$fitted, col = 2, lwd = 2)
开头有一些异常值,因此斜率低于应有的值。我们现在考虑移除异常值并重新调整模型以提高稳健性。
异常值具有较大的残差。一起来看看:
plot(fit$resi, main = "residuals")
我们可以标记并删除那些异常值。看起来 0.5
是一个足够好的阈值来过滤那些异常值。
exclude <- fit$resi > 0.5
n <- n[!exclude]
t <- t[!exclude]
现在我们重新拟合线性模型并作图:
robust_fit <- lm(log(t) ~ log(n))
#Coefficients:
#(Intercept) log(n)
# -11.197 1.039
plot(log(n), log(t), main = "scatter plot (outliers removed)")
lines(log(n), robust_fit$fitted, col = 2, lwd = 2)
哇,太好了,我们是金色的!斜率估计为 1。因此 O(n)
复杂度是合理的!
如果我的时间序列有 n
个成员并且我想拟合 ARIMA(1,0,1)
模型,那么大 O
表示法的复杂度是多少?
在下面的例子中我需要知道第二行代码的复杂度:
series <- arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=1000)
result <- arima(series, order=c(1,0,1))
谢谢。
它是 O(n)
复杂度。一个故事吗?见下文。
正如我在评论中所说,我们可以通过回归模型来衡量它。作为玩具演示,请考虑以下数据收集和回归过程。
我们首先定义一个函数来测量ARMA(1,1)
模型(或ARIMA(1,0,1)
)的模型拟合时间。 (请注意,我在这里使用基本的 system.time()
。您可以考虑使用包 microbenchmark
中的 microbenchmark()
。但在下文中,我将使用相当大的 n
, 以降低时间测量的灵敏度。)
t_arma <- function(N) {
series <- arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=N)
as.numeric((system.time(arima(series, order=c(1,0,1)))))[3]
}
我们需要收集,比如 100 条数据。我们尝试 100 个越来越大 n
,并测量模型拟合时间 t
。
k <- 100; t <- numeric(k)
n <- seq(20000, by = 1000, length = k) ## start from n = 20000 (big enough)
for (i in 1:k) t[i] <- t_arma(n[i])
现在,如果我们假设复杂度是:a * (n ^ b)
,或O(n ^ b)
,我们可以通过回归模型估计a
,b
:
log(t) ~ log(a) + b * log(n)
我们对斜坡估算特别感兴趣:b
。
所以让我们调用lm()
fit <- lm(log(t) ~ log(n))
#Coefficients:
#(Intercept) log(n)
# -9.2185 0.8646
我们还可以画出log(n)
v.s的散点图。 log(t)
,以及拟合线:
plot(log(n), log(t), main = "scatter plot")
lines(log(n), fit$fitted, col = 2, lwd = 2)
开头有一些异常值,因此斜率低于应有的值。我们现在考虑移除异常值并重新调整模型以提高稳健性。
异常值具有较大的残差。一起来看看:
plot(fit$resi, main = "residuals")
我们可以标记并删除那些异常值。看起来 0.5
是一个足够好的阈值来过滤那些异常值。
exclude <- fit$resi > 0.5
n <- n[!exclude]
t <- t[!exclude]
现在我们重新拟合线性模型并作图:
robust_fit <- lm(log(t) ~ log(n))
#Coefficients:
#(Intercept) log(n)
# -11.197 1.039
plot(log(n), log(t), main = "scatter plot (outliers removed)")
lines(log(n), robust_fit$fitted, col = 2, lwd = 2)
哇,太好了,我们是金色的!斜率估计为 1。因此 O(n)
复杂度是合理的!