LSTM Seq2Seqモデル実装

概要

株価、為替、天気、動画など時系列データの予測でよく使われるディープラーニングの代表的手法RNN(再帰型ニューラルネットワーク)の拡張バージョンに、LSTM(Long short-term memory)と呼ばれるモデルがある。今回は、LSTM Seq2Seq、Many to Manyモデルを実装して、円の第4象限の一部(訓練未実施)に対して、時系列データの予測を行ってみる。

環境

keras/LSTM, Google Colab CPU

モデル

lstm_seq2seq_model
lstm_seq2seq_model

以下モデルの概要を説明する。
・Sequentialは、あるレイヤの全ノードと、次レイヤの全ノードをつなぐRNNのモデル。
・RNNの第1レイヤとして、LSTMを追加する。第1引数は出力の次元数、activationは活性化関数で、input_shapeは入力データのフォーマット(in_time_steps x in_features)となる。
・RepeatVectorにより、出力を入力として繰り返す。ここでの繰り返し回数は予測範囲(out_vectors)となる。
・再度のLSTM、ただし、ここではシーケンス出力まさにmany outputつまり、return_sequencesをTrueに指定する。
・TimeDistributedを指定し、かつ、Dense(out_features)で、出力の次元数を指定する。
・最後にcompileメソッドで、学習時の最適化手法や、損失関数を指定する。ここでは最適化手法としてAdamを、損失関数としてMSE(Mean Squared Errorつまり平均2乗誤差)に指定する。

実装

# -*- coding: utf-8 -*-
# Brief: A Seq2Seq Model Implementation by keras Recurrent Neural Network LSTM.
# Author: Tateo_YANAGI @soarcloud.com
#
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from keras.models import Sequential
from keras.layers.recurrent import LSTM
from keras.layers.core import Dense, Activation
from keras.layers.wrappers import TimeDistributed
from keras.layers import RepeatVector

PI = 3.1415926535
EPOCHS = 100
TRAIN_SPLIT = 0.98

#
# Class Prediction for lstm model.
#
class Prediction:
  def __init__(self):
    # Array dimension of input&output:(batch_size, time_steps or vectors, features).
    self.batch_size = 48
    self.in_time_steps = 6
    self.in_features = 2
    self.out_vectors = 3
    self.out_features = 2
    self.hidden_neurons = 100 # Neurons of hidden layer.
    self.epochs = EPOCHS # Iteration times.
    self.train_split = TRAIN_SPLIT # Partition of learning dataset.
    self.activation_hidden = 'tanh' # Activation function for hidden units.
    self.activation_output = 'linear' # Activation function for output units.
    self.optimizer= 'adam' # Optimization function.
    self.loss = 'mse' # Loss(cost) function.

