[PYTHON] Vorhersage des heißen Sommers mit linearem Regressionsmodell
Einführung
Meteorologische Prognostiker erstellen saisonale Temperaturvorhersagen aus Wetterkarteninformationen und Meerwassertemperaturkarten im Zusammenhang mit dem Ernino-Phänomen. Dieses Mal untersuchte ich anhand vorheriger Temperaturinformationsdaten, wie ich vorhersagen kann, ob es ein "heißer Sommer" sein wird, ohne auf diesem Expertenwissen zu beruhen.
In Bezug auf den Basistemperaturdatensatz hatte die Meteorologische Agentur zunächst Daten zu "langfristigen Änderungen der Temperatur und des Niederschlags in der Welt und in Japan" als Informationen zur "globalen Erwärmung", weshalb ich darauf Bezug nahm. Zusätzlich zur jährlichen Durchschnittstemperatur in Japan gab es Daten zur saisonalen Durchschnittstemperatur, daher habe ich mich für diese entschieden.
Die formatierte CSV-Datei enthält den folgenden Inhalt. Die saisonale Durchschnittstemperatur besteht aus Daten für Frühling (März bis Mai), Sommer (Juni bis August), Herbst (September bis November) und Winter (Dezember bis Februar).
# Data source : www.data.go.jp,,,,,,
# column1 -Jährliche durchschnittliche Temperaturabweichung in Japan (℃),,,,,,
# column2 -Jährliche durchschnittliche Niederschlagsabweichung in Japan (mm),,,,,,
# column3..6 -Japans saisonale durchschnittliche Temperaturabweichung (℃) 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
(Weggelassen)
Ab 1898 wurden Daten für mehr als 100 Jahre gesammelt und veröffentlicht. In Bezug auf die Temperaturdaten (Niederschlagsdaten) wurde anstelle der Rohdaten der Temperatur der Wert der "Abweichung vom kumulierten Durchschnittswert" angegeben.
http://www.data.jma.go.jp/cpdinfo/temp/index.html
Zunächst wurde die Temperatur (Abweichung) jeder Jahreszeit in chronologischer Reihenfolge aufgezeichnet.
** Fig. Seasonal Temperature Trend **
Es scheint, dass es von Jahr zu Jahr und von Saison zu Saison fein schwankt. Zum besseren Verständnis haben wir mit pandas.rolling_mean () einen gleitenden Durchschnitt von durchschnittlich 10 Jahren berechnet und geplottet. (Grundstück nur im Frühjahr und Sommer.)
** Fig. Seasonal Temperature Trend (moving average) **
Der Erwärmungszustand kann von der ansteigenden Linie erfasst werden. Da die beiden Linien ähnliche Bewegungen haben, kann ich auch erwarten, dass ein wenig "einen heißen Sommer vorwegnimmt".
Linear Regression Model (OLS, 1 variable)
Zunächst haben wir das folgende Modell betrachtet. ** (Temperaturabweichung in diesem Sommer) ~ (Temperaturabweichung im vorherigen Frühjahr) ** Die Regressionsanalyse wurde unter Verwendung der Statistikmodelle OLS () durchgeführt.
# 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)
** Abb. Federtemperaturabweichung (x-Achse) vs. Sommertemperaturabweichung (y-Achse) **
Dies ist ein Ergebnis mit einem geringeren Korrelationsgrad als erwartet. Die angepasste Linie zeigt vorerst eine positive Korrelation. Die zusammenfassenden () Informationen waren wie folgt.
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
==============================================================================
Wie oben erwähnt, ist die Korrelation mit R-Quadrat = 0,282 schwach.
Linear Regression Model (OLS, 2 vars, multiple)
Immerhin dachte ich, dass es nicht nur mit der Temperatur des Frühlings, sondern auch des Winters davor zusammenhängen könnte, also dachte ich über das folgende Modell nach. ** (Sommertemperatur) ~ (Wintertemperatur (vor 6 Monaten)) & (Frühlingstemperatur) **
Dies ist ein Modell für die multiple Regression. Da jede Zeile von DataFrame die Temperaturabweichung desselben Jahres enthält, wird shift () für (Wintertemperatur) so verarbeitet, dass sie sich auf die des Vorjahres bezieht.
# 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)
** Abb. Winter- und Frühlingstemperaturabweichung (x-Achse) vs. Sommertemperaturabweichung (y-Achse) **
(OLS: summer ~ b0 + b1 x winter.shift(1) + b2 x spring)
Mit R-Quadrat = 0,303 ist die Änderung gegenüber dem Vorgängermodell vernachlässigbar. Die Regressionsparameter sind wie folgt. (Auszug aus der Ausgabe von summary ().) Nur der P-Wert (0,065) von x1 ist geringfügig groß.
==============================================================================
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
==============================================================================
Aus den R-Quadrat-Werten wurde festgestellt, dass es ziemlich schwierig ist, den Temperaturtrend mit dieser Methode mit hoher Genauigkeit vorherzusagen.
PyMC3 Trial (Linear Regression)
Bis zu dem oben Gesagten können wir fast das Spiel der "Sommertemperaturvorhersage" sehen, aber die Schlussfolgerung der "schwachen Korrelation" ist langweilig. Versuchen Sie also die Bayes'sche Schätzung der Regressionsparameter mit PyMC3, der Python-Implementierung von MCMC. Ich habe es gemacht. (Es ist eine Weile her, aber ich hatte große Probleme bei der Installation von PyMC3.)
Die Problemeinstellungen sind wie folgt.
- Bayes-Schätzung des linearen Regressionsmodells. Das Modell hier ist das zweite Modell oben. (summer ~ beta0(intercept) + beta1 x winter + beta2 x spring)
- Bei Anwendung des Konzepts der Nichtinformations-Vorverteilung ergibt der Anfangswert des Varianz-Sigma eine angemessene Verteilung.
Unter Bezugnahme auf die Tutorial-Seite von PyMC3 wurde der folgende Code erstellt und simuliert.
# 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)
Für PyMC3 werden die Quelle usw. auf Github veröffentlicht, aber leider fehlt das Dokument zu diesem Zeitpunkt (August 2015) definitiv. Wenn dies mit meinem eigenen Unverständnis der Bayes'schen statistischen Theorie kombiniert wird, ist es eine Situation von doppeltem Schmerz und dreifachem Schmerz. Was ich diesmal gelernt habe ist. .. ..
--PyMC3 unterscheidet sich (in Bezug auf die Kompatibilität) von der vorherigen Version (PyMC 2.x.x). Daher ist die PyMC-Dokumentation nicht sehr hilfreich.
- Suchen Sie nach einem Beispiel-Python-Code für das, was Sie tun möchten. Auf github gibt es einen Beispielordner (obwohl nicht sehr häufig), der Codes für typische Probleme enthält.
- Wenn Sie die Details wie die Argumente von Klasse und Funktion kennen möchten, müssen Sie sich die Quelle ansehen. (Es wird gespeichert, dass die Erklärung im Header der Datei angegeben ist und die Struktur aufgrund der Funktionen von Python leicht zu erkennen ist.)
Die Ergebnisse der erhaltenen Simulation sind wie folgt.
** Fig. Traceplot **
(Ich bin nicht sicher, warum der Status von "sigma_log" ausgegeben wird ...)
** Zusammenfassung (Trace) Ausgabe **
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
(Weggelassen)
Hier wird die 95% höchste hintere Dichte (HPD) anstelle des 95% -Konfidenzintervalls (C.I.) ausgegeben, das von OLS ausgegeben wurde. (High-Density-Posterior-Distribution-Bereich, Credit-Bereich) Ich persönlich möchte eine Ausgabe mit einem etwas engeren Zeilenabstand (in einem engen Format) anfordern.
Aus diesem Simulationsergebnis wurde ein Streudiagramm erstellt.
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)
** Abb. Winter- und Frühlingstemperaturabweichung (x-Achse) vs. Sommertemperaturabweichung (y-Achse) **
Vergleich der Ergebnisse zwischen OLS und 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 |
Die Regressionsparameter waren nahezu gleich.
Zusammenfassung und zukünftige Ausgaben
Die Situation war ein wenig enttäuschend für den ursprünglichen Zweck, "einen heißen Sommer vorwegzunehmen". Die Ursache scheint zu sein, dass dieses Modell zu einfach war, aber ich denke, dass dies durch die Einführung eines Modells (AR-Modell usw.) verbessert werden kann, das Elemente der Zeitreihenanalyse enthält.
Ich habe gerade angefangen, PyMC3 zu verwenden, und ich muss noch die Theorie der Bayes'schen Statistik und die Verwendung von Modulen studieren. Es gibt viele Dinge zu tun, wie z. B. das Anpassen von Berechnungsparametern und das Bestimmen des Konvergenzstatus. (In guter Weise ist es ein unordentliches Werkzeug (Spielzeug).)
Darüber hinaus wartet der Mangel an Informationen zu PyMC3 auf die Erstellung von Dokumenten durch nahe stehende Personen. Ich möchte jedoch mein Verständnis der Modellierungs- und Analysemethoden vertiefen und Beispiele wie BUGS und Stan untersuchen, die bereits bewiesen wurden.
Referenzen (Website)
- 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/
- Meteorologische Agentur (Globale Erwärmung) http://www.data.jma.go.jp/cpdinfo/index_temp.html
- Einführung in die statistische Modellierung zur Datenanalyse (Iwanami Shoten) https://www.iwanami.co.jp/cgi-bin/isearch?isbn=ISBN978-4-00-006973-1
Recommended Posts