[PYTHON] Comparaison des modèles de régression - ARMA vs Random Forest Regression

(Une addition) Dans la seconde moitié de cet article, j'ai publié un article de correction sur la régression aléatoire de Rorest. Voir aussi ici. Ce que vous ne devez pas faire lors du processus d'analyse des données de séries chronologiques (y compris la réflexion)

Lors de l'analyse de données de séries chronologiques avec Python, le module d'analyse des séries temporelles (tsa) fourni par StatsModels pour la modélisation statistique peut être la première option. De plus, la bibliothèque d'apprentissage automatique Scikit-learn a une régression d'arbre de décision et une régression de forêt aléatoire, qui effectuent des régressions basées sur différents concepts, j'ai donc voulu les essayer également et les comparer.

(Environnement de programmation: Python 2.7.11, Pandas 0.18.0, StatsModels 0.6.1, Scikit-learn 0.17.1.)

Un groupe de modèles traditionnels d'auto-retour

Lorsque vous ouvrez la documentation sur l'analyse des données de séries chronologiques, les explications apparaissent souvent dans l'ordre suivant.

Pour plus de détails sur chaque modèle, reportez-vous à Wikipedia, etc. Fondamentalement, ** ARMA ** est un modèle contenant ** AR ** et ** MA **, et ** ARIMA ** est * Puisqu'il s'agit d'un modèle qui inclut * ARMA **, si ** ARIMA ** est pris en charge en tant que bibliothèque, les quatre modèles ci-dessus peuvent être pris en charge. Cependant, StatsModels fournit ** AR **, ** ARMA **, ** ARIMA ** en tant qu'API.

http://statsmodels.sourceforge.net/stable/tsa.html

Cette fois, "Nile" a été repris comme données à analyser.

Nile River flows at Ashwan 1871-1970 This dataset contains measurements on the annual flow of the Nile as measured at Ashwan for 100 years from 1871-1970. There is an apparent changepoint near 1898.

Ce sont les données pour chaque année du débit du Nil en Afrique, et sont les données d'échantillon préparées dans la bibliothèque StatsModel.

Dans le modèle d'auto-retour, la «stationnarité» est souvent une condition préalable à l'application du modèle, mais comme le graphique ci-dessus ne montre pas d'augmentation monotone ou de diminution monotone, nous avons décidé d'appliquer le modèle ** ARMA **. Faire.

Sélection de l'ordre du modèle ARMA

Puisqu'il s'agit de données de séries chronologiques pour environ 100 ans, nous allouerons 70 ans pour la première moitié comme données pour l'ajustement (pour la formation) et 30 ans pour la seconde moitié comme données pour les tests de modèles. Il est nécessaire de sélectionner l'ordre (p, q) du modèle ARMA. Commençons par tracer le tracé ACF (Autocorrelation function plot) et PACF (Partial autocorrelation plot), également appelés cholérogrammes.

def load_data():
    df_read = sm.datasets.nile.load_pandas().data
    s_date = pd.Series(
        [pd.to_datetime(str(int(y_str))) for y_str in df_read['year']]
    )
    df = df_read.set_index(s_date)
    df = df.drop('year', axis=1)
    
    return df

df_nile = load_data()

# Plot Time Series Data
ax = df_nile['volume'].plot(figsize=(8,4), grid=True)
ax.set_xlabel('year')
ax.set_ylabel('volume x 10^8 m^3')
plt.show()

# Data Split (70: 30)
# ACF, PACF
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_train.values, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_train.values, lags=40, ax=ax2)

Fig. ACF and PACF plot of "Nile" data