#
# Create dataset.
#
  def create_dataset(self, train_split=0.8, in_time_steps=10, out_vectors=4, in_features=2, out_features=2):
    a = np.array([[np.cos(0.0), np.sin(0.0)]])
    for i in range(1, int(2*PI*10)):
      a = np.append(a, np.array([[np.cos(i*0.1), np.sin(i*0.1)]]), axis=0)
    print(f'[Debug-0]\nlen(a):{len(a)}\na:\n{a}')

    # Create dataset for train and test
    # Initialize x_train and y_train.
    x_total = a[0:in_time_steps, 0:in_features]
    y_total = a[in_time_steps:in_time_steps+out_vectors, 0:in_features]

    # Calculate length for x_total and y_total.
    total_len = a.shape[0] - in_time_steps - out_vectors + 1
    print(f'[Debug-1]\ntotal_len:{total_len}')

    # Fill out x_total and y_total.
    for i in range(1, total_len):
      x_total = np.append(x_total, a[i:i+in_time_steps], axis=0)
      y_total = np.append(y_total, a[i+in_time_steps:i+in_time_steps+out_vectors], axis=0)
    print(f'[Debug-2]\nx_total.shape[0]):{x_total.shape[0]}\nx_total:\n{x_total}\ny_total.shape[0]:{y_total.shape[0]}\ny_total:\n{y_total}')

    # Reshape x_toal and y_total from 2D to 3D.
    x_total = np.reshape(x_total[0:x_total.shape[0]], (-1, in_time_steps, in_features))
    y_total = np.reshape(y_total[0:y_total.shape[0]], (-1, out_vectors, out_features))

    # Split dataset for train and test
    x_train = x_total[0:int(x_total.shape[0]*train_split)]
    y_train = y_total[0:int(x_total.shape[0]*train_split)]
    #x_test  = x_total[x_train.shape[0]:]
    #y_test  = y_total[y_train.shape[0]:]
    x_test  = x_total[-2:]
    y_test  = y_total[-2:]
    print(f'[Debug-3]\nx_train.shape[0]:{x_train.shape[0]}\nx_train:\n{x_train}\ny_train.shape[0]:{y_train.shape[0]}\ny_train:\n{y_train}')
    print(f'[Debug-4]\nx_test.shape[0]:{x_test.shape[0]}\nx_test:\n{x_test}\ny_test.shape[0]:{y_test.shape[0]}\ny_test:\n{y_test}')

    return x_train, y_train, x_test, y_test

#
# Create lstm model.
#
  def create_model(self) :
    model = Sequential()
    # Encoder
    model.add(LSTM(self.hidden_neurons, activation=self.activation_hidden, input_shape=(self.in_time_steps, self.in_features)))
    # Output used as Input.
    model.add(RepeatVector(self.out_vectors))
    # Decoder, output sequence
    model.add(LSTM(self.hidden_neurons, activation=self.activation_hidden, return_sequences=True))
    model.add(TimeDistributed(Dense(self.out_features, activation=self.activation_output)))
    #model.add(Activation(self.activation_output))   
    model.compile(optimizer=self.optimizer, loss=self.loss, metrics=['accuracy'])
    return model

  def train_model(self,model, x_train, y_train) :
    model.fit(x_train, y_train, epochs=self.epochs, verbose=2, batch_size=self.batch_size)
    return model

#
# main()
#
if __name__ == "__main__":
  predict = Prediction()
  train_in_data = None
  # Create dataset
  x_train, y_train, x_test, y_test = predict.create_dataset(predict.train_split, predict.in_time_steps,\
                                                            predict.out_vectors, predict.in_features, predict.out_features)
  # Create model
  model = predict.create_model()
  # Train model
  model = predict.train_model(model, x_train, y_train)
  print(model.summary())
  # Test
  predicted = model.predict(x_test, verbose=1)
  print(f'[info_1]predicted:\n{predicted}')
  # Plot result
  predicted = np.reshape(predicted, (-1, predict.out_features))
  x_test = np.reshape(x_test, (-1, predict.in_features))
  y_test = np.reshape(y_test, (-1, predict.out_features))
  x_train = np.reshape(x_train, (-1, predict.in_features))
  y_train = np.reshape(y_train, (-1, predict.out_features))
  plt.figure(figsize=(8, 8))
  plt.scatter(x_train[:,0], x_train[:,1])
  plt.scatter(x_test[:,0], x_test[:,1])
  plt.scatter(y_test[:,0], y_test[:,1])
  plt.scatter(predicted[:,0], predicted[:,1])
  plt.show()

※ソースコードは、
https://github.com/soarbear/lstm_seq2seq_model_predictionへ公開した。

結果

lstm_seq2seq_model_prediction
lstm_seq2seq_model_prediction

今回、epochsとneuronsを増やすことで、上図のとおり精度において納得できる結果を得った。これでコスト関数、活性化関数の選定はじめ、ハイパーパラメータの調整が精度にいかに重要なことを分かる。熱伝播、振動、柔軟ロボット関節の駆動力など元々フーリエ変換、ラプラス変換、Z変換しか解けない偏微分方程式が、LSTMモデルを適用できたらどうかと今度試してみる。

