<jupyter_start><jupyter_text>Introduction 複数の時系列データの時系列予測を行います。  各種ライブラリの取得<jupyter_code>%matplotlib inline 
import numpy as np
import pylab
from pandas import read_csv
from random import gauss
from statsmodels.tsa.ar_model import AR
from sklearn.metrics import roc_curve, auc
import matplotlib.pyplot as plt
import statsmodels.api as sm
import scipy.stats as st
import sys
from matplotlib import animation as ani
import sklearn.metrics as mt
import keras
from keras import backend as K
from keras.models import Model
from keras.layers import Input, LSTM, GRU, RepeatVector
from keras.layers.core import Dense, Lambda
from keras import objectives
from sklearn.preprocessing import MinMaxScaler
from copy import copy 
import time
from contextlib import contextmanager
import tensorflow as tf
from keras.callbacks import TensorBoard
from pathlib import Path

plt.style.use('ggplot')<jupyter_output><empty_output><jupyter_text>時間計測のための関数<jupyter_code>@contextmanager
def timer(title):
    t0 = time.time()
    yield
    print("{} - done in {:.10f}s".format(title, time.time() - t0))<jupyter_output><empty_output><jupyter_text>気温データの取得直接HTMLから読み込んでpandas形式にしてくれる関数です。https://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_csv.html<jupyter_code>wave_data = read_csv('https://raw.githubusercontent.com/jamesrobertlloyd/gpss-research/master/data/tsdlr_5050/daily-minimum-temperatures-in-me-train.csv', header=None, names=["Date", "Temp"])
wave_data = wave_data.sort_values(by=['Date'])
max(wave_data["Temp"].values)
production_of_gas_data = read_csv('https://raw.githubusercontent.com/jamesrobertlloyd/gpss-research/master/data/tsdlr_5050/monthly-production-of-gas-in-aus-train.csv', header=None, names=["Date", "production-of-gas"])
production_of_gas_data = production_of_gas_data.sort_values(by=['Date'])
X_orig = np.nan_to_num(np.log(wave_data["Temp"].values + 1))
X_day = wave_data["Date"].values

def show_graph(day, data, label, color="b", ylabel='temperature'):
    pylab.figure(figsize=(14, 8))
    pylab.subplot(211)
    pylab.xlabel('time')
    pylab.ylabel(ylabel)
    pylab.plot(day, data, color=color, label=label)
    pylab.legend(loc='upper right')
    pylab.show()
    
show_graph(X_day, X_orig, 'real')
X_orig_second = np.nan_to_num(np.log(production_of_gas_data["production-of-gas"].values + 1))
X_day_second = production_of_gas_data["Date"].values

show_graph(X_day_second, X_orig_second, 'real', ylabel='production-of-gas')
print(X_orig_second.shape)
print(X_orig.shape)<jupyter_output>(238,)
(500,)<jupyter_text>データの分割- 学習データ- テストデータ<jupyter_code>X_train = X_orig[:150]
X_test = X_orig[151:238]
X_train_day = X_day[:150]
X_test_day = X_day[151:238]

X_train_second = X_orig_second[:150]
X_test_second = X_orig_second[151:238]
X_train_day_second = X_day_second[:150]
X_test_day_second = X_day_second[151:238]<jupyter_output><empty_output><jupyter_text>学習データ<jupyter_code>show_graph(X_train_day, X_train, 'temperature_train')
show_graph(X_train_day_second, X_train_second, 'production_of_gas_train')<jupyter_output><empty_output><jupyter_text>テストデータ<jupyter_code>show_graph(X_test_day, X_test, 'temperature_test')
show_graph(X_test_day_second, X_test_second, 'production_of_gas_train')<jupyter_output><empty_output><jupyter_text>前処理- `get_data`: データを前処理して深層学習で扱える形に変換する- `transform_data`: データを正規化して学習が容易な形に変更する。データを正規化後の値に戻して予測結果が適切か確認する<jupyter_code>X_train_joint = np.vstack((X_train, X_train_second))
X_test_joint = np.vstack((X_test, X_test_second))

X_train_joint.shape
def get_data(data, time_steps: int=3):
    dataX = []
    print(data.shape)
    dataX = np.zeros((data.shape[0], data.shape[1], time_steps))
    
    for i in range(data.shape[0]):
        for j in range(data.shape[1] - time_steps - 1):
            dataX[i][j] = data[i, j:(j + time_steps)].T
    return np.array(dataX)

