Implemented and tested bootstrap idea.. may need to fix issue of outlier predictions.

Author: shz9

RNN_experiments/bootstrap/__init__.py

@@ -7,9 +7,7 @@
 """
 
 from sklearn.gaussian_process import GaussianProcessRegressor
-from sklearn.gaussian_process.kernels import ExpSineSquared
 import numpy as np
-import matplotlib.pylab as plt
 
 
 def bootstrap_data(X, Y, kernel, samples=50):
@@ -28,40 +26,19 @@ def bootstrap_data(X, Y, kernel, samples=50):
     m.fit(X, Y)
 
     boot_samples = []
-    # num_records = len(X)
+    num_records = len(X)
 
     for si in range(samples):
         # - Start of the sequence is a normal distribution centered on X[0]
         # and with standard deviation equal to 5% the value of X[0].
         # - Same logic applies to the end of the sequence.
         # - The number of steps in the sequence is equal to the number
         # of records in X +/- 10% of that number.
-        xx = np.atleast_2d(np.linspace(np.random.normal(loc=X[0], scale=.05 * X[0]),
-                                       np.random.normal(loc=X[-1], scale=.05 * X[-1]),
-                                       len(X)*2)).T
-                                       #num_records + np.random.randint(low=-int(.1 * num_records),
-                                       #                                high=int(.1 * num_records)))).T
+        xx = np.atleast_2d(np.linspace(np.random.normal(loc=X[0], scale=.0001 * X[0]),
+                                       np.random.normal(loc=X[-1], scale=.0001 * X[-1]),
+                                       num_records + np.random.randint(low=-int(.1 * num_records),
+                                                                       high=int(.1 * num_records)))).T
         yy = m.predict(xx)
-        #print yy
         boot_samples.append((xx, yy))
 
     return boot_samples
-
-
-def test_bootstrap():
-
-    def f(x):
-        return 2. * x * np.sin(x)
-
-    X_ = np.atleast_2d(np.linspace(1., 50., 20)).T
-    Y_ = f(X_).ravel()
-
-    boots = bootstrap_data(X_, Y_, ExpSineSquared(), samples=15)
-
-    for bt in boots:
-        plt.plot(bt[0], bt[1])
-
-    plt.show()
-
-
-test_bootstrap()

RNN_experiments/bootstrap/gp_bootstrap.py

@@ -1,13 +1,15 @@
-# src : https://machinelearningmastery.com/time-series-forecasting-long-short-term-memory-network-python/
 import pandas as pd
-from sklearn.metrics import mean_squared_error
 from sklearn.preprocessing import MinMaxScaler
 from keras.models import Sequential
 from keras.layers import Dense
 from keras.layers import LSTM
-from math import sqrt
-from matplotlib import pyplot
-import numpy
+from keras.models import load_model, save_model
+from matplotlib import pyplot as plt
+import numpy as np
+import glob
+from joblib import Parallel, delayed
+from sklearn.gaussian_process.kernels import RBF, WhiteKernel, RationalQuadratic, ExpSineSquared
+from RNN_experiments.bootstrap.gp_bootstrap import bootstrap_data
 
 
 # frame a sequence as a supervised learning problem
@@ -19,6 +21,7 @@ def timeseries_to_supervised(data, lag=1):
     df.fillna(0, inplace=True)
     return df
 
+
 # scale train and test data to [-1, 1]
 def scale(train, test):
     # fit scaler
@@ -36,7 +39,7 @@ def scale(train, test):
 # inverse scaling for a forecasted value
 def invert_scale(scaler, X, value):
     new_row = [x for x in X] + [value]
-    array = numpy.array(new_row)
+    array = np.array(new_row)
     array = array.reshape(1, len(array))
     inverted = scaler.inverse_transform(array)
     return inverted[0, -1]
@@ -61,71 +64,118 @@ def fit_lstm(train, batch_size, nb_epoch, neurons):
 def forecast_lstm(model, batch_size, X):
     X = X.reshape(1, 1, len(X))
     yhat = model.predict(X, batch_size=batch_size)
-    return yhat[0,0]
-
-
-# load dataset
-ex_dataset = pd.read_csv('../../data/mauna-loa-atmospheric-co2.csv',
-                         header=None)
-ex_dataset.columns = ['CO2Concentration', 'Time']
-
-print ex_dataset['CO2Concentration'].diff()
-
-# transform data to be stationary
-raw_values = series.values
-diff_values = difference(raw_values, 1)
-# transform data to be supervised learning
-supervised = timeseries_to_supervised(diff_values, 1)
-supervised_values = supervised.values
-
-# split data into train and test-sets
-train, test = supervised_values[0:-228], supervised_values[-228:]
-# transform the scale of the data
-# scaler, train_scaled, test_scaled = scale(train, test)
-
-# fit the model
-# lstm_model = fit_lstm(train_scaled, 1, 1000, 4)
-lstm_model = fit_lstm(train, 1, 1000, 4)
-
-# forecast the entire training dataset to build up state for forecasting
-train_reshaped = train[:, 0].reshape(len(train), 1, 1)
-lstm_model.predict(train_reshaped, batch_size=1)
-
-# walk-forward validation on the test data
-predictions = list()
-prev = None
-prev_history = [raw_values[len(test)+1]]  # initial
-for i in range(len(test)):
-    # make one-step forecast
-    if prev is None:
-        prev = test[i, 0:-1]
-    print('Start---------------------')
-    print('prev=', prev)
-    yhat = forecast_lstm(lstm_model, 1, prev)
-    prev = yhat
-    print('forcast=', yhat)
-    # reshape
-    prev = numpy.array([prev])
-    # invert scaling
-    # yhat = invert_scale(scaler, X, yhat)
-    # print('yhat after inverse scale: ', yhat)
-    # invert differencing
-    # yhat = inverse_difference(prev_history, yhat, i)
-    yhat = yhat + prev_history[i]
-    print('value added=', prev_history[i])
-    print('prediction=', yhat)
-    prev_history.append(yhat)
-    # store forecast
-    predictions.append(yhat)
-    expected = raw_values[len(train) + i +1]
-    print('End---------------')
-    print('Month=%d, Predicted=%f, Expected=%f' % (i+1, yhat, expected))
-
-# report performance
-rmse = sqrt(mean_squared_error(raw_values[-228:], predictions))
-print('Test RMSE: %.3f' % rmse)
-# line plot of observed vs predicted
-pyplot.plot(raw_values[-228:])
-pyplot.plot(predictions)
-pyplot.savefig('test.eps')
-# pyplot.show()
+    return yhat[0, 0]
+
+
+def train_model(ith_dataset, idx):
+    diff_values = ith_dataset['CO2Concentration'].diff()
+
+    # transform data to be supervised learning
+    supervised = timeseries_to_supervised(diff_values, 1)
+
+    lstm_model = fit_lstm(supervised.values, 1, 50, 4)
+
+    save_model(lstm_model, "./models/" + str(idx) + ".kmodel")
+
+
+def make_model_predictions(model, train, all_y, init_val):
+
+    # forecast the entire training dataset to build up state for forecasting
+    train_reshaped = train[:, 0].reshape(len(train), 1, 1)
+
+    model.predict(train_reshaped, batch_size=1)
+
+    # walk-forward validation on the test data
+    predictions = []
+    prev = None
+    prev_history = [init_val]  # initial
+    for i in range(len(all_y)):
+        # make one-step forecast
+        if prev is None:
+            prev = all_y[i, 0:-1]
+
+        yhat = forecast_lstm(model, 1, prev)
+        prev = yhat
+        # reshape
+        prev = np.array([prev])
+        # invert scaling
+        # yhat = invert_scale(scaler, X, yhat)
+        # print('yhat after inverse scale: ', yhat)
+        # invert differencing
+        # yhat = inverse_difference(prev_history, yhat, i)
+        yhat = yhat + prev_history[i]
+        prev_history.append(yhat)
+        # store forecast
+        predictions.append(yhat)
+
+    return predictions
+
+
+def main(retrain=True):
+
+    # load dataset
+    ex_dataset = pd.read_csv('../../data/mauna-loa-atmospheric-co2.csv',
+                             header=None)
+    ex_dataset.columns = ['CO2Concentration', 'Time']
+
+    train_data = ex_dataset.loc[ex_dataset.Time <= 1980, ['CO2Concentration', 'Time']]
+
+    if retrain:
+        bootstrapped_dataset = bootstrap_data(train_data['Time'].reshape(-1, 1),
+                                              train_data['CO2Concentration'].reshape(-1, 1),
+                                              34.4**2 * RBF(length_scale=41.8) +
+                                              3.27**2 * RBF(length_scale=180) * ExpSineSquared(length_scale=1.44,
+                                                                                               periodicity=1) +
+                                              0.446**2 * RationalQuadratic(alpha=17.7, length_scale=0.957) +
+                                              0.197**2 * RBF(length_scale=0.138) + WhiteKernel(noise_level=0.0336),
+                                              samples=10)
+        # Need to run this in parallel:
+
+        # t_pool = ThreadPool(20)
+        Parallel(n_jobs=10)(delayed(train_model)(pd.DataFrame({'Time': np.ravel(dat[0]),
+                                                               'CO2Concentration': np.ravel(dat[1])}),
+                                                 idx)
+                            for idx, dat in enumerate(bootstrapped_dataset))
+
+    """for x, y in bootstrapped_dataset:
+        temp_df = pd.DataFrame({'Time': np.ravel(x),
+                                'CO2Concentration': np.ravel(y)})
+        rnn_models.append(train_model(temp_df))"""
+
+    rnn_models = []
+    for mod_path in glob.glob("./models/*.kmodel"):
+        rnn_models.append(load_model(mod_path))
+
+    diff_values = ex_dataset['CO2Concentration'].diff()
+    # transform data to be supervised learning
+    supervised = timeseries_to_supervised(diff_values, 1)
+    supervised_values = supervised.values
+
+    # split data into train and test-sets
+    train, test = supervised_values[0:-228], supervised_values[-228:]
+
+    preds = []
+    for mod in rnn_models:
+        preds.append(make_model_predictions(mod, train, test, ex_dataset['CO2Concentration'][len(test)+1]))
+
+    # line plot of observed vs predicted
+    rnn_means = np.array([])
+    rnn_conf = np.array([])
+
+    for k in range(len(preds[0])):
+        step_vals = [el[k] for el in preds]
+        rnn_means = np.append(rnn_means, np.mean(step_vals))
+        rnn_conf = np.append(rnn_conf, np.std(step_vals))
+
+    plt.plot(ex_dataset['CO2Concentration'][-228:])
+    plt.plot(rnn_means)
+    plt.fill_between(list(range(len(rnn_means))),
+                     rnn_means - rnn_conf,
+                     rnn_means + rnn_conf,
+                     color="gray", alpha=0.2)
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main(False)