LSTM で正弦波を予測する

あゆたでデータサイエンティストとして働く佐々木です。

今回は LSTM を用いて正弦波の予測をしてみます。
正弦波のように明確な規則性のある時系列データは、LSTM で予測を行うのに向いているテーマだと言えます。そのため、LSTM の原理や実装を理解するための練習問題としては都合がいいんじゃないでしょうか。


LSTM とは

ディープラーニングで時系列データの学習を行う際に有効なものとして、
ニューラルネットワークを時間方向に展開した RNN (Recurrent newral network) があります。

しかし、従来の RNN では勾配消失問題があり、長期間の時間依存性が保持されないという問題点がありました。LSTM は勾配消失問題を解決するために改良された RNN の一種で、自然言語処理・機械翻訳などの分野で成果を上げています。ちなみに LSTM は Long Short Term Memory の略です。

原理などの詳しい説明はこちらの記事に分かりやすくまとまっているので、興味のある人は目を通してみてください。


コードを書いてみよう!

それでは、実際にコードを書いていきましょう!
今回はフレームワークとして Keras を用いて実装していきます。 Keras は Tensorflow, thenano のラッパーライブラリで、非常に簡潔な記述でモデルを構築することができます。Tensorflow 1.0 から本家の方にも導入されたので、今後はより一層使いやすくなるかもしれませんね!

定数を定義

まずはプログラム中で使う定数を宣言しておきます。
python では慣習的に、値を変えない変数は大文字で書くことになっています。

BATCH_SIZE  = 50  
HIDDEN_SIZE = 150  
NB_EPOCH    = 100  

上記の定数は上からそれぞれ、バッチサイズ・LSTMユニットのセル数・学習のエポック数です。


入力データの準備

次は入力データを作成します。 正弦波のデータとして、1周期を50分割、100周期分のデータを作成します。

def generate_sine_data():  
    STEPS_PAR_CYCLE  = 50
    NUMBER_OF_CYCLES = 100

    df = pd.DataFrame(np.arange(STEPS_PAR_CYCLE * NUMBER_OF_CYCLES + 1), columns= ['x'])
    df['sine_x'] =  df.x.apply(lambda x: math.sin(x * (2 * math.pi / STEPS_PAR_CYCLE)))
    return df['sine_x']

dataframe = generate_sine_data()  
dataset   = dataframe.values.astype('float32')  


正規化

LSTM は入力として、値の範囲を 0 から 1 の間になるよう変換したものを入れるとうまくいきます。この変換のことを正規化と呼びます。
実装は scikit-learn の MinMaxScaler を使えば簡単です。

scaler  = MinMaxScaler( feature_range=(0, 1) )  
dataset = scaler.fit_transform(dataset)  


データを学習用とテスト用に分ける

データをトレーニング用とテスト用に分割します。
ここでは前半の67%を学習用、後半の33%をテスト用として分割します。

LENGTH      = len(dataset)  
train_size  = int(LENGTH * 0.67)  
test_size   = LENGTH - train_size  
train, test = dataset[:train_size], dataset[train_size:]  


教師データを作成する

ここでは、教師用のデータを作成していきます。

def create_dataset(dataset, time_steps= 1):  
    dataX, dataY = [], []
    for i in range(len(dataset) - time_steps - 1):
        a = dataset[i : (i + time_steps)]
        dataX.append(a)
        dataY.append(dataset[i + time_steps])
    return np.array(dataX), np.array(dataY)

TIME_STEPS = 3  
trainX, trainY = create_dataset(train, TIME_STEPS)  
testX,  testY  = create_dataset(test,  TIME_STEPS)  

上記のコードで、TIME_STEPS は LSTM で何ステップ分まで過去の出力を考慮に入れて予測するかを決める変数になります。論文などでは window size や num steps などと呼ばれていることもあります。
ここでは簡単に TIME_STEPS は 3ステップとします。


データの個数を調整する

Keras の LSTM モジュールでは、バッチ処理を行う場合には入力データの個数をバッチサイズの整数倍にしなければならないという制約があります。

trainX = trainX[ len(trainX) % batch_size: ]  

上記の処理を同様に trainY, testX, testY に対しても行います。


形状を変換する

Keras の LSTM セルに入れる入力データは、 [サンプル数, time steps, 特徴量の数]の形式の三次元配列として渡す必要があります。さきほど作成した trainX, testX は [サンプル数、time steps] という二次元配列になっているので、これを3次元に変換します。