def transform_data(original_data: np.array, 
                   inverse_option: bool, 
                   scaler: object,
                   variable_number: int,
                  ):
    data_shape = original_data.shape
    print(original_data.shape)
    data = original_data.reshape(-1, variable_number)
    print(data.shape)
    if inverse_option is True:
        print('before max {}'.format(max(data[0])))
        print('Inverse')
        data = scaler.inverse_transform(data)
        print('after max {}'.format(max(data[0])))
    else:
        print('before max {}'.format(max(data[0])))
        print('Normalize')
        data = scaler.fit_transform(data)
        print('after max {}'.format(max(data[0])))
    data = data.reshape(data_shape)
    return data, scaler

def prepare_data(original_data, time_steps, variable_number):
    copy_data = original_data.copy()
    scaler = MinMaxScaler(feature_range=(0, 1), copy=False)
    data, scaler = transform_data(original_data=copy_data, 
                              inverse_option=False, scaler=scaler, variable_number=variable_number)
    data = np.asarray(data)
    x = get_data(data, time_steps=time_steps)
    x = np.swapaxes(x, 0, 1)
    x = np.swapaxes(x, 1, 2)
    return x, scaler

time_steps = 3
variable_number = 2

x, scaler = prepare_data(X_train_joint, time_steps, variable_number)


input_dim = x.shape[-1]
timesteps = x.shape[1]
print(x.shape)
print('input_dim: ', input_dim)
print('timesteps: ', timesteps)<jupyter_output>(2, 150)
(150, 2)
before max 3.126760535960395
Normalize
after max 0.3982653990092887
(2, 150)
(150, 3, 2)
input_dim:  2
timesteps:  3<jupyter_text>モデルの作成- `create_model`でモデルを作成しています。指定できるのは   - timestepsの幅   - ノードの数   - 自身の状態をリセットするかしないかを指定する`stateful`(周期的な波形は`True`の方が若干良くなります。)   - `model_option`でGRUもしくはLSTMを選択できます。   - `optimizer`で設定したい`optimizer`を指定できます。デファオルトでは`adam`が入っています。<jupyter_code>def create_model(input_dim,
                 time_steps,
                 latent_dim,
                 # データが一つしかないので1しか選べない
                 batch_size=1,
                 model_option='lstm',
                 optimizer='adam',
                ):
    with tf.name_scope('Model'):
        x = Input(shape=(time_steps, input_dim,))

        if model_option == 'lstm':
            with tf.name_scope('LSTM'):
                h = LSTM(latent_dim, stateful=False, return_sequences=True)(x)
        elif model_option == 'gru':
            with tf.name_scope('GRU'):
                h = GRU(latent_dim, stateful=False, return_sequences=True)(x)

        print('input_dim:', input_dim)
        with tf.name_scope('Dense'):
            out = Dense(input_dim)(h)

        model = Model(x, out)
        model.summary()

        with tf.name_scope('ModelCompile'):
            model.compile(optimizer=optimizer, loss='mean_squared_error', metrics=['mse'])

    return model<jupyter_output><empty_output><jupyter_text>予測関数<jupyter_code>def predict_model_show_graph(day, x, scaler, model, variable_number):
    prepare_value = x.copy()
    with timer('predict'):
        preds = model.predict(prepare_value)
    
    print('prepare_value: ', prepare_value.shape)
    print('preds: ', preds.shape)
       
    x_scale, scaler = transform_data(original_data=prepare_value, inverse_option=True, 
                                     scaler=scaler, variable_number=variable_number)

    predict_scale, scaler = transform_data(original_data=preds, inverse_option=True, 
                                           scaler=scaler, variable_number=variable_number)
    
    for i in range(prepare_value.shape[-1]):
        plt.figure(figsize=(14, 8))
        print('x_scale: ', prepare_value.shape)
        plt.plot(prepare_value[:, 0, i], color='r', label='data')
        plt.plot(preds[:, 0, i], color='b', label='predict')
        plt.legend()
        plt.show()
    
    print('model evaluate mse:', model.evaluate(preds, prepare_value))
    return preds, prepare_value
    # return predict_scale, x_scale<jupyter_output><empty_output><jupyter_text>モデルの学習<jupyter_code>predict_list = []
var_list = []

print(x.shape)
model = create_model(input_dim, 
                     time_steps=time_steps,
                     latent_dim=120,
                     model_option='lstm',
                    )

log_dir = 'simple_multi_model'
p = Path(log_dir)
p.mkdir(parents=True, exist_ok=True)
tensorboard = TensorBoard(log_dir=log_dir,
                          write_graph=True,
                          embeddings_freq=0,
                          )

with timer('train model simple'):
    model.fit(x, x, epochs=400, callbacks=[tensorboard])
