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 **
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) **
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) **
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) **
(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 appliquant le concept de distribution a priori non informative, la valeur initiale de la variance sigma donne une distribution appropriée.
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.
- Recherchez un exemple de code python de ce que vous voulez faire. Il existe un dossier d'exemple sur github (bien que pas très abondant), qui fournit des codes pour les problèmes typiques.
- Si vous voulez connaître les détails tels que les arguments de la classe et de la fonction, vous devez regarder la source. (Il est enregistré que l'explication est donnée dans l'en-tête du fichier et que la structure est facile à voir grâce aux fonctionnalités de python.)
Les résultats de la simulation obtenue sont les suivants.
** Fig. Traceplot **
(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) **
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)
- PyMC3 tutorial https://pymc-devs.github.io/pymc3/getting_started/
- PyMC (previous version) documentation http://pymc-devs.github.io/pymc/index.html
- StatsModels documentation http://statsmodels.sourceforge.net/stable/ --Agence météorologique (réchauffement climatique) http://www.data.jma.go.jp/cpdinfo/index_temp.html --Introduction à la modélisation statistique pour l'analyse des données (Iwanami Shoten) https://www.iwanami.co.jp/cgi-bin/isearch?isbn=ISBN978-4-00-006973-1
Recommended Posts