suzuzusu日記

(´・ω・`)

RustのHashMapの再現性を確保する

Rustの標準のHashMapのハッシュアルゴリズムでは再現性を確保することができないので,別のハッシュアルゴリズムを使って再現性を確保する方法をメモしておく.これはHashSetの場合も同様に再現性を確保することができます.

HashMapのランダム性

By default, HashMap uses a hashing algorithm selected to provide resistance against HashDoS attacks. The algorithm is randomly seeded, and a reasonable best-effort is made to generate this seed from a high quality, secure source of randomness provided by the host without blocking the program.

doc.rust-lang.org

上記の引用のように,標準のHashMapのハッシュアルゴリズムではHashDoS攻撃を防ぐためにランダム性があり,同じコードでも再現不可能となる場合があります.

再現不可能なコード

例えば以下のようなコードは実行ごとに結果が異なります.

use std::collections::HashMap;

fn main() {
    let mut map = HashMap::new();
    for i in 0..10 {
        map.insert(i, i);
    }
    let _ = map.iter().map(|x| print!("{:?},", x)).collect::<Vec<_>>();
    println!("");
}

実行1

(1, 1),(7, 7),(8, 8),(6, 6),(3, 3),(4, 4),(0, 0),(9, 9),(2, 2),(5, 5),

実行2

(0, 0),(3, 3),(5, 5),(6, 6),(2, 2),(4, 4),(7, 7),(8, 8),(9, 9),(1, 1),

rust-fnv

fnvはkeyのサイズが小さい場合に効率的に動作するハッシュアルゴリズムです. このライブラリには,標準のHashMapをfnvハッシュアルゴリズムに代替してエイリアスされたものがあります.今回はこれを使用して再現性を確保します.

再現可能なコード

use fnv::FnvHashMap;

fn main() {
    let mut map = FnvHashMap::default();
    for i in 0..10 {
        map.insert(i, i);
    }
    let _ = map.iter().map(|x| print!("{:?},", x)).collect::<Vec<_>>();
    println!("");
}

実行1

(5, 5),(4, 4),(7, 7),(6, 6),(1, 1),(0, 0),(3, 3),(2, 2),(9, 9),(8, 8),

実行2

(5, 5),(4, 4),(7, 7),(6, 6),(1, 1),(0, 0),(3, 3),(2, 2),(9, 9),(8, 8),

参考

時系列予測フレームワークfableを試す

fableとは?

https://fable.tidyverts.org/reference/figures/logo.png

ETS,ARIMA,Random Walk,Neural Network autoregressionなどの時系列予測手法群の評価や予測結果の描画をしてくれるフレームワークです. 短いコード量でそれぞれのモデルを手軽に実行できるのが特徴です.使い方を忘れないようにメモしておきます. 作者の中にはRのforecastなどを開発している Rob Hyndman先生もいます.

インストール方法

install.packages("fable")

モデル一覧

現時点では以下のような手法が実装されています.

  • ARIMA
  • ETS(Exponential smoothing state space models)
  • TSLM(Time series linear models)
  • Simple forecasting methods(移動平均モデル,ランダムウォークなど)
  • Neural network autoregression
  • Vector autoregression

Function reference • fable

使い方

データフレームのtsibbleに変換してmodelに突っ込んでforecastするだけです. modelに予測したいモデルを複数入れることが可能です.

tsibble_data %>%
  model(
    # value は予測したいデータ
    hoge = hoge_model(value),
    fuga = fuga_model(log(value)) # 対数変換やBox-Cox変換も可能
  ) %>%
  forecast()

AirPassengersの時系列予測

f:id:suzuzusu:20191215234315p:plain
AirPassengers

おなじみのAirPassengersで試してみます.以下が,前処理も含めた12種類ほどの手法で比較した実装です.

library(fable)
library(tsibble)
library(lubridate)
library(dplyr)

AirPassengers %>%
  as_tsibble %>%
  model(
    ets = ETS(value),
    ets_log = ETS(log(value)),
    ets_box_cox = ETS(box_cox(value, 0.3)),
    arima = ARIMA(value),
    arima_log = ARIMA(log(value)),
    arima_box_cox = ARIMA(box_cox(value, 0.3)),
    snaive = SNAIVE(value),
    snaive_log = SNAIVE(log(value)),
    snaive_box_cox = SNAIVE(box_cox(value, 0.3)),
    nnetar = NNETAR(value),
    nnetar_log = NNETAR(log(value)),
    nnetar_box_cox = NNETAR(box_cox(value, 0.3))
  ) %>%
  forecast(h = "2 years") %>% 
  autoplot(filter(as_tsibble(AirPassengers), year(index) > 1950), level = NULL)

f:id:suzuzusu:20191215234429p:plain
結果

f:id:suzuzusu:20191215234441p:plain
結果(拡大)

参考

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なモデルの方が便利です.

入力は以下のような三角波です.

f:id:suzuzusu:20191213011127p:plain

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()

f:id:suzuzusu:20191213011913p:plain

Colab

gist.github.com

参考

多項式カーネルに対応する特徴ベクトルを可視化する

多項式カーネルによって写像される特徴空間を可視化してカーネル法の挙動を把握します.

線形分離不可能なデータ

以下のように半径が異なる円周上に存在する2クラスのデータを生成します.

f:id:suzuzusu:20191204190233p:plain

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クラスのデータは直線を描いて分離することができない線形分離不可能なデータです.このようなデータを識別するためには直線を描いて分離できるような空間に変換したほうが都合がよいです.次のカーネル法でその問題を解決します.

カーネル法

先程のデータを線形分離可能な特徴空間に写像するための写像関数を\phi(\boldsymbol{x})とします.カーネル法では明示的に\phi(\boldsymbol{x})を定義せずに特徴空間の内積カーネル関数 K(\boldsymbol{x}, \boldsymbol{y})=\phi(\boldsymbol{x}) ^T \phi(\boldsymbol{y})によって定義することで写像関数\phi(\boldsymbol{x})の特徴空間で分離して識別することが可能となる手法の総称です.今回はそのカーネル関数によってどのような特徴空間に写像されて識別しているのかを可視化します.

多項式カーネル

(式が揃っていないのは許して...分からなかった...)

多項式カーネルはパラメータc, dによって次の定義されるカーネル関数です.

 K(\boldsymbol{x}, \boldsymbol{y})=(\boldsymbol{x}^T \boldsymbol{y} + c)^d

\boldsymbol{x} = [ x_1, x_2 ]として,c=0,d=1のときの特徴ベクトルは次のようになります.


\begin{eqnarray}
K(\boldsymbol{x}, \boldsymbol{y}) = (\boldsymbol{x}^T \boldsymbol{y})^1 \\
= x_1 y_1 + x_2 y_2 \\
= [ x_1, x_2 ] ^T [ y_1, y_2 ] \\
= \phi(\boldsymbol{x}) ^T \phi(\boldsymbol{y})
\end{eqnarray}

よって特徴ベクトル\phi(\boldsymbol{x}) = [ x_1, x_2 ]となる.

同様に,c=0,d=2のときの特徴ベクトルは


\begin{eqnarray}
K(\boldsymbol{x}, \boldsymbol{y}) = (\boldsymbol{x}^T \boldsymbol{y})^2 \\
= x_1 ^2 y_1 ^2   + 2 x_1 x_2 y_1 y_2 + x_2 ^2 y_2 ^2 \\
= [ x_1 ^2  + \sqrt{2} x_1 x_2 + x_2 ^2 ] ^T [ y_1 ^2  + \sqrt{2} y_1 y_2 + y_2 ^2 ] \\
= \phi(\boldsymbol{x}) ^T \phi(\boldsymbol{y})
\end{eqnarray}

よって特徴ベクトル\phi(\boldsymbol{x}) = [ x_1 ^2 , \sqrt{2} x_1 x_2, x_2 ^2 ]となる.

 dが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としてプロットしてみます.

f:id:suzuzusu:20191204193527p:plain

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(次数)が偶数の場合

f:id:suzuzusu:20191204195346p:plain f:id:suzuzusu:20191204195349p:plain f:id:suzuzusu:20191204195354p:plain f:id:suzuzusu:20191204195359p:plain f:id:suzuzusu:20191204195407p:plain

偶数の場合は一部(右上,左下)の特徴空間を見ると線形分離可能な空間に写像されていることがよくわかります.

d(次数)が奇数の場合

f:id:suzuzusu:20191204195909p:plain f:id:suzuzusu:20191204195912p:plain f:id:suzuzusu:20191204195915p:plain f:id:suzuzusu:20191204195921p:plain f:id:suzuzusu:20191204195900p:plain

奇数の場合はどの写像した空間を見ても線形分離できないことが確認できます.

参考

幾何ブラウン運動をARIMAで予測

忘備録としてまとめておく。

幾何ブラウン運動

以下のような確率微分方程式(SDE)によって定義されるものを幾何ブラウン運動と言います。金融工学のモデルとしてよく用いられます。

d S_t = \mu S_t dt + \sigma S_t dB_t

初期値をS_0とすると解は次のようになることが知られています。

 S_t = S_0 \exp( (\mu - \frac{\sigma ^2}{2})t + \sigma B_t) )

シミュレーションすると次のようになります。

# 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")

f:id:suzuzusu:20191124020747p:plain
幾何ブラウン運動

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"))

f:id:suzuzusu:20191124021654p:plain

こちらが拡大図です。

# 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"))

f:id:suzuzusu:20191124021711p:plain
拡大図

実装まとめ

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分布が外れ値に頑健な回帰をすることが可能であることを確認します.

線形回帰

今回は以下のパラメータ\alpha,\betaの単純な線形回帰モデルを使用します.

f(x)=\alpha+\beta x

近似させた \widehat{y} \approx f(x) に対する尤度関数を教師データy_tとしてガウス分布の場合は\mathcal{N}(y_t|\widehat{y}, \sigma),スチューデントのt分布の場合は\mathcal{St}(y_t|\widehat{y}, \nu, \sigma)とします.

データ

求めたい関数を以下に定義して,ノイズと外れ値を含むデータを作成します.

f(x)=100 + 2x

f:id:suzuzusu:20191117151710p:plain
サンプルデータ

上記の画像のようなデータを使用します.

結果

RStanを使用してMCMCサンプリングをしてパラメータ\alpha,\betaの事後分布の平均を使用すると以下のような回帰曲線が得られます.

y_{gaussian} = 295.9909 + 1.904833 x

y_{st} = 103.6868 + 1.985594 x

f:id:suzuzusu:20191117151736p:plain

ガウス分布では外れ値に引っ張られた回帰をしているのに対して,スチューデントの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波形にノイズを付加したデータを使用します.

データ

f:id:suzuzusu:20191108033743p:plain
raw

分解

低周波

f:id:suzuzusu:20191108033803p:plain
low

高周波

f:id:suzuzusu:20191108033816p:plain
high

まとめ

f:id:suzuzusu:20191108033830p:plain
summary

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にあげておきました.

github.com

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)

参考