[PYTHON] Prédire l'été chaud avec un modèle de régression linéaire

introduction

Les prévisionnistes météorologiques font des prévisions de température saisonnière à partir d'informations cartographiques météorologiques et de cartes de température de l'eau de mer liées au phénomène Ernino. Cette fois, j'ai étudié comment prédire si ce sera un "été chaud" à partir des données de température passées, non basées sur de telles connaissances d'experts.

Tout d'abord, en ce qui concerne l'ensemble de données de base sur la température, l'Agence météorologique disposait de données sur "les changements à long terme des températures et des précipitations dans le monde et au Japon" en tant qu'informations sur le "réchauffement planétaire", alors j'y ai fait référence. En plus de la température moyenne annuelle au Japon, il y avait des données sur la température moyenne saisonnière, alors j'ai décidé de l'utiliser.

Le fichier csv formaté a le contenu suivant. La température moyenne saisonnière est constituée de données pour le printemps (mars à mai), l'été (juin à août), l'automne (septembre à novembre) et l'hiver (décembre à février).

# Data source : www.data.go.jp,,,,,, 
# column1 -Écart annuel moyen de température au Japon (℃),,,,,, 
# column2 -Écart annuel moyen des précipitations au Japon (mm),,,,,, 
# column3..6 -Écart saisonnier de la température moyenne au Japon (℃) Mar-May, Jun-Aug, Sep-Nov, Dec-Feb,,,,,, 
Year,temp_yr,rain_yr,spring,summer,autumn,winter 
1898,-0.75,15.1,-0.98,-1.7,-0.53,-0.67
1899,-0.81,199.2,-0.27,-0.21,-0.6,-2.12
1900,-1.06,-43.3,-1.35,-1.09,-1.11,-0.41
1901,-1.03,48.6,-0.74,-0.89,-1.26,-0.87
1902,-1.03,154.7,-1.65,-0.87,-2.2,-0.43
1903,-0.77,266.2,0.73,-0.19,-1.5,-0.93
1904,-0.86,-48.6,-1.37,-0.75,-0.16,-1.68
1905,-0.95,256.9,-0.55,-1.45,-1.39,-0.86
(Omis)

Les données ont été accumulées et rendues publiques pendant plus de 100 ans à partir de 1898. En ce qui concerne les données de température (précipitations), la valeur de «l'écart par rapport à la valeur moyenne cumulée» a été divulguée à la place des données brutes de la température.

http://www.data.jma.go.jp/cpdinfo/temp/index.html

Tout d'abord, la température (écart) de chaque saison a été tracée par ordre chronologique. ** Fig. Seasonal Temperature Trend ** temp_TL1.png

Il semble qu'il fluctue finement d'année en année et de saison en saison. Pour faciliter la compréhension, nous avons calculé une moyenne mobile de 10 ans en moyenne avec pandas.rolling_mean () et l'avons tracée. (Terrain seulement au printemps et en été.)

** Fig. Seasonal Temperature Trend (moving average) ** temp_TL2.png

L'état de réchauffement peut être saisi par la ligne montante. Aussi, comme les deux lignes ont des mouvements similaires, je sens que je peux m'attendre un peu à "anticiper un été chaud".

Linear Regression Model (OLS, 1 variable)

Tout d'abord, nous avons considéré le modèle suivant. ** (Écart de température cet été) ~ (Écart de température au printemps précédent) ** L'analyse de régression a été réalisée à l'aide de statsmodels OLS ().

# Regression Analysis (OLS)

import statsmodels.api as sm 
x1 = tempdf['spring'].values 
x1a = sm.add_constant(x1) 
y1 = tempdf['summer'].values 

model1 = sm.OLS(y1, x1a) 
res1 = model1.fit() 
print res1.summary() 

plt.figure(figsize=(5,4)) 
plt.scatter(tempdf['spring'], tempdf['summer'], c='b', alpha=0.6)
plt.plot(x1, res1.fittedvalues, 'r-', label='OLS(model1)') 
plt.grid(True)

** Fig. Écart de température au printemps (axe x) par rapport à l'écart de température en été (axe y) ** temp_trend_scatter_model1.png