参考文献

1、LSTMネットワークを理解する(英文原稿)、Christopher Olah氏
2、Keras公式資料

1+

ロボット・ドローン部品お探しなら
ロボット翔・電子部品ストア

深層学習から上がる株みつかる

概要

株価、為替、天気、動画など時系列データの予測でよく使われるディープラーニングの代表的手法RNN(再帰型ニューラルネットワーク)の拡張バージョンに、LSTM(Long short-term memory)と呼ばれるモデルがある。今回はLSTM Many to Oneモデルを実装して、複数銘柄(例:10銘柄)の株価から翌日の上がる株を探ってみる。

環境

keras 2.2.5 LSTM
Google Colab CPU/GPU/TPU
Ubuntu 18.04.3 LTS
Python 3.6.8
Numpy 1.17.3
Pandas 0.25.2
sklearn 0.21.3

実装

Many-to-Oneモデルの例として、以下2銘柄から翌日の上がる株を探る。

# -*- coding: utf-8 -*-
import numpy
import pandas
import matplotlib.pyplot as plt

from sklearn import preprocessing
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.recurrent import LSTM

class Predict:
  def __init__(self):
    self.length_of_sequences = 10
    self.in_out_neurons = 1
    self.hidden_neurons = 300
    self.batch_size = 32
    self.epochs = 100
    self.percentage = 0.8

  # データ用意
  def load_data(self, data, n_prev):
    x, y = [], []
    for i in range(len(data) - n_prev):
      x.append(data.iloc[i:(i+n_prev)].values)
      y.append(data.iloc[i+n_prev].values)
    X = numpy.array(x)
    Y = numpy.array(y)
    return X, Y

  # モデル作成
  def create_model(self) :
    Model = Sequential()
    Model.add(LSTM(self.hidden_neurons, batch_input_shape=(None, self.length_of_sequences, self.in_out_neurons), return_sequences=False))
    Model.add(Dense(self.in_out_neurons))
    Model.add(Activation("linear"))
    Model.compile(loss="mape", optimizer="adam")
    return Model

  # 学習
  def train(self, x_train, y_train) :
    Model = self.create_model()
    Model.fit(x_train, y_train, self.batch_size, self.epochs)
    return Model

if __name__ == "__main__":

  predict = Predict()
  nstocks = 2;

  # 銘柄毎に学習、予測、表示
  for istock in range(1, nstocks + 1):
    # データ準備
    data = None
    data = pandas.read_csv('/content/drive/My Drive/LSTM/csv/' + str(istock) + '_stock_price.csv')
    data.columns = ['date', 'open', 'high', 'low', 'close']
    data['date'] = pandas.to_datetime(data['date'], format='%Y-%m-%d')

    # 終値のデータを標準化
    data['close'] = preprocessing.scale(data['close'])
    data = data.sort_values(by='date')
    data = data.reset_index(drop=True)
    data = data.loc[:, ['date', 'close']]

    # 割合で学習、試験データ分割
    split_pos = int(len(data) * predict.percentage)
    x_train, y_train = predict.load_data(data[['close']].iloc[0:split_pos], predict.length_of_sequences)
    x_test, y_test = predict.load_data(data[['close']].iloc[split_pos:], predict.length_of_sequences)

    # 学習
    model = predict.train(x_train, y_train)

    # 試験
    predicted = model.predict(x_test)
    result = pandas.DataFrame(predicted)
    result.columns = [str(istock) + '_predict']
    result[str(istock) + '_actual'] = y_test

    # 表示
    result.plot()
    plt.show()

    # 翌日株価比較
    current = result.iloc[-1][str(istock) + '_actual']
    predictable = result.iloc[-1][str(istock) + '_predict']
    if (predictable - actual) > 0:
      print(f'{istock} stock price of the next day INcreases: {predictable-actual:.2f}, predictable:{predictable:.2f}, current:{current:.2f}')
    else:
      print(f'{istock} stock price of the next day DEcreases: {actual-predictable:.2f}, predictable:{predictable:.2f}, current:{current:.2f}')