window = time_steps
x_test, scaler = prepare_data(X_test_joint, time_steps, variable_number)
predict_test, x_scale_test = predict_model_show_graph(X_test_day[window + 1:], 
                                                       x_test, scaler, model,
                                                       variable_number)<jupyter_output>Epoch 1/400
150/150 [==============================] - 1s 3ms/step - loss: 0.3210 - mean_squared_error: 0.3210
Epoch 2/400
150/150 [==============================] - 0s 166us/step - loss: 0.2473 - mean_squared_error: 0.2473
Epoch 3/400
150/150 [==============================] - 0s 164us/step - loss: 0.1829 - mean_squared_error: 0.1829
Epoch 4/400
150/150 [==============================] - 0s 171us/step - loss: 0.1308 - mean_squared_error: 0.1308
Epoch 5/400
150/150 [==============================] - 0s 164us/step - loss: 0.0961 - mean_squared_error: 0.0961
Epoch 6/400
150/150 [==============================] - 0s 163us/step - loss: 0.0828 - mean_squared_error: 0.0828
Epoch 7/400
150/150 [==============================] - 0s 171us/step - loss: 0.0809 - mean_squared_error: 0.0809
Epoch 8/400
150/150 [==============================] - 0s 170us/step - loss: 0.0744 - mean_squared_error: 0.0744
Epoch 9/400
150/150 [==============================] - 0s 172us/step - loss: 0.0659 - mean_squared[...]<jupyter_text>RMSEを計算<jupyter_code>rmse_list = []
for i in range(predict_test.shape[-1]):
    print
    rmse = np.sqrt(np.average(np.subtract(x_scale_test[:, :, i], predict_test[:, :, i]) ** 2))
    print('rmse:', rmse)
    rmse_list.append(rmse)

print(np.average(rmse_list))<jupyter_output>rmse: 0.03410091276360929
rmse: 0.09107017469006462
0.06258554372683696<jupyter_text>個別入力型モデル<jupyter_code>def create_model_individual(
     input_dim,
     time_steps,
     latent_dim_list,
     # データが一つしかないので1しか選べない
     batch_size=1,
     model_option='lstm',
     optimizer='adam',
):
    input_list = []
    output_list = []
    
    with tf.name_scope('Model'):
        for i in range(input_dim):
            with tf.name_scope('Input' + str(i)):    
                x = Input(shape=(time_steps, 1,))

            if model_option == 'lstm':
                with tf.name_scope('LSTM' + str(i)):
                    h = LSTM(latent_dim_list[i], stateful=False, return_sequences=True)(x)
            elif model_option == 'gru':
                with tf.name_scope('GRU' + str(i)):
                    h = GRU(latent_dim_list[i], stateful=False, return_sequences=True)(x)
            
            with tf.name_scope('Dense' + str(i)):
                out = Dense(1)(h)
            input_list.append(x)
            output_list.append(out)

        model = Model(inputs=input_list, outputs=output_list)
        model.summary()
        with tf.name_scope('ModelCompile'):
            model.compile(optimizer=optimizer, loss='mean_squared_error', metrics=['mse'])

    return model<jupyter_output><empty_output><jupyter_text>個別入力型のためのデータ準備<jupyter_code>def prepare_data_individual(original_data, time_steps, variable_number):
    copy_data = original_data.copy()
    scaler = MinMaxScaler(feature_range=(0, 1), copy=False)
    data, scaler = transform_data(original_data=copy_data, 
                              inverse_option=False, scaler=scaler, variable_number=variable_number)
    data = np.asarray(data)
    x = get_data(data, time_steps=time_steps)
    x_reshape = []
    for i in range(x.shape[0]):
        x_reshape.append([x[i]])
    x_reshape = np.array(x_reshape)
    x_reshape = np.swapaxes(x_reshape, 1, 2)
    x_reshape = np.swapaxes(x_reshape, 2, 3)
    return list(x_reshape), scaler
   