Il s'agit d'un résultat avec un degré de corrélation plus faible que prévu. La ligne ajustée montre une corrélation positive pour le moment. Les informations résumées () étaient les suivantes.

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.282
Model:                            OLS   Adj. R-squared:                  0.276
Method:                 Least Squares   F-statistic:                     45.26
Date:                Mon, 10 Aug 2015   Prob (F-statistic):           7.07e-10
Time:                        16:39:14   Log-Likelihood:                -107.55
No. Observations:                 117   AIC:                             219.1
Df Residuals:                     115   BIC:                             224.6
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         -0.3512      0.068     -5.148      0.000        -0.486    -0.216
x1             0.4489      0.067      6.727      0.000         0.317     0.581
==============================================================================
Omnibus:                        0.159   Durbin-Watson:                   1.571
Prob(Omnibus):                  0.924   Jarque-Bera (JB):                0.075
Skew:                           0.062   Prob(JB):                        0.963
Kurtosis:                       2.990   Cond. No.                         1.88
==============================================================================

Comme mentionné ci-dessus, la corrélation est faible avec R-carré = 0,282.

Linear Regression Model (OLS, 2 vars, multiple)

Après tout, je pensais que cela pouvait être lié à la température non seulement au printemps mais aussi en hiver avant cela, alors j'ai envisagé le modèle suivant. ** (Température d'été) ~ (Température d'hiver (il y a 6 mois)) & (Température de printemps) **

Ceci est un modèle de régression multiple. Étant donné que chaque ligne de DataFrame contient l'écart de température de la même année, shift () est traité pour (température d'hiver) pour faire référence à l'année précédente.

# Regression Analysis (OLS: summer ~ b0 + b1*winter.shift(1) + b2*spring)
tempdf['winter_last'] = tempdf['winter'].shift(1)
X2 = tempdf[['winter_last', 'spring']].iloc[1:,:].values
X2a = sm.add_constant(X2)
y2 = tempdf['summer'].iloc[1:,].values

model2 = sm.OLS(y2, X2a)
res2 = model2.fit()
print res2.summary()

x2_model2 = res2.params[1] * tempdf['winter_last'] + res2.params[2] *tempdf['spring']
plt.figure(figsize=(5,4))
plt.scatter(x2_model2[1:], y2, c='g', alpha=0.6)
plt.plot(x2_model2[1:], res2.fittedvalues, 'r-', label='OLS(model2)')
plt.grid(True)

** Fig. Écart de température en hiver et au printemps (axe x) vs écart de température en été (axe y) ** temp_trend_scatter_model2.png

(OLS: summer ~ b0 + b1 x winter.shift(1) + b2 x spring)

Avec R-carré = 0,303, le changement par rapport au modèle précédent est négligeable. Les paramètres de régression sont les suivants. (Extrait de la sortie summary ().) Seule la valeur P (0,065) de x1 est légèrement élevée.

==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         -0.2927      0.073     -4.007      0.000        -0.437    -0.148
x1             0.1609      0.086      1.863      0.065        -0.010     0.332
x2             0.3984      0.070      5.674      0.000         0.259     0.537
==============================================================================


À partir des valeurs R-carré, il a été constaté qu'il est assez difficile de prédire la tendance de la température avec une grande précision par cette méthode.

PyMC3 Trial (Linear Regression)

Jusqu'à ce qui précède, nous pouvons presque voir le jeu de la "prévision de la température estivale", mais la conclusion de la "faible corrélation" est ennuyeuse, alors essayez l'estimation bayésienne des paramètres de régression à l'aide de PyMC3, qui est l'implémentation Python de MCMC. Je l'ai fait. (C'était il y a quelque temps, mais j'ai eu beaucoup de mal à installer PyMC3.)

Les paramètres du problème sont les suivants. --Estimation de Bayes du modèle de régression linéaire. Le modèle ici est le deuxième modèle ci-dessus. (summer ~ beta0(intercept) + beta1 x winter + beta2 x spring)

En référence à la page Tutorial de PyMC3, le code suivant a été préparé et simulé.

# using PyMC3

import pymc3 as pm

model3 = pm.Model()

myX1 = tempdf['winter_last'].iloc[1:,].values
myX2 = tempdf['spring'].iloc[1:,].values
myY = y2

# model definition (step.1)
with model3:
    # Priors for unknown model parameters
    intercept = pm.Normal('intercept', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10, shape=2)
    sigma = pm.HalfNormal('sigma', sd=1)
    # Expected value of outcome
    mu = intercept + beta[0]*myX1 + beta[1]*myX2
    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=myY)

# run simulation (step.2)
from scipy import optimize
with model3:
    # obtain starting values via MAP
    start = pm.find_MAP(fmin=optimize.fmin_powell)
    # instantiate sampler
    step = pm.NUTS(scaling=start)
    # draw 500 posterior samples
    trace = pm.sample(500, step, start=start)

pm.summary(trace)
pm.traceplot(trace)