またソースコード、csvファイルは、https://github.com/soarbear/stocks-lstm-kerasへ公開した。

結果

keras lstmで複数銘柄の株価予測
keras lstmで複数銘柄の株価予測

【追記】再現性

毎回学習してから予測の結果が変わるので、再現性、生産性がない。ユニット重みや偏りの初期値など【Randomness in Initialization, Regularization, Layers, Optimization】がRandomであり、以下の例として、SEED_IDを固定することによって、予測結果の再現が確認できた。ただしepochsは、SEED_IDを固定しないRandomのほうよりlossがいかに最小まで収束するかを確かめる必要がある。

import numpy as np
import tensorflow as tf
import random as rn
import os
os.environ['PYTHONHASHSEED'] = '0'
from keras import backend as K
session_conf = tf.ConfigProto(intra_op_parallelism_threads=1, inter_op_parallelism_threads=1)
sess = tf.Session(graph=tf.get_default_graph(), config=session_conf)
K.set_session(sess)
......
def fix_seed():
  np.random.seed(SEED_ID)
  rn.seed(SEED_ID)
  tf.set_random_seed(SEED_ID)
......
def create_model(self) :
  Model = Sequential()
  fix_seed()
  Model.add(LSTM(self.hidden_neurons, batch_input_shape=(None, self.length_of_sequences, self.in_out_neurons), return_sequences=False))
  fix_seed()
  Model.add(Dense(self.in_out_neurons))
  Model.add(Activation("linear"))
  Model.compile(loss="mape", optimizer="adam")
  return Model
......

また、2のKeras公式資料もあわせて確認しておく。

【追記】学習状況をスマホに記録

学習状況をスマホでモニタしたい際、アプリHyperdashが使える。以下コードを実行する。学習状況の履歴がスマホに残る。

!pip install hyperdash
from hyperdash import monitor_cell
!hyperdash signup --github
%%monitor_cell 'xxx'

参考文献

1、LSTMネットワークを理解する(英文原稿)、Christopher Olah氏
2、開発中にKerasを用いて再現可能な結果を得るには?、Keras公式資料

3+

ロボット・ドローン部品お探しなら
ロボット翔・電子部品ストア

流体力学ベルヌーイの定理

概要

ベルヌーイの定理(Bernoulli’s theroem)は流体のエネルギ保存則(the conversation law of energy)で、固体のエネルギ保存則との相違点が圧力エネルギが含まれている。以下にまとめてみる。
物体のエネルギ保存則:
$$\frac{1}{2}mu^2+mgz = 一定\space [J] $$
ベルヌーイの定理(流体のエネルギ保存則):
単位体積あたりのエネルギ
$$ \frac{1}{2}ρu^2+ρgz+p = 一定\space [J/m^3 = Pa] $$
単位質量あたりのエネルギ
$$ \frac{1}{2}u^2+gz+\frac{p}{ρ} = 一定\space [J/kg] $$
単位重量あたりのエネルギ
$$ \frac{1}{2g}u^2+z+\frac{p}{ρg} = 一定\space [J/N = m] $$
ただし、\(\frac{1}{2}ρu^2、\frac{1}{2}u^2、\frac{1}{2g}u^2 \) = 運動エネルギ、\( ρgz、gz、z \) = 位置エネルギ、\( p、\frac{p}{ρ}、\frac{p}{ρg} \) = 圧力エネルギ、p = 圧力\([Pa]\)、ρ = 流体密度\([kg/m^3]\)、u = 平均流速\([m/s]\)、g = 重力加速度\([9.8m/s^2]\)

参考書籍

はじめての流体力学、田村恵万氏著、科学図書出版

1+

ロボット・ドローン部品お探しなら
ロボット翔・電子部品ストア