trainX = np.reshape(trainX, (trainX.shape[0], trainX.shape[1], 1))  
testX  = np.reshape(testX,  (testX.shape[0],  testX.shape[1],  1))  

コードだけみるとわかりづらいところなのですが、time step 個の要素をもつリストだったところを、特徴量の数だけの要素をもつリスト のリストにします。
文章で説明するとややこしいので、図で示してみます。

---- 変換前 ----
[
    [t1, t2, t3],
    [t2, t3, t4], 
      -- 中略 --
    [tn-2, tn-1, tn],
]

---- 変換後 ----
[
    [[t1], [t2], [t3]],
    [[t2], [t3], [t4]], 
      -- 中略 --
    [[tn-2], [tn-1], [tn]],
]

ちなみに特徴量の数が複数 (複数の時系列データを学習に使う) の場合には以下のようにすればおkです。次の例では特徴量が2つで、それぞれ A, B とします。

[
    [
        [A_t1, B_t1], 
        [A_t2, B_t2], 
        [A_t3, B_t3]
    ],
      -- 以下略 --
]


モデルを構築

さて、いよいよ本題のモデル構築をしていきます。 これに関しては Keras を用いればほんの数行で書くことができます。

model = Sequential()  
model.add( LSTM(HIDDEN_SIZE, batch_input_shape=(BATCH_SIZE, TIME_STEPS, 1)) )  
model.add( Dense(1) )  
model.compile(loss= 'mean_squared_error', optimizer= 'adam')  


LSTM()

HIDDEN_SIZE とは LSTMセルの数(下記の図における A の個数)のことで、この値が多いほどモデルの表現力が高くなると言われています。
(source: http://colah.github.io/posts/2015-08-Understanding-LSTMs/)

・batch_input_shape
これは LSTM に入力するデータの形状を指定するもので、[バッチサイズ、time steps、特徴量の数]をタプルにして入れます。なお、バッチ処理を行わない場合、ここは input_shapeという引数に変わります。その場合は入力するデータの形状をタプルで渡します。

Dense()

Dense() は単純にニューロンの数を調整するだけのレイヤーです。 ここでは、LSTM の出力を 1つの値にしています。 Dense() は擬似コードで書くと以下のような処理を行っているものです。行列の計算を行って形状を変換しているだけですね。

x : 直前の層の出力  
W : shape= [直前の層の出力数、Dense()の出力数] の重み  
b : shape= [Dense() の出力数] のバイアス

Dense() の出力 = x・W + b   # ・は行列の積  


compile()

compile() は誤差の計算や、パラメータ更新のアルゴリズムをしていするものです。
・loss
誤差の計算法を指定する引数です。ここでは、MSE(Mean Squared Error: 二乗平均誤差)を用いています。MSE は数式で示すと次のようになります。

t_iが正解の値、y_iが LSTM の出力です。真の値との2乗の差の平均を求めます。

・optimizer
パラメータ調整のアルゴリズムを指定するものです。Tensorflow でいうところの、tf.train.Adamoptimizer()等に対応するものです。ここでは、ADAM を指定しています。

学習

モデルの学習は、以下のように書けばおkです。

model.fit(trainX, trainY, nb_epoch= NB_EPOCH, batch_size= BATCH_SIZE)  

動かしてみよう!

正常に動けば、以下のような出力があるはずです。エポック毎に誤差が小さくなっており、学習が進んでいるのがわかりますね!


結果

最終的に MSE はここまで下がりました。
予測結果をグラフにすると、下の画像のようになります。


グラフを見た感じ、概ねうまくいっているように見えます。
あまりにもピッタリ重なっていて区別がつかなくなっていたので、Predict の方は意図的に1ステップ分ずらしてあります。

さて、LSTM では一見正常に学習できているように見えて、実は前のステップとほぼ同じ値を出力しているだけ。ということがあります。そこで、正常に1ステップ先を予測するように学習できているか確かめてみましょう!
※ちなみにこの問題については、このあたりで議論がされています。

これを確認するため、LSTM の予測結果を入力データとして入れ直し、再帰的に予測を繰り返すことで未来予測をしてみます。この時、正常に学習されていれば正弦波、ただ前のステップと同じ値を出しているだけならば、横ばいのグラフになるはずです。

という訳で、
学習済みのモデルを使って400ステップ先まで予測してみた結果がこちら↓ ちゃんと正弦波の未来予測ができていますね!


参考文献

[1] Time Series Prediction with LSTM Recurrent Neural Networks in Python with Keras
[2] Understanding LSTM Networks