Pour PyMC3, les sources etc. sont publiées sur github, mais malheureusement, le document fait définitivement défaut pour le moment (août 2015). Si cela se combine avec mon propre manque de compréhension de la théorie statistique bayésienne, c'est une situation de double douleur et de triple douleur. Ce que j'ai appris cette fois, c'est. .. ..

--PyMC3 est différent (en termes de compatibilité) de la version précédente (PyMC 2.x.x). Par conséquent, la documentation PyMC n'est pas très utile.

Les résultats de la simulation obtenue sont les suivants.

** Fig. Traceplot ** temp_trend_traceplot1.png

(Je ne sais pas pourquoi le statut de "sigma_log" est affiché ...)

** Sortie récapitulative (trace) **

intercept:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  -0.286           0.080            0.005            [-0.431, -0.127]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  -0.435         -0.341         -0.286         -0.235         -0.129


beta:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.161            0.095            0.005            [-0.024, 0.347]
  0.407            0.077            0.004            [0.266, 0.568]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  -0.017         0.095          0.159          0.222          0.361
  0.252          0.362          0.405          0.459          0.559

(Omis)

Ici, la densité postérieure la plus élevée de 95% (HPD) est sortie au lieu de l'intervalle de confiance de 95% (C.I.) qui a été produit par OLS. (Section distribution postérieure haute densité, section crédit) Personnellement, je voudrais demander une sortie avec un espacement des lignes un peu plus étroit (dans un format serré).

À partir de ce résultat de simulation, un nuage de points a été dessiné.

myX = trace['beta'][-1][0] *myX1 + trace['beta'][-1][1] *myX2

plt.figure(figsize=(5,4))
plt.scatter(myX, myY, c='g', alpha=0.6)
plt.grid(True)

for i in range(len(trace)):
    plt.plot(myX, trace['beta'][i][0] *myX1 + trace['beta'][i][1] * myX2, color='blue', alpha=0.005)

** Fig. Écart de température en hiver et au printemps (axe x) vs écart de température en été (axe y) ** temp_trend_scatter_model3.png

Comparaison des résultats entre OLS et MCMC (PyMC3)

intercept beta1 beta2
OLS -0.293 0.161 0.398
MCMC(500step) -0.286 0.161 0.407
MCMC(10000step) -0.291 0.163 0.397

Les paramètres de régression étaient presque les mêmes.

Résumé et futurs numéros

La situation était un peu décevante dans le but initial «d'anticiper un été chaud». La cause semble être que ce modèle était trop simple, mais je pense que cela peut être amélioré en introduisant un modèle (modèle AR, etc.) qui incorpore des éléments d'analyse de séries chronologiques.

Je viens de commencer à utiliser PyMC3, et je dois encore étudier la théorie des statistiques bayésiennes et comment utiliser les modules. Il y a de nombreuses choses à faire, comme l'ajustement des paramètres de calcul et la détermination de l'état de convergence. (Dans le bon sens, c'est un outil (jouet) en désordre.)

De plus, le manque d'informations sur PyMC3 attend la préparation de documents par des parties liées, mais je souhaite approfondir ma compréhension des méthodes de modélisation et d'analyse en examinant des exemples tels que BUGS et Stan qui ont déjà fait leurs preuves.

Références (site Web)

Recommended Posts

Prédire l'été chaud avec un modèle de régression linéaire
Régression avec un modèle linéaire
Régression linéaire avec statsmodels
Implémenter un modèle de régression logistique en temps discret avec stan
[Python] Régression linéaire avec scicit-learn
Régression linéaire robuste avec scikit-learn
Régression linéaire avec distribution t de Student
Créer un itérateur de modèle avec PySide
Essayez de déduire à l'aide d'un modèle de régression linéaire sur Android [PyTorch Mobile]
<Cours> Machine learning Chapitre 1: Modèle de régression linéaire
Régression linéaire
Implémenter un modèle avec état et comportement
Prédire les épidémies de maladies infectieuses avec le modèle SIR
Essayez TensorFlow RNN avec un modèle de base
Un modèle qui identifie la guitare avec fast.ai
Optimisation apprise avec OR-Tools [Planification linéaire: modèle en plusieurs étapes]
Modèle de régression multivariée avec scikit-learn - J'ai essayé de comparer et de vérifier SVR
Simulez une bonne date de Noël avec un modèle optimisé Python
Implémentation python de la classe de régression linéaire bayésienne
Introduction à l'hypothèse Tensorflow-About et au coût de la régression linéaire
Créer une visionneuse de modèle 3D avec PyQt5 et PyQtGraph