在 R 中使用 arima 循环
Looping with arima in R
我正在尝试使用 for 函数执行多个 arimas。
到目前为止我的尝试是这样的。
for(p in 0:20){
for(q in 0:20){
for (d in 0:3) {
fit <- arima(y, order=c(p,d,q),method="ML")
acc <- accuracy(fit)
print(p);print(d);print(q)
}
}
}
我想获得一个数据集中的所有精度向量,对于每个 arima,有 3 个额外的列,p、d、q。
然后我想挽救模型和 AIC 的对数似然
所以最终输出应该是一个看起来像这样的数据框,其中每一行都是一个模型。
ME RMSE MAE MPE MAPE MASE ACF1 loglikeli AIC p d q
Training set x x x x x x x x x x x x
Training set w w w w w w w w w w w w
Training set y y y y y y y y y y y y
<span class="math-container">```</span>
你快到了。最简单的方法是迭代地将每个 Arima 对象的摘要添加到数据框中。
下面的代码可以满足您的需求(我已经减少了迭代次数,否则会花费太长时间)。
# load package
library(forecast)
# create some data
y <- rnorm(100)
# object to store arima summary in
model_smry <- data.frame()
# loop to store
for(p in 1:2){
for(q in 1:2){
for (d in 1:2) {
fit <- arima(y, order=c(p,d,q),method="ML")
acc <- accuracy(fit)
# gather everything into a single data frame
acc_ext <- data.frame(# information from accuracy function
acc,
# goodness of fit
loglikeli = logLik(fit),
AIC = AIC(fit),
# arima order
p,
q,
d)
# add arima summary
model_smry <- rbind(model_smry, acc_ext)
}
}
}
# show summary
model_smry
输出:
ME RMSE MAE MPE MAPE MASE ACF1 loglikeli AIC p q d
Training set 0.03590650 0.8270888 0.6536260 61.96955 124.9386 0.6845685 -0.006412280 -124.4806 254.9612 1 1 1
Training set1 -0.03384812 0.9791048 0.7565525 -540.90163 825.9039 0.7923675 -0.129324621 -140.8802 287.7604 1 1 2
Training set2 0.03709185 0.8225973 0.6502470 95.99749 134.0075 0.6810295 0.026597486 -123.9961 255.9921 1 2 1
Training set3 -0.04914004 0.8317765 0.6574596 -51.59546 250.4728 0.6885836 -0.013578522 -129.2061 266.4123 1 2 2
Training set4 0.03698832 0.8239479 0.6516438 26.90046 162.6580 0.6824924 0.001452607 -124.0094 256.0188 2 1 1
Training set5 -0.04342442 0.9527430 0.7233051 -39.01621 319.2866 0.7575462 -0.050439230 -138.4554 284.9108 2 1 2
Training set6 0.03606286 0.8227565 0.6522152 -30.19270 220.3092 0.6830908 -0.003839680 -123.8827 257.7654 2 2 1
Training set7 -0.05099161 0.8291406 0.6503652 -91.41328 289.5055 0.6811533 -0.004315754 -128.5307 267.0613 2 2 2
除了原始问题和@ralph 的答案中的嵌套循环,还可以使用一个函数,然后用 apply
将其矢量化。这种风格模块化更好,更容易调试:
library(forecast)
# create some data
set.seed(456)
y <- rnorm(100)
# create a matrix of all desired orders
pqd <- expand.grid(p=1:2, q=1:2, d=1:2)
# a function that does the analysis for a single case
fit_arima <- function(ord) {
fit <- arima(y, order = ord, method = "ML")
acc <- accuracy(fit)
c(acc = acc, loglikeli = logLik(fit), AIC = AIC(fit))
}
# test a single case
fit_arima(c(p=1, q=1, d=1))
# run all
ret <- apply(pqd, 1, fit_arima)
# bind input and results together
cbind(pqd, t(ret))
我正在尝试使用 for 函数执行多个 arimas。
到目前为止我的尝试是这样的。
for(p in 0:20){
for(q in 0:20){
for (d in 0:3) {
fit <- arima(y, order=c(p,d,q),method="ML")
acc <- accuracy(fit)
print(p);print(d);print(q)
}
}
}
我想获得一个数据集中的所有精度向量,对于每个 arima,有 3 个额外的列,p、d、q。
然后我想挽救模型和 AIC 的对数似然
所以最终输出应该是一个看起来像这样的数据框,其中每一行都是一个模型。
ME RMSE MAE MPE MAPE MASE ACF1 loglikeli AIC p d q
Training set x x x x x x x x x x x x
Training set w w w w w w w w w w w w
Training set y y y y y y y y y y y y
<span class="math-container">```</span>
你快到了。最简单的方法是迭代地将每个 Arima 对象的摘要添加到数据框中。
下面的代码可以满足您的需求(我已经减少了迭代次数,否则会花费太长时间)。
# load package
library(forecast)
# create some data
y <- rnorm(100)
# object to store arima summary in
model_smry <- data.frame()
# loop to store
for(p in 1:2){
for(q in 1:2){
for (d in 1:2) {
fit <- arima(y, order=c(p,d,q),method="ML")
acc <- accuracy(fit)
# gather everything into a single data frame
acc_ext <- data.frame(# information from accuracy function
acc,
# goodness of fit
loglikeli = logLik(fit),
AIC = AIC(fit),
# arima order
p,
q,
d)
# add arima summary
model_smry <- rbind(model_smry, acc_ext)
}
}
}
# show summary
model_smry
输出:
ME RMSE MAE MPE MAPE MASE ACF1 loglikeli AIC p q d
Training set 0.03590650 0.8270888 0.6536260 61.96955 124.9386 0.6845685 -0.006412280 -124.4806 254.9612 1 1 1
Training set1 -0.03384812 0.9791048 0.7565525 -540.90163 825.9039 0.7923675 -0.129324621 -140.8802 287.7604 1 1 2
Training set2 0.03709185 0.8225973 0.6502470 95.99749 134.0075 0.6810295 0.026597486 -123.9961 255.9921 1 2 1
Training set3 -0.04914004 0.8317765 0.6574596 -51.59546 250.4728 0.6885836 -0.013578522 -129.2061 266.4123 1 2 2
Training set4 0.03698832 0.8239479 0.6516438 26.90046 162.6580 0.6824924 0.001452607 -124.0094 256.0188 2 1 1
Training set5 -0.04342442 0.9527430 0.7233051 -39.01621 319.2866 0.7575462 -0.050439230 -138.4554 284.9108 2 1 2
Training set6 0.03606286 0.8227565 0.6522152 -30.19270 220.3092 0.6830908 -0.003839680 -123.8827 257.7654 2 2 1
Training set7 -0.05099161 0.8291406 0.6503652 -91.41328 289.5055 0.6811533 -0.004315754 -128.5307 267.0613 2 2 2
除了原始问题和@ralph 的答案中的嵌套循环,还可以使用一个函数,然后用 apply
将其矢量化。这种风格模块化更好,更容易调试:
library(forecast)
# create some data
set.seed(456)
y <- rnorm(100)
# create a matrix of all desired orders
pqd <- expand.grid(p=1:2, q=1:2, d=1:2)
# a function that does the analysis for a single case
fit_arima <- function(ord) {
fit <- arima(y, order = ord, method = "ML")
acc <- accuracy(fit)
c(acc = acc, loglikeli = logLik(fit), AIC = AIC(fit))
}
# test a single case
fit_arima(c(p=1, q=1, d=1))
# run all
ret <- apply(pqd, 1, fit_arima)
# bind input and results together
cbind(pqd, t(ret))