x, scaler = prepare_data_individual(X_train_joint, time_steps, variable_number)
print(len(x))<jupyter_output>(2, 150)
(150, 2)
before max 3.126760535960395
Normalize
after max 0.3982653990092887
(2, 150)
2<jupyter_text>個別入力型モデルの予測関数<jupyter_code>def predict_model_show_graph_individual(day, x, scaler, model, variable_number):
    prepare_value = x.copy()
    with timer('predict'):
        preds = model.predict(prepare_value)
    prepare_value = np.array(prepare_value)
    preds = np.array(preds)
    
    print('prepare_value: ', prepare_value.shape)
    print('preds: ', preds.shape)
       
    x_scale, scaler = transform_data(original_data=prepare_value, inverse_option=True, 
                                     scaler=scaler, variable_number=variable_number)

    predict_scale, scaler = transform_data(original_data=preds, inverse_option=True, 
                                           scaler=scaler, variable_number=variable_number)
    
    for i in range(prepare_value.shape[0]):
        plt.figure(figsize=(14, 8))
        print('x_scale: ', prepare_value.shape)
        plt.plot(prepare_value[i, :, 0, 0], color='r', label='data')
        plt.plot(preds[i, :, 0, 0], color='b', label='predict')
        plt.legend()
        plt.show()
    
    print('model evaluate mse:', model.evaluate(list(preds), list(prepare_value)))
    return preds, prepare_value
    # return predict_scale, x_scale<jupyter_output><empty_output><jupyter_text>個別入力モデルの定義<jupyter_code>predict_list = []
var_list = []

model = create_model_individual(input_dim, 
                     time_steps=time_steps,
                     latent_dim_list=[120, 150],
                     model_option='lstm',
                    )<jupyter_output>__________________________________________________________________________________________________
Layer (type)                    Output Shape         Param #     Connected to                     
==================================================================================================
input_17 (InputLayer)           (None, 3, 1)         0                                            
__________________________________________________________________________________________________
input_18 (InputLayer)           (None, 3, 1)         0                                            
__________________________________________________________________________________________________
lstm_17 (LSTM)                  (None, 3, 120)       58560       input_17[0][0]                   
__________________________________________________________________________________________________
lstm_18 (LSTM)                  (None, 3, 150)       91200       input_18[0][0]                   
__________[...]<jupyter_text>個別入力モデルの学習<jupyter_code>log_dir = 'model_individual'
p = Path(log_dir)
p.mkdir(parents=True, exist_ok=True)
tensorboard = TensorBoard(log_dir=log_dir,
                          write_graph=True,
                          embeddings_freq=0,
                          )

with timer('train model individual'):
    model.fit(x, x, epochs=400, callbacks=[tensorboard])

window = time_steps
x_test, scaler = prepare_data_individual(X_test_joint, time_steps, variable_number)
predict_test, x_scale_test = predict_model_show_graph_individual(X_test_day[window + 1:], 
                                                       x_test, scaler, model,
                                                       variable_number)<jupyter_output>Epoch 1/400
150/150 [==============================] - 1s 6ms/step - loss: 0.7434 - dense_17_loss: 0.0684 - dense_18_loss: 0.6750 - dense_17_mean_squared_error: 0.0684 - dense_18_mean_squared_error: 0.6750
Epoch 2/400
150/150 [==============================] - 0s 227us/step - loss: 0.5646 - dense_17_loss: 0.0478 - dense_18_loss: 0.5168 - dense_17_mean_squared_error: 0.0478 - dense_18_mean_squared_error: 0.5168
Epoch 3/400
150/150 [==============================] - 0s 232us/step - loss: 0.4149 - dense_17_loss: 0.0332 - dense_18_loss: 0.3818 - dense_17_mean_squared_error: 0.0332 - dense_18_mean_squared_error: 0.3818
Epoch 4/400
150/150 [==============================] - 0s 257us/step - loss: 0.2932 - dense_17_loss: 0.0247 - dense_18_loss: 0.2685 - dense_17_mean_squared_error: 0.0247 - dense_18_mean_squared_error: 0.2685
Epoch 5/400
150/150 [==============================] - 0s 245us/step - loss: 0.2116 - dense_17_loss: 0.0223 - dense_18_loss: 0.1893 - dense_17_mean_squared_error: 0.0223 [...]<jupyter_text>精度の評価<jupyter_code>rmse_list = []
for i in range(predict_test.shape[0]):
    print
    rmse = np.sqrt(np.average(np.subtract(x_scale_test[i, :, :, :], predict_test[i, :, :, :]) ** 2))
    print('rmse:', rmse)
    rmse_list.append(rmse)

print(np.average(rmse_list))<jupyter_output>rmse: 0.017414759006381165
rmse: 0.0732443803902204
0.04532956969830078
{'fraction_non_alphanumeric': 0.14388787036489442, 'fraction_numerical': 0.04814004376367615, 'mean_word_length': 3.425543051557184, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '331180', 'license': 'Unlicense', 'n_tokens_mistral': 6579, 'n_tokens_neox': 5881, 'n_words': 938, 'path': 'advanced/time_series_anomaly_detect_keras_multivariate.ipynb', 'repo_name': 'SnowMasaya/time_series_anomaly_detect_hands_on', 'size': '16909'}