suzuzusu日記

(´・ω・`)

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

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

線形分離不可能なデータ

以下のように半径が異なる円周上に存在する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)

参考

Hankel Singular Value Decompositionによる時系列データの前処理

HSVD

Hankel Singular Value Decomposition(HSVD)は時系列データを低周波と高周波の2つの成分に分解する手法です.HSVDを使用して時系列データの前処理をします.

Embedding

まず最初に,N個の時系列データ  X(t) を以下のような窓サイズLのハンケル行列H(正確にはハンケル行列は正方行列なのでちょっと違います...)に変換します.

 K = N - L + 1

f:id:suzuzusu:20191107220308p:plain

Decomposition

得られたハンケル行列H特異値分解して低ランク近似をして低周波成分を抽出します.

H = U \Sigma V^{\ast}

低ランク近似したハンケル行列をH'とします.

f:id:suzuzusu:20191107223125p:plain

Unembedding

低周波成分C_Lを以下のように抽出します.

C_L = [ H'(1,1), H'(1, 2), \cdots, H'(1, K), H'(2, K), H'(3, K), \cdots, H'(L, K) ]

以下の画像のように赤い要素が低周波成分になります.

f:id:suzuzusu:20191107222527p:plain

高周波成分C_Hは元のデータから低周波成分を引いたものになります.

C_H(t) = X(t) - C_L(t)

例としてsin波形にノイズを加えたデータをHSVDで前処理を施します.

データ

f:id:suzuzusu:20191107222404p:plain
raw

分解

低周波

f:id:suzuzusu:20191107222419p:plain
low

高周波

f:id:suzuzusu:20191107222435p:plain
high

まとめ

f:id:suzuzusu:20191107222453p:plain
summary

上の画像から元のデータを低周波と高周波に分解できていることが分かると思います. 元のデータ = 低周波 + 高周波 なのでそれぞれの成分を時系列予測アルゴリズムにかけるなどの応用が可能です.

実装

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にあげました

github.com

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で獲得する

データの作成

以下の画像から,しょぼん(´・ω・`)基底を作成します.

f:id:suzuzusu:20191106035235p:plain
オリジナル画像

まず,画像を読み込みます.

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)

プロットすると以下のような,しょぼん(´・ω・`)基底が得られたことが確認できます.

f:id:suzuzusu:20191106035727p:plain
しょぼん(´・ω・`)基底

高次元変換

ガウス分布に従う重みを用いて,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)

結果

得られた潜在変数をプロットすると以下のようになります.

f:id:suzuzusu:20191106040414p:plain
潜在変数

しょぼん(´・ω・`)基底と似たものが獲得できたので成功ですね!!

実装

github.com

参考

Higher-Order-SVDで画像の低ランク近似

HOSVD

SVDは行列の低ランク近似手法であり,Higher-Order-SVD(HOSVD)はSVDをテンソルに拡張した手法です. HOSVDではタッカー分解というコアテンソルとそれぞれのモードの行列に分解する方法を使用します.以下にその概要を図示しました.

f:id:suzuzusu:20191104153601p:plain
タッカー分解

方法

m \times n \times l の3階テンソルを例にしてタッカー分解をする方法を説明します.まず最初に,m \times n \cdot ln \times l \cdot ml \times m \cdot nのそれぞれの行列に分解しSVDをします.

A = U \Sigma V^{\ast}

上記の式から,分解されたユニタリ行列UからコアテンソルSを以下のように定義します.

S :=\{ U_{m}^{T}, U_{n}^{T}, U_{l}^{T} \} \cdot M

元のテンソルに復元するためには,以下のように再度ユニタリ行列Uを使います.

M :=\{ U_{m}, U_{n}, U_{l} \} \cdot S

上記のユニタリ行列Uを任意のランクに切り取ることで,低ランク近似が可能となります.

3階テンソルであるカラー画像を使って,画像の低ランク近似をしてみます.使用する画像は画像処理でよく使われるレナ画像を使用します.

オリジナル画像

f:id:suzuzusu:20191104160723p:plain

Rank1の画像

f:id:suzuzusu:20191104160855p:plain

Rank5の画像

f:id:suzuzusu:20191104160906p:plain

Rank10の画像

f:id:suzuzusu:20191104160914p:plain

Rank20の画像

f:id:suzuzusu:20191104160924p:plain

Rank30の画像

f:id:suzuzusu:20191104160934p:plain

Rank50の画像

f:id:suzuzusu:20191104160942p:plain

Rank100の画像

f:id:suzuzusu:20191104160950p:plain

実装コード

gist.github.com

github.com

参考