Dans ARMA (p, q), q = [0, 1, 2, 3] peut être sélectionné dans le tracé ACF (en haut), et p = [0, 1] peut être sélectionné dans le tracé PACF (en bas) ( D'une certaine manière) Je comprends. Pour le moment, lorsque ARMA (p = 1, q = 3) est pris en compte et que les informations sont fournies par l'utilitaire de sélection d'ordre, le résultat est le suivant.

info_criteria = sm.tsa.stattools.arma_order_select_ic(
                    df_train.values, ic=['aic', 'bic']
                )
print(info_criteria.aic_min_order)
print(info_criteria.bic_min_order)

ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
>>> (1, 1)
>>> (1, 0)

Notez qu'il y a un avertissement concernant la convergence. Sur la base de l'avertissement, nous avons obtenu des informations selon lesquelles ARMA (1, 1) est bon à partir de l'AIC (valeur de référence de sélection de modèle) et AMMA (1, 0) est bon à partir de BIC (valeur de référence de sélection de modèle). Par conséquent, adaptons ces deux modèles aux données.

arma_10 = sm.tsa.ARMA(df_train, (1, 0)).fit()
arma_11 = sm.tsa.ARMA(df_train, (1, 1)).fit()

Cela a été traité sans erreur ni avertissement particulier. Les valeurs de retour arma_10 et arma_11 sont des objets modèle (classe ARMAResults) après ajustement. Pour le moment, j'ai également essayé des modèles d'ordre supérieur ARMA (1, 2) et ARMA (1, 3).

arma_12 = sm.tsa.ARMA(df_train, (1, 2)).fit()
arma_13 = sm.tsa.ARMA(df_train, (1, 3)).fit()

(Omis, informations de traçage, etc.)

ValueError: The computed initial AR coefficients are not stationary
You should induce stationarity, choose a different model order, or you can
pass your own start_params.

L'erreur de valeur ci-dessus s'est produite. (Bien que le message d'erreur mentionne stationnaire ...) Il y a un problème de convergence lors du calcul des paramètres, et il est difficile de l'appliquer à ARMA (1, 2) et ARMA (1, 3). y a-t-il. À l'avenir, cette cause peut être étudiée en profondeur, mais cette fois, le thème est ARMA vs Random Forest Regressor, nous envisagerons donc d'utiliser ARMA (1, 0) et ARMA (1, 1) qui correspondent. Procéder.

Validation du modèle ARMA

Vérifiez d'abord la valeur de l'AIC.

print('ARMA(1,0): AIC = %.2f' % arma_10.aic)
print('ARMA(1,1): AIC = %.2f' % arma_11.aic)
>>> ARMA(1,0): AIC = 910.94
>>> ARMA(1,1): AIC = 908.92

La valeur AIC est plus petite dans le modèle ARMA (1, 1).

Ensuite, les estimations du modèle pour la section des données de formation (70 ans) et les estimations du modèle pour la section d'essai suivante (30 ans) sont calculées.

# in-sample predict
arma_10_inpred = arma_10.predict(start='1871-01-01', end='1940-01-01')
arma_11_inpred = arma_11.predict(start='1871-01-01', end='1940-01-01')

# out-of-sample predict
arma_10_outpred = arma_10.predict(start='1941-01-01', end='1970-01-01')
arma_11_outpred = arma_11.predict(start='1941-01-01', end='1970-01-01')

Le calcul AIC et le calcul de la valeur estimée peuvent être obtenus par la méthode de la classe de résultats du modèle ARMA (classe ARMAResults). Tracez ensemble les valeurs des données d'origine et les valeurs estimées.

# plot data and predicted values
def plot_ARMA_results(origdata, pred10in, pred10out, pred11in, pred11out):
    px = origdata.index
    py1 = origdata.values
    plt.plot(px, py1, 'b:', label='orig data')
    
    px_in = pred10in.index
    plt.plot(px_in, pred10in.values, 'g')
    plt.plot(px_in, pred11in.values, 'c')
    
    px_out = pred10out.index
    plt.plot(px_out, pred10out.values, 'g', label='ARMA(1,0)')
    plt.plot(px_out, pred11out.values, 'c', label='ARMA(1,1)')
    
    plt.legend()
    plt.grid(True)
    plt.show()

plot_ARMA_results(df_nile, arma_10_inpred, arma_10_outpred, 
    arma_11_inpred, arma_11_outpred)

Fig. Nile data - original data and ARMA estimated data

La section de formation (dans l'échantillon) est de 1870 à 1940, et la section de test (hors échantillon) est après cela. C'est un peu difficile à comprendre, mais l'ARMA (1,1) semble être plus proche des données réelles. L'EQM (erreur quadratique moyenne) est obtenue à partir des résidus pour une comparaison ultérieure du modèle.

# Residue (mse for train) 
arma_10_mse_tr = np.array([r ** 2 for r in arma_10.resid]).mean()
arma_11_mse_tr = np.array([r ** 2 for r in arma_11.resid]).mean()

# Residue (mse for test)
arma_10_mse_te = np.array([(df_test.values[i] - arma_10_outpred[i]) **2 
                            for i in range(30)]).mean()
arma_11_mse_te = np.array([(df_test.values[i] - arma_11_outpred[i]) **2 
                            for i in range(30)]).mean()

Retour aléatoire dans la forêt

La base de la régression de forêt aléatoire est la régression par arbre de décision. Il semble que l'algorithme d'arbre de décision lui-même doive être compris grossièrement, mais dans cet article, l'explication est omise et la fonction de Scikit-learn est traitée comme une boîte noire. Il y a un exemple dans le document Scikit-learn, mais le résultat de son essai est montré ci-dessous.

Fig. Decision Tree Regression Example

La caractéristique est qu'il s'insère dans une ligne droite en escalier. Cette fois, ce sont des données de séries chronologiques univariées, mais je les ai essayées comme modèle pour estimer la valeur actuelle en utilisant des données passées. En particulier, Volume_current ~ Volume_Lag1 + Volume_Lag2 + Volume_Lag3 Voici le modèle. Je voulais que les données d'entraînement et les données de test soient identiques à celles de la dernière fois (70, 30), mais comme NaN apparaît dans le processus de calcul de la valeur de décalage, dropna () la première partie et définissez la longueur des données sur (67, 30). fait.

Tout d'abord, les données sont prétraitées.

df_nile['lag1'] = df_nile['volume'].shift(1) 
df_nile['lag2'] = df_nile['volume'].shift(2)
df_nile['lag3'] = df_nile['volume'].shift(3)

df_nile = df_nile.dropna()

X_train = df_nile[['lag1', 'lag2', 'lag3']][:67].values
X_test = df_nile[['lag1', 'lag2', 'lag3']][67:].values

y_train = df_nile['volume'][:67].values
y_test = df_nile['volume'][67:].values

Utilisez Random Forest Regressor de Scikit-learn.

from sklearn.ensemble import RandomForestRegressor
r_forest = RandomForestRegressor(
            n_estimators=100,
            criterion='mse',
            random_state=1,
            n_jobs=-1
)
r_forest.fit(X_train, y_train)
y_train_pred = r_forest.predict(X_train)
y_test_pred = r_forest.predict(X_test)

Les valeurs estimées sont indiquées dans la figure ci-dessous.

Fig. Nile data - original data and Random Forest Regression results

Enfin, calculez MSE.

# check residue (mse)
train_resid = y_train - y_train_pred
RFR_mse_train = np.array([r ** 2 for r in train_resid]).mean()
test_resid = y_test - y_test_pred
RFR_mse_test = np.array([r ** 2 for r in test_resid]).mean()

Comparaison de la précision du modèle

Comparez la MSE (Mean Squared Error) de chaque modèle.

(MSE : smaller ... better)

En comparant des modèles ARMA avec des ordres différents, on peut voir que ARMA (1,1) a un MSE plus petit et une meilleure précision à la fois dans la section d'entraînement et dans la section de test, ce qui est cohérent avec la comparaison des valeurs AIC obtenues précédemment. Il y a.

Dans la comparaison entre le modèle ARMA et le modèle de régression Random Forest, la MSE de la régression Random Forest est petite avec une grande différence dans la section d'apprentissage. La régression Random Forest est également petite dans la section test, mais la différence est d'environ 7%. La régression Random Forest est avantageuse en termes de précision, mais il faut tenir compte du fait que le modèle ARMA auquel le modèle statistique est appliqué peut obtenir des informations sur l'intervalle de confiance (par exemple, intervalle de confiance à 95%) en plus de la valeur médiane prédite. Est fait. Je veux l'utiliser correctement selon le but.

(En comparant les graphiques de régression ARMA et Random Forest, j'ai été surpris de nouveau que les lignes semblent très différentes. Notez également qu'aucun des modèles n'a peut-être eu d'ajustements de paramètres suffisants. S'il te plait donne moi.)

Références (site Web)

• StatsModels Documentation http://statsmodels.sourceforge.net/stable/index.html
• Scikit-learn Documentation http://scikit-learn.org/stable/
• Analyse de données chronologiques utilisables immédiatement sur le terrain (Daisuke Yokouchi (Auteur), Yoshimitsu Aoki (Auteur)) http://www.amazon.co.jp/ebook/dp/B00KNRL068/
• Python Machine Learning - O'reilly media (By Sebastian Raschka) http://shop.oreilly.com/product/9781783555130.do