オイラー角~ジンバルロック~クォータニオン

オイラー角

オイラー角とは、三次元ユークリッド空間中の2つの直交座標系(デカルト座標系 )の関係を表現する方法の一つである。ロボット、ドローンの姿勢(ポーズ、向き)を直観に表現するには、オイラー角を用いることがある。オイラー角は3つの角度の組に、3座標軸まわりの回転順序(計12通りもある)で表されるので、オイラー角の変数は4つである。ロボットが姿勢Aから姿勢Bへ変わる際に、例えばx軸-y軸-z軸の順にx-y-z系オイラー角 \( (\alpha,\beta,\gamma) \) で回転することで姿勢の転換を表すことが可能である。姿勢B座標 = 回転行列・姿勢A座標より、姿勢A座標から回転後の姿勢B座標が求められる。姿勢A座標、姿勢B座標が分かれば、回転行列からオイラー角が求められる場合ある。積 \( R_{z}( \gamma ) R_{y}( \beta )R_{x}( \alpha ) \) は、x-y-z系で表したときのオイラー角が \( (\alpha,\beta,\gamma) \) であるような回転を表す回転行列\( R(\alpha,\beta,\gamma) \)である。

また、回転軸がワールド座標軸ローカル座標軸により、オイラー角が以下の2種類ほどある。
・ワールド座標系の3つの軸を中心とした回転、座標軸は静止しているので、静的オイラー角と呼ばれる。
・ローカル座標系の3つの軸を中心とした回転で、回転中に座標軸がロボットとともに回転するため、動的オイラー角と呼ばれる。

ロボットやドローンが移動の際、ワールド座標系を参照できないと、姿勢の転換を表す動的オイラー角を使うしかない場合がある。

ジンバルロック現象、他の問題点

オイラー角の問題点として、動的オイラー角を使用すると、オイラー角そのものの定義より、ジンバルロック現象を起こしてしまう場合がある。但し、静的オイラー角にはジンバルロック現象が起きない。オイラー角の2番目の回転角度が±π/2であれば、1番目と3番目の回転が同じ回転軸になることで、1自由度が失われて、また回転行列より、\( \alpha,\gamma \)の解は無限にあることをジンバルロック現象と呼ばれる。また他にジンバルロック現象を起こさないオイラー角\( (\alpha,\beta,\gamma) \)に対して、\( (\alpha+π,π−\beta,\gamma+π) \)の回転で同じ姿勢になる問題点もある。

よって、以下の条件付きであれば、ジンバルロック現象が除けて、オイラー角でロボットの姿勢を表すのが一意になる。言い換えれば、回転行列がオイラー角と一対一に対応することになる。

・ \( \beta ∈ (-π/2, π/2), \alpha,\gamma ∈ [-π, π] \)

しかし、オイラー角\( (\alpha,\beta,\gamma) \)を制限するのは無理の場合、オイラー角からクォータニオン・四元数の登場となる。

また、オイラー角の2番目の回転軸が1番目の回転軸 x 3番目の回転軸の外積の方向にある。つまり、2番目の回転軸が1番目と3番目の回転軸と直交する。1番目と3番目の回転軸の外積が0になるのが、ジンバルロック現象が発生する場合である。

クォータニオン、四元数

姿勢の表現には、オイラー角の3回転よりも、単純に回転軸と回転軸まわりの1回転で済む場合、クォータニオンが用いられる。クォータニオンの変数の個数がオイラー角と同様に4つであり、クォータニオンを見るだけでは姿勢転換のイメージが難しい。回転行列からオイラー角を求めるのと、逆に回転順とオイラー角から回転行列を求めるのが面倒だが、クォータニオンで回転を表現するとさっぱりとなる。

例えば以下のように、ベクトルxで表すロボットの姿勢Aを原点を通る単位ベクトルu回転軸として、右ねじが進む方向に角θだけ回転させて、ベクトルx’で表すロボットの姿勢Bになると、以下に優雅なる表現がある。$$ x’=q x \bar{q} \\ q=cos \frac{θ}{2} + u sin \frac{θ}{2} \\ \bar{q}=cos \frac{θ}{2} – u sin \frac{θ}{2} $$
但し、\( u=iu_1+ju_2+ku_3 \)

また、回転 \(q_1\)に引き続き、回転 \(q_2\)を行う場合、次のように書くと良い、
$$ x’=(q_{2}q_{1})x(\bar{q_{1}}\bar{q_{2}}) $$

参考文献

・wiki: オイラー角
・wiki: 四元数

関連記事

SLAM~拡張カルマンフィルタ
SLAM~Unscentedカルマンフィルタ
粒子フィルタ
6軸IMU~拡張カルマンフィルタ
研究開発用 台車型ロボット キット

1+

ロボット・ドローン部品お探しなら
ロボット翔・電子部品ストア

ブラシレスDCモータ~Hブリッジ回路

概要

DCブラシレスモータはブラシや整流子に依存しなくなり、代わりに整流用の半導体デバイスを使用する。同期モータの一種で、特性はDCモータと同様。速度はモータ電圧に比例し、トルクはモータ電流に比例する。下図のように、回転部としてのロータは外側にある、アウターロータ型モータと呼ばれるのを例に説明してみる。回転子の内周には磁石が配置されており、U、V、W相コイルの位置は120°ずれている。三相コイルは中心(中性点)で互いに接続されている。U、V、W相コイルの外側にホールセンサが配置され、出力信号はプルアップされて制御プロセッサに入力される。アウターロータ(回転子、磁石)のN極がホールセンサに近いときはH、S極がホールセンサに近いときはLとなる。

三相ブラシレスDCモータの例
三相ブラシレスDCモータの例

Hブリッジ回路でモータ制御

ブラシレスDCモータを駆動するHブリッジ回路の一例は下図のように示される。

BLDCM_H-bridge
BLDCM_H-bridge

下図のように、120°矩形波のU、V、W相ホールセンサと、U、V、W相電圧のH、Lの対応関係が分かる。
BLDCM_Pullup
BLDCM_Pullup

よって、アウターロータが一周360°回転の場合、U、V、W相ホールセンサのUH、UL、VH、VL、WH、WLがそれぞれ180°、U、V、W相電圧のUH、UL、VH、VL、WH、WLがそれぞれ120°と分かる。
ブラシレスDCモータには、FOC(ベクトル周波数変換、磁界ベクトル方向制御とも呼ばれる)、方形波制御(台形波制御、120°制御、6ステップ整流制御とも呼ばれる)、正弦波制御の3つの主な制御方法がある。

矩形波制御

矩形波制御は、ホールセンサまたは無誘導推定アルゴリズムを使用してモータの回転子の位置を取得し、360°の電気サイクルで回転子の位置に応じて6回の転流(60°の転流ごと)を実行する。各転流位置モータは特定の方向に力を出力するので、矩形波制御の位置精度は電気的に60°であると言える。この制御では、モータの相電流波形は方形波に近いため、矩形波制御と呼ばれている。

正弦波制御

正弦波制御方式はSVPWM波を使用し、出力は三相正弦波電圧であり、対応する電流も正弦波電流。矩形波制御と比較してトルク変動が少なく、高調波が少なく、制御時の「細かい」感じが明らかだが、制御器の性能要件は矩形波制御よりわずかに高く、モータ効率が発揮できない。

FOC制御

正弦波制御は、電圧ベクトルの制御を実現し、間接的に電流の大きさの制御を実現するが、電流の方向を制御することはできない。 FOC制御は、電流ベクトルの制御、すなわちモータの固定子磁界のベクトル制御を実現する正弦波制御の改良版と見なすことができる。

1+

ロボット・ドローン部品お探しなら
ロボット翔・電子部品ストア