kerasのLSTMのstatelessとstatefulの切り替え
kerasのLSTMなどのRNN系のモデルは状態を保持するstatefulなモデルと状態を保持しないstatelessなモデルがあります.その切り替え方法をメモしておきます.
切り替え方法
あらかじめ同じネットワークのstatelessなモデルとstatefulなモデルを別々に作成する.
from keras.models import Sequential from keras.layers import Dense, LSTM stateless_model = Sequential() stateless_model.add(LSTM(hidden_unit, input_shape=(seq_len, 1))) stateless_model.add(Dense(1)) stateful_model = Sequential() stateful_model.add(LSTM(hidden_unit, batch_input_shape=(1, 1, 1), stateful=True)) stateful_model.add(Dense(1))
以下のようにネットワークが同じなので重みの共有が可能となる.これによって切り替えることができる.
# stateless to stateful stateful_model.load_weights(stateless_model.get_weights()) # hdf5に保存して重みの共有も可能 stateful_model.load_weights('model.hdf5')
Example
statelessなモデルでトレーニングをしてstatefulなモデルでテストをするコードを以下に示します.再帰的に予測したいときはstatefulなモデルの方が便利です.
入力は以下のような三角波です.
import numpy as np import random from keras.models import Sequential from keras.layers import Dense, LSTM from keras.callbacks import ModelCheckpoint import matplotlib.pyplot as plt def hankel_matrix(x,seq_len): n = x.shape[0] stride = x.strides[0] return np.lib.stride_tricks.as_strided(x, shape=(n-seq_len+1, seq_len), strides=(stride,stride)).copy() def triangle_wave(t, m=10000): arange = np.arange(1, m+1) return 8/(np.pi**2)*np.sum(np.sin(arange*np.pi/2)*np.sin(arange*t)/(arange**2)) seed = 0 hidden_unit = 100 random.seed(seed) np.random.seed(seed) # dataset N = 500 x_n = [] ts = np.linspace(0, 500, N) for t in ts: tmp = triangle_wave(t) x_n.append(tmp) x_n = np.asarray(x_n) seq_len = 100 H = hankel_matrix(x_n, seq_len=seq_len) X = H[:-1,:].reshape(-1, seq_len, 1) Y = H[1:,-1] # plot data plt.plot(x_n[:50]) plt.savefig('data.png') plt.show() cp_cb = ModelCheckpoint(filepath = 'model.hdf5', monitor='val_loss', verbose=1, save_best_only=True, mode='auto') # stateless train_model = Sequential() train_model.add(LSTM(hidden_unit, input_shape=(seq_len, 1))) train_model.add(Dense(1)) train_model.compile(loss='mean_squared_error', optimizer='adam') train_model.fit(X, Y, epochs=20, callbacks = [cp_cb], validation_split=0.3, shuffle=True) # stateful predict_model = Sequential() predict_model.add(LSTM(hidden_unit, batch_input_shape=(1, 1, 1), stateful=True)) predict_model.add(Dense(1)) predict_model.load_weights('model.hdf5') input_len = 30 x_n = [] xs = [] ts = np.linspace(500, 1000, N) for t in ts: tmp = triangle_wave(t) x_n.append(tmp) xs.append(tmp) xs = xs[:input_len] # test input for o in xs: p = predict_model.predict(np.asarray([o]).reshape(1, 1, 1))[0,0] # predict predict_len = 100 xs.append(p) for i in range(predict_len-1): p = predict_model.predict(np.asarray([p]).reshape(1, 1, 1))[0,0] xs.append(p) # plot predict plt.plot(xs[:input_len], label='input value') plt.plot(np.arange(predict_len) + input_len, x_n[input_len:][:predict_len], label='true value') plt.plot(np.arange(predict_len) + input_len, xs[input_len:], label='lstm predict value') plt.xlabel('time') plt.legend(loc='upper right') plt.savefig('fig.png') plt.show()
Colab
参考
多項式カーネルに対応する特徴ベクトルを可視化する
多項式カーネルによって写像される特徴空間を可視化してカーネル法の挙動を把握します.
線形分離不可能なデータ
以下のように半径が異なる円周上に存在する2クラスのデータを生成します.
import seaborn as sns import pandas as pd import matplotlib.pyplot as plt import numpy as np def toy_dataset(n1, n2): # class 1 data set rs = np.random.random(n1) r = 1.0 X1 = np.c_[r * np.cos(2*np.pi*rs), r * np.sin(2*np.pi*rs)] Y1 = ['class1' for _ in range(n1)] # class 2 data set rs = np.random.random(n2) r = 2.0 X2 = np.c_[r * np.cos(2*np.pi*rs), r * np.sin(2*np.pi*rs)] Y2 = ['class2' for _ in range(n2)] # concat X = np.r_[X1, X2] Y = np.r_[Y1, Y2] x_df = pd.DataFrame(data = X, columns = ['x', 'y']) y_series = pd.Series(Y, name = 'class') df = pd.concat([x_df, y_series], axis=1) return df # seed np.random.seed(0) # dataset num = 100 df = toy_dataset(n1=num, n2=num) # dataset scatter c1_df = df[df['class'] == 'class1'] c2_df = df[df['class'] == 'class2'] sns.scatterplot(c1_df['x'], c1_df['y'], label='class1') sns.scatterplot(c2_df['x'], c2_df['y'], label='class2') plt.show()
上のような2クラスのデータは直線を描いて分離することができない線形分離不可能なデータです.このようなデータを識別するためには直線を描いて分離できるような空間に変換したほうが都合がよいです.次のカーネル法でその問題を解決します.
カーネル法
先程のデータを線形分離可能な特徴空間に写像するための写像関数をとします.カーネル法では明示的にを定義せずに特徴空間の内積をカーネル関数によって定義することで写像関数の特徴空間で分離して識別することが可能となる手法の総称です.今回はそのカーネル関数によってどのような特徴空間に写像されて識別しているのかを可視化します.
多項式カーネル
(式が揃っていないのは許して...分からなかった...)
多項式カーネルはパラメータによって次の定義されるカーネル関数です.
]として,のときの特徴ベクトルは次のようになります.
よって特徴ベクトル ]となる.
同様に,のときの特徴ベクトルは
よって特徴ベクトル ]となる.
が3以上の場合は同じように計算すればいいので省略します.
SVMによる識別
さきほどの線形分離不可能なデータをSVMによって識別してみます.その際に多項式カーネルのd(次数)を変えて識別性能を比較します.
トレーニング,テストデータの準備
from sklearn.svm import SVC from sklearn.metrics import accuracy_score np.random.seed(0) # train dataset train_num = 100 df = toy_dataset(n1=train_num, n2=train_num) # test dataset test_num = 30 test_df = toy_dataset(n1=test_num, n2=test_num)
d(次数)が偶数の場合
print('even degree') for d in range(2, 11, 2): clf = SVC(gamma = 'auto', kernel='poly', degree=d) clf.fit(df[['x', 'y']], df['class']) y_pred = clf.predict(test_df[['x', 'y']]) print('degree =', d, 'acc:', accuracy_score(test_df['class'], y_pred))
出力
even degree degree = 2 acc: 1.0 degree = 4 acc: 1.0 degree = 6 acc: 1.0 degree = 8 acc: 1.0 degree = 10 acc: 1.0
d(次数)が奇数の場合
print('odd degree') for d in range(1, 10, 2): clf = SVC(gamma = 'auto', kernel='poly', degree=d) clf.fit(df[['x', 'y']], df['class']) y_pred = clf.predict(test_df[['x', 'y']]) print('degree =', d, 'acc:', accuracy_score(test_df['class'], y_pred))
出力
odd degree degree = 1 acc: 0.5166666666666667 degree = 3 acc: 0.6833333333333333 degree = 5 acc: 0.6833333333333333 degree = 7 acc: 0.6666666666666666 degree = 9 acc: 0.7
まとめ
横軸をd(次数),縦軸をAccuracyとしてプロットしてみます.
dが偶数の場合は識別可能に対して,奇数の場合は識別不可能であることが分かります.なぜこのような結果になったのかを次に特徴ベクトルを可視化して調べてみます.
特徴ベクトルの可視化
特徴ベクトル
多項式カーネルの特徴ベクトルは次のように二項定理から計算することが可能です.
def poly_feature(x, d): n = x.shape[0] Z = np.zeros((n, d+1)) for i in range(d+1): # 二項定理 a = np.sqrt(comb(d, i, exact=True)) Z[:,i] = a * (x[:,0]**(d-i)) * (x[:,1]**(i)) return Z
これによって計算できた特徴ベクトルをseabornのペアプロット図を使って可視化してみます.
from scipy.special import comb from itertools import combinations for d in range(1, 11): z = poly_feature(df[['x', 'y']].values, d) columns = [ 'feature' + str(i) for i in range(d+1)] feature_df = pd.DataFrame(data = z, columns = columns) feature_df = pd.concat([feature_df, df['class']], axis=1) g = sns.pairplot(feature_df, hue='class', vars=feature_df.columns[:-1]) g.fig.suptitle('degree = ' + str(d)) plt.show()
d(次数)が偶数の場合
偶数の場合は一部(右上,左下)の特徴空間を見ると線形分離可能な空間に写像されていることがよくわかります.
d(次数)が奇数の場合
奇数の場合はどの写像した空間を見ても線形分離できないことが確認できます.
参考
幾何ブラウン運動をARIMAで予測
忘備録としてまとめておく。
幾何ブラウン運動
以下のような確率微分方程式(SDE)によって定義されるものを幾何ブラウン運動と言います。金融工学のモデルとしてよく用いられます。
初期値をとすると解は次のようになることが知られています。
シミュレーションすると次のようになります。
# geometric Brownian motion T = 1.0 N = 10000 dt = T/N sigma = 2 mu = 1 X0 = 1 t = seq(from = 0.0, to = T, length.out = N) dw = sqrt(dt) * rnorm(N) dw[1] = 0.0 W = cumsum(dw) X = X0 * exp((sigma - 0.5 * mu^2) * t + (mu*W)) plot(t, X, type="l")
ARIMA
ARIMAはRのforecastパッケージの auto.arima
を使います。
library(forecast) # forecast step h = 1000 # ARIMA model <- auto.arima(X[1:(N-h)]) summary(model)
ARIMA(2, 1, 2)が選択されました。以下がsummaryの出力です。
ARIMA(2,1,2) Coefficients: ar1 ar2 ma1 ma2 0.7916 -0.7906 -0.8089 0.8230 s.e. 0.0604 0.0804 0.0559 0.0746 sigma^2 estimated as 0.001461: log likelihood=16609.56 AIC=-33209.12 AICc=-33209.11 BIC=-33173.59 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set 0.0005448627 0.03820814 0.02565005 0.0147103 0.7971145 0.9991512 ACF1 Training set 0.0009177808
結果
予測の結果を次に示します。
# small plot plot(forecast(model, h=h), xlim=c(0,N), ylim=c(0,15)) par(new=T) plot(c((N-h):N), X[(N-h):N], xlim=c(0,N), ylim=c(0,15), ann=F, type='l', col = "red") legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))
こちらが拡大図です。
# large plot plot(forecast(model, h=h), xlim=c(N-h-200,N), ylim=c(0,15)) par(new=T) plot(c((N-h):N), X[(N-h):N], xlim=c(N-h-200,N), ylim=c(0,15), ann=F, type='l', col = "red") legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))
実装まとめ
library(forecast) set.seed(0) # geometric Brownian motion T = 1.0 N = 10000 dt = T/N sigma = 2 mu = 1 X0 = 1 t = seq(from = 0.0, to = T, length.out = N) dw = sqrt(dt) * rnorm(N) dw[1] = 0.0 W = cumsum(dw) X = X0 * exp((sigma - 0.5 * mu^2) * t + (mu*W)) plot(t, X, type="l") # forecast step h = 1000 # ARIMA model <- auto.arima(X[1:(N-h)]) # small plot plot(forecast(model, h=h), xlim=c(0,N), ylim=c(0,15)) par(new=T) plot(c((N-h):N), X[(N-h):N], xlim=c(0,N), ylim=c(0,15), ann=F, type='l', col = "red") legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue")) # large plot plot(forecast(model, h=h), xlim=c(N-h-200,N), ylim=c(0,15)) par(new=T) plot(c((N-h):N), X[(N-h):N], xlim=c(N-h-200,N), ylim=c(0,15), ann=F, type='l', col = "red") legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))
参考
尤度関数におけるガウス分布とスチューデントのt分布の比較
RStanを使って尤度関数がガウス分布と場合とスチューデントのt分布の場合でどのように異なるのかを比較します. 最終的にガウス分布に比べてスチューデントのt分布が外れ値に頑健な回帰をすることが可能であることを確認します.
線形回帰
今回は以下のパラメータの単純な線形回帰モデルを使用します.
近似させた に対する尤度関数を教師データとしてガウス分布の場合は,スチューデントのt分布の場合はとします.
データ
求めたい関数を以下に定義して,ノイズと外れ値を含むデータを作成します.
上記の画像のようなデータを使用します.
結果
RStanを使用してMCMCサンプリングをしてパラメータの事後分布の平均を使用すると以下のような回帰曲線が得られます.
ガウス分布では外れ値に引っ張られた回帰をしているのに対して,スチューデントのt分布の方が頑健な回帰をしていることが分かります.
実装
library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) library(ggplot2) # データ作成 set.seed(1) N <- 100 x = sort(sample(1:1000, N)) alpha = 100 beta = 2 y = alpha + beta*x y = y + rnorm(N, sd=30) ris = sample(1:100, 30) for (ri in ris) { # 外れ値 y[ri] <- y[ri] + runif(1, 0, 1000) } dat <- list(N = N, x = x, y = y) # 尤度関数がスチューデントのt分布 fit_t <- stan(file = 'student_t.stan', data = dat) # 尤度関数がガウス分布 fit_g <- stan(file = 'gaussian.stan', data = dat) x_tst = c(1:1000) # パラメータの平均値を取得 m <- get_posterior_mean(fit_g) alpha_g = m[1] beta_g = m[2] sigma_g = m[3] y_mean_g = beta_g*x_tst + alpha_g m <- get_posterior_mean(fit_t) alpha_t = m[1] beta_t = m[2] sigma_t = m[3] nu_t = m[4] y_mean_t = beta_t*x_tst + alpha_t p <- data.frame(x = x, y = y, x_tst = x_tst, y_mean_g = y_mean_g, y_mean_t = y_mean_t) g <- ggplot(data=p) + geom_point(aes(x=x, y=y, colour="")) + labs(colour="sample") plot(g) ggsave(file = "sample.png", plot = g) g <- ggplot(data=p) + geom_point(aes(x=x, y=y)) + geom_line(aes(x=x_tst, y=y_mean_g, colour="red")) + geom_line(aes(x=x_tst, y=y_mean_t, colour="darkblue")) + scale_color_discrete(name = "likelihood", labels = c("Student t", "Gaussian")) plot(g) ggsave(file = "summary.png", plot = g)
gaussian.stan
data { int<lower=0> N; vector[N] x; vector[N] y; } parameters { real alpha; real beta; real<lower=0> sigma; } model { y ~ normal(alpha + beta * x, sigma); }
student_t.stan
data { int<lower=0> N; vector[N] x; vector[N] y; } parameters { real alpha; real beta; real nu; real<lower=0> sigma; } model { y ~ student_t(nu, alpha + beta * x, sigma); }
参考
Hankel Non-negative Matrix Factorizationによる時系列データの前処理
HNMF
先日ブログで時系列データをHankel Singular Value Decompositionという低周波成分と高周波成分に分解するという前処理のアルゴリズムを紹介しました. 今回は特異値分解(SVD)の低ランク近似の部分を非負値行列因子分解(NMF)に変えても分解できるのかを試しました. とりあえずHankel Non-negative Matrix Factorization(HNMF)と呼びます.
詳しいアルゴリズムは先日のブログを見てください. 差分はDecompositionのSVDをNMFに変えただけです.
例
NMFを使用するので,時系列データは非負値であるという制約があります. 以下のようなsin波形にノイズを付加したデータを使用します.
データ
分解
低周波
高周波
まとめ
HSVDと同様に低周波成分と高周波成分に分解できたことが確認できました.
実装
import numpy as np from sklearn.decomposition import NMF def moving_window_matrix(x,window_size): # Fork from https://qiita.com/bauer/items/48ef4a57ff77b45244b6 n = x.shape[0] stride = x.strides[0] return np.lib.stride_tricks.as_strided(x, shape=(n-window_size+1, window_size), strides=(stride,stride) ).copy() def hnmf(x, window, rank, random_state = 0): assert(np.min(x) > 0.0) m = moving_window_matrix(x, window) model = NMF(n_components=rank, init='random', random_state=random_state) w = model.fit_transform(m) h = model.components_ hankel = w @ h c = hankel[0,:] c = np.append(c, hankel[1:,-1]) return c, x-c
githubにあげておきました.
Usage
import numpy as np from tfilter import hnmf N = 500 x = np.sin(np.arange(N) * np.pi/50.0) x = x + np.random.normal(0, 0.3, size=N) x = x + 2.0 assert(np.min(x) > 0.0) window = 100 rank = 3 low_freq, high_freq = hnmf(x, window, rank)
参考
Hankel Singular Value Decompositionによる時系列データの前処理
HSVD
Hankel Singular Value Decomposition(HSVD)は時系列データを低周波と高周波の2つの成分に分解する手法です.HSVDを使用して時系列データの前処理をします.
Embedding
まず最初に,個の時系列データ を以下のような窓サイズのハンケル行列H(正確にはハンケル行列は正方行列なのでちょっと違います...)に変換します.
Decomposition
得られたハンケル行列を特異値分解して低ランク近似をして低周波成分を抽出します.
低ランク近似したハンケル行列をとします.
Unembedding
低周波成分を以下のように抽出します.
]
以下の画像のように赤い要素が低周波成分になります.
高周波成分は元のデータから低周波成分を引いたものになります.
例
例としてsin波形にノイズを加えたデータをHSVDで前処理を施します.
データ
分解
低周波
高周波
まとめ
上の画像から元のデータを低周波と高周波に分解できていることが分かると思います. 元のデータ = 低周波 + 高周波 なのでそれぞれの成分を時系列予測アルゴリズムにかけるなどの応用が可能です.
実装
def moving_window_matrix(x,window_size): # Fork from https://qiita.com/bauer/items/48ef4a57ff77b45244b6 n = x.shape[0] stride = x.strides[0] return np.lib.stride_tricks.as_strided(x, shape=(n-window_size+1, window_size), strides=(stride,stride) ).copy() def hsvd(x, window, rank): m = moving_window_matrix(x, window) u, s, vh = np.linalg.svd(m) h = u[:,:rank] @ np.diag(s[:rank]) @ vh[:rank,:] c = h[0,:] c = np.append(c, h[1:,-1]) return c, x-c
実装をgithubにあげました
Usage
import numpy as np from tfilter import hsvd N = 500 x = np.sin(np.arange(N) * np.pi/50.0) x = x + np.random.normal(0, 0.3, size=N) window = 100 rank = 2 low_freq, high_freq = hsvd(x, window, rank)
参考
しょぼん(´・ω・`)基底をGPLVMで獲得する
データの作成
以下の画像から,しょぼん(´・ω・`)基底を作成します.
まず,画像を読み込みます.
import numpy as np from PIL import Image image = np.asarray(Image.open('./data/shobon.png').convert('L'))
次に,2値化を行いしょぼん(´・ω・`)基底を作成します.
# binarization a = np.where(image < 240) shobon = np.c_[a[1], np.flipud(a[0])] # shobon.shape = (896, 2)
プロットすると以下のような,しょぼん(´・ω・`)基底が得られたことが確認できます.
高次元変換
ガウス分布に従う重みを用いて,100次元の高次元にしょぼん(´・ω・`)基底を変換します. そしてノイズを加えておきます.
dim = 100 w = np.random.normal(size=2*dim).reshape(2, dim) high = shobon @ w high = high + np.random.normal(size=high.shape[0]*high.shape[1]).reshape(high.shape[0], high.shape[1]) # high.shape = (896, 100)
GPLVMによる低次元埋め込み
Gaussian Process Latent Variable Models(GPLVM)とはガウス過程を利用して潜在変数空間に変換する教師なし学習のアルゴリズムです. 応用例として以下の動画のように人間のポーズを低次元埋め込みするなどの試みがあります.
今回はGPLVMを使うためにGPyを使用します.
GPLVMを使って100次元のデータを2次元の低次元埋め込みをします.
# latent gplvm = GPLVM(high, 2) latent = gplvm.X # latent.shape = (896, 2)
結果
得られた潜在変数をプロットすると以下のようになります.
しょぼん(´・ω・`)基底と似たものが獲得できたので成功ですね!!