本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程
R小盐准备介绍R语言机器学习与预测模型的学习笔记
你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】
“预测非常困难，特别是关于未来”。丹麦物理学家尼尔斯·波尔（Neils Bohr）
01 ARIMA模型
自回归移动平均模型(ARIMA, Autoregressive Integrated Moving Average)。
ARIMA是由博克思(Box)和詹金斯(Jenkins)于70年代初提出的一著名时间序列预测方法，所以又称为box-jenkins模型、博克思-詹金斯法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归, p为自回归项; MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。
ARIMA(p，d，q)模型可以表示为：
其中L是滞后算子（Lag operator），
ARIMA模型含有三个参数：p,d,q。
- p--代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项
- d--代表时序数据需要进行几阶差分化，才是稳定的，也叫Integrated项。
- q--代表预测模型中采用的预测误差的滞后数(lags)，也叫做MA/Moving Average项
ARIMA模型结合了三种基本方法： - 自回归（AR） - 在自回归的一个给定的时间序列数据在他们自己的滞后值，这是由在模型中的“P”值表示回归的值。
- 差分（I-for Integrated） - 这涉及对时间序列数据进行差分以消除趋势并将非平稳时间序列转换为平稳时间序列。这由模型中的“d”值表示。如果d = 1，则查看两个时间序列条目之间的差分，如果d = 2，则查看在d = 1处获得的差分的差分，等等。
- 移动平均线（MA） - 模型的移动平均性质由“q”值表示，“q”值是误差项的滞后值的数量。
02 ARIMA模型的R语言实现
##ARIMA模型
skirtsseriesdiff1 <- diff(skirtsseries, differences=1)
plot.ts(skirtsseriesdiff1)
skirtsseriesdiff2 <- diff(skirtsseries, differences=2)
plot.ts(skirtsseriesdiff2)
kingtimeseriesdiff1 <- diff(kingstimeseries, differences=1)
plot.ts(kingtimeseriesdiff1)
##选择合适的ARIMA模型
acf(kingtimeseriesdiff1, lag.max=20) # plot a correlogram
acf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # get the autocorrelation values
pacf(kingtimeseriesdiff1, lag.max=20) # plot a partial correlogram
pacf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # get the partial autocorrelation values
volcanodust <- scan("http://robjhyndman.com/tsdldata/annual/dvi.dat", skip=1)
volcanodustseries <- ts(volcanodust,start=c(1500))
plot.ts(volcanodustseries)
acf(volcanodustseries, lag.max=20) # plot a correlogram
acf(volcanodustseries, lag.max=20, plot=FALSE) # get the values of the autocorrelations
pacf(volcanodustseries, lag.max=20)
pacf(volcanodustseries, lag.max=20, plot=FALSE)
##使用ARIMA模型预测
kingstimeseriesarima <- arima(kingstimeseries, order=c(0,1,1)) # fit an ARIMA(0,1,1) model
kingstimeseriesarima
library("forecast") # load the "forecast" R library
kingstimeseriesforecasts <- forecast(kingstimeseriesarima, h=5)
kingstimeseriesforecasts
plot(kingstimeseriesforecasts)
acf(kingstimeseriesforecasts$residuals, lag.max=20)
Box.test(kingstimeseriesforecasts$residuals, lag=20, type="Ljung-Box")
plot.ts(kingstimeseriesforecasts$residuals) # make time plot of forecast errors
plotForecastErrors(kingstimeseriesforecasts$residuals) # make a histogram
volcanodustseriesarima <- arima(volcanodustseries, order=c(2,0,0))
volcanodustseriesarima
volcanodustseriesforecasts <- forecast(volcanodustseriesarima, h=31)
volcanodustseriesforecasts
plot(volcanodustseriesforecasts)
acf(volcanodustseriesforecasts$residuals, lag.max=20)
Box.test(volcanodustseriesforecasts$residuals, lag=20, type="Ljung-Box")
plot.ts(volcanodustseriesforecasts$residuals) # make time plot of forecast errors
plotForecastErrors(volcanodustseriesforecasts$residuals) # make a histogram
mean(volcanodustseriesforecasts$residuals)
####
参考资料：
关注R小盐，关注科研私家菜(VX_GZH: SciPrivate)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型
