Skip to content

Instantly share code, notes, and snippets.

@SnowMasaya
Last active February 8, 2018 09:35
Show Gist options
  • Save SnowMasaya/5eb8878c9348cc218994e90d24b88c9d to your computer and use it in GitHub Desktop.
Save SnowMasaya/5eb8878c9348cc218994e90d24b88c9d to your computer and use it in GitHub Desktop.
時系列データの予測ライブラリ--PyFlux-- ref: https://qiita.com/GushiSnow/items/437dde3293f6d77bfa58
import numpy as np
import pandas as pd
import pyflux as pf
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv')
data.index = data['time'].values
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['sunspot.year'])
plt.ylabel('Sunspots')
plt.title('Yearly Sunspot Data')
model = pf.ARIMA(data=data, ar=4, ma=4, target='sunspot.year', family=pf.Normal())
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/MASS/drivers.csv')
data.index = data['time']
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['drivers'])
plt.ylabel('Driver Deaths')
plt.title('Deaths of Car Drivers in Great Britain 1969-84')
plt.plot()
data_d = data['drivers']
data.loc[(data['time']>=1983.05), 'seat_belt'] = 1;
data.loc[(data['time']<1983.05), 'seat_belt'] = 0;
data.loc[(data['time']>=1974.00), 'oil_crisis'] = 1;
data.loc[(data['time']<1974.00), 'oil_crisis'] = 0;
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['drivers'])
plt.ylabel('Driver Deaths')
plt.title('Deaths of Car Drivers in Great Britain 1969-84')
plt.plot()
model = pf.ARIMAX(data=data, formula='drivers~1+seat_belt',
ar=1, ma=1, family=pf.Normal())
model = pf.ARIMAX(data=data, formula='drivers~1+oil_crisis',
ar=1, ma=1, family=pf.Normal())
model = pf.ARIMAX(data=data, formula='drivers~1+seat_belt+oil_crisis',
ar=1, ma=1, family=pf.Normal())
y_{t} = Z_t\alpha_t+\epsilon_t
\alpha_{t} = T_t\alpha_{t-1}+\eta_t
\epsilon_{t} \sim N(0, \Sigma^{2}_{\epsilon})
\eta_{t} \sim N(0, \Sigma_{\eta})
\alpha_{t+1} = \alpha_{t} + K_t(y_t - Z_t\alpha_t)
x = model.fit('MLE')
y_{t} = \mu_t+\epsilon_t
\mu_{t} = \mu_{t-1}+\eta_t
\epsilon_{t} \sim N(0, \sigma^{2}_{\epsilon})
\eta_{t} \sim N(0, \sigma^{2}_{\eta})
nile = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/Nile.csv')
nile.index = pd.to_datetime(nile['time'].values,format='%Y')
plt.figure(figsize=(15,5))
plt.plot(nile.index,nile['Nile'])
plt.ylabel('Discharge Volume')
plt.title('Nile River Discharge');
plt.show()
model = pf.LocalLevel(data=nile, target='Nile', family=pf.Normal())
x = model.fit()
x.summary()
LLEV
======================================================= ==================================================
Dependent Variable: Nile Method: MLE
Start Date: 1871-01-01 00:00:00 Log Likelihood: -641.5238
End Date: 1970-01-01 00:00:00 AIC: 1287.0476
Number of observations: 100 BIC: 1292.258
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
Sigma^2 irregular 15098.5410
Sigma^2 level 1469.13678
==========================================================================================================
model.plot_fit(figsize=(15,10))
y_{t} = \sum^{P}_{i=1}\phi_{i,t}y_{t-i}+\epsilon_t
\phi_{i,t} = \phi_{i,t-1}+\eta_t
x.summary()
\epsilon_{t} \sim N(0, \sigma^{2}_{\epsilon})
\eta_{t} \sim N(0, \sigma^{2}_{\eta})
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv')
model= pf.DAR(data=data, ar=9, integ=0, target='sunspot.year')
x = model.fit('MLE')
x.summary()
y_{t+1} = \phi y_t + \theta(y_t - \mu_t)
(y_t - \mu_t)
p(y_t) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(y_t-\mu)^2}{2\sigma^2})
\log p(y_t) = \frac{1}{2}\log(2\pi\sigma^2)-\frac{1}{2}\frac{(y_t-\mu)^2}{\sigma^2}
\nabla_{\mu_t}\log p(y_t) = \frac{(y_t-\mu)^2}{\sigma^2}
\nabla^{2}_{\mu_t}\log p(y_t) = \frac{1}{\sigma^2}
\frac{\nabla_{\mu_t}\log p(y_t)}{-\nabla^{2}_{\mu_t}\log p(y_t)} = (y_t-\mu)
Normal ARIMA(4,0,4)
======================================================= ==================================================
Dependent Variable: sunspot.year Method: MLE
Start Date: 1704 Log Likelihood: -1191.0131
End Date: 1988 AIC: 2402.0263
Number of observations: 285 BIC: 2438.5512
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
Constant 11.5047 3.7053 3.1049 0.0019 (4.2422 | 18.7672)
AR(1) 1.4764 0.0738 20.0147 0.0 (1.3318 | 1.621)
AR(2) -0.2544 0.2076 -1.2257 0.2203 (-0.6613 | 0.1524)
AR(3) -0.9033 0.1831 -4.932 0.0 (-1.2622 | -0.5443)
AR(4) 0.452 0.0761 5.9379 0.0 (0.3028 | 0.6013)
MA(1) -0.2985 0.0727 -4.1047 0.0 (-0.441 | -0.156)
MA(2) -0.497 0.1008 -4.9307 0.0 (-0.6946 | -0.2995)
MA(3) 0.1824 0.0988 1.8468 0.0648 (-0.0112 | 0.376)
MA(4) 0.329 0.0837 3.9315 0.0001 (0.165 | 0.4931)
Normal Scale 15.8645
==========================================================================================================
\lambda_t = \exp(\theta_t)
p(y_t) = \frac{\exp(\theta_t)^{y_t}\exp^{-\exp(\theta_t)}}{y_t!}
\log p(y_t) = y_t\theta_t-\exp(\theta_t)-\log(y_t!)
\nabla_{\theta_t}\log p(y_t) = y_t-\exp(\theta_t)
\nabla^2_{\theta_t}\log p(y_t) = -\exp(\theta_t)
\frac{\nabla_{\mu_t}\log p(y_t)}{-\nabla^{2}_{\mu_t}\log p(y_t)} = \frac{1}{\exp(\theta_t)}(y_t-\exp(\theta_t))
leicester = pd.read_csv('http://www.pyflux.com/notebooks/leicester_goals_scored.csv')
leicester.columns= ["Time","Goals","Season2"]
plt.figure(figsize=(15,5))
plt.plot(leicester["Goals"])
plt.ylabel('Goals Scored')
plt.title('Leicester Goals Since Joining EPL');
plt.show()
model = pf.LocalLevel(data=leicester, target='Goals', family=pf.Poisson())
x = model.fit(iterations=5000)
x.summary()
10% done : ELBO is -251.57569112887904, p(y,z) is -166.3745900900033, q(z) is 85.20110103887575
20% done : ELBO is -251.0404488657637, p(y,z) is -163.08148512110355, q(z) is 87.95896374466014
30% done : ELBO is -248.3826857773143, p(y,z) is -160.95781417531765, q(z) is 87.42487160199666
40% done : ELBO is -243.78188227888018, p(y,z) is -158.15537328287218, q(z) is 85.62650899600798
50% done : ELBO is -240.84398713641178, p(y,z) is -156.96138472006595, q(z) is 83.88260241634582
60% done : ELBO is -238.7879147438848, p(y,z) is -156.14707335243622, q(z) is 82.6408413914486
70% done : ELBO is -237.27773593445806, p(y,z) is -155.37226587768643, q(z) is 81.90547005677165
80% done : ELBO is -234.66469485120308, p(y,z) is -153.9445088318954, q(z) is 80.72018601930766
90% done : ELBO is -234.21050603276024, p(y,z) is -153.86353460095933, q(z) is 80.34697143180091
100% done : ELBO is -232.41596530920148, p(y,z) is -153.02493554086158, q(z) is 79.3910297683399
Final model ELBO is -232.6451770092436
Poisson Local Level Model
======================================================= ==================================================
Dependent Variable: Goals Method: BBVI
Start Date: 0 Unnormalized Log Posterior: -153.3266
End Date: 74 AIC: 308.65314788761225
Number of observations: 75 BIC: 310.97063600114853
==========================================================================================================
Latent Variable Median Mean 95% Credibility Interval
======================================== ================== ================== =========================
Sigma^2 level 0.0404 0.0404 (0.0353 | 0.0462)
==========================================================================================================
model.plot_z(figsize=(15, 5))
model.plot_fit(figsize=(15,10))
model.plot_fit(figsize=(15, 10))
model.plot_predict_is(h=50, figsize=(15, 10))
model.plot_predict(h=20, past_values=20, figsize=(15,5))
Unnamed: 0 time drivers seat_belt oil_crisis
time
1969.000000 1 1969.000000 1687 0.0 0.0
1969.083333 2 1969.083333 1508 0.0 0.0
1969.166667 3 1969.166667 1507 0.0 0.0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment