2020-08-04

スチューデントのt分布の最尤推定

技術系機械学習系

スチューデントのt分布は外れ値のあるデータに対してもロバストな推定を可能とする分布として知られている。忘備録としてこの記事では、そのt分布の最尤推定をする。

スチューデントのt分布

一次元のt分布は、パラメータ $\mu, \lambda, \nu$ を用いて次のような確率密度関数で表すことができる。

f:id:suzuzusu:20200803061703p:plain

EMアルゴリズム

最尤推定には、EMアルゴリズムを用いる。 $X$ をデータ、 $Z$ を潜在変数、 $\Theta$ をパラメータとする。次式で尤度関数と対数尤度関数を示す。

f:id:suzuzusu:20200803062219p:plain

f:id:suzuzusu:20200803063006p:plain

Eステップ

$Z$ の事後分布を計算する。

f:id:suzuzusu:20200804072030p:plain

$Z$ の事後分布による $\eta, \log \eta$ の期待値は以下のようになる。

f:id:suzuzusu:20200804073114p:plain

f:id:suzuzusu:20200804073050p:plain

Mステップ

対数尤度関数 $Q(\Theta, \Theta_0)$ の期待値を最大化するように、それぞれのパラメータを更新する。

f:id:suzuzusu:20200804073748p:plain

f:id:suzuzusu:20200804074107p:plain

f:id:suzuzusu:20200804074343p:plain

f:id:suzuzusu:20200804074804p:plain

ただし、パラメータ $\nu$ に関しては、上の式のように解析解が得られないのでニュートン法による数値解析によってパラメータを更新する。

f:id:suzuzusu:20200804075200p:plain

f:id:suzuzusu:20200804075429p:plain

このようにEステップ、Mステップを繰り返すことによって尤度関数を最大化させて最尤推定をする。

実装

juliaを用いて実装をする。外れ値による影響を正規分布と比較するために、 $\mathcal{N}(0.0, 1.0)$ から10のデータのサンプルに $\mathcal{N}(100.0, 1.0)$ から生起されるデータを1つ加えたデータセットでそれぞれを最尤推定して比較する。

using Plots
using StatsPlots
using Random
using Distributions
using SpecialFunctions
using Statistics
using StatsBase
using Optim
using ReverseDiff: gradient

Random.seed!(1234)

function e_eta(x, nu, lambda_, mu)
    return (nu + 1)/(nu + lambda_ * (x - mu)^2)
end

function e_etas(X, nu, lambda_, mu)
    return [ e_eta(x, nu, lambda_, mu) for x in X]
end

function e_log_eta(x, nu, lambda_, mu)
    return digamma((nu + 1.0)/2.0) - log((nu + lambda_ * (x - mu)^2)/2.0)
end

function e_log_etas(X, nu, lambda_, mu)
    return [ e_log_eta(x, nu, lambda_, mu) for x in X]
end

function fit_t(X)
    n = length(X)
    
    # init
    mu = median(X)
    lambda_ = iqr(X)/2.0
    nu = 1.0
   
    for i in 1:1000
        # E step
        e_etas_ = e_etas(X, nu, lambda_, mu)
        e_log_etas_ = e_log_etas(X, nu[1], lambda_, mu)
        
        
        # M step
        # mu
        mu = (X' * e_etas_) / sum(e_etas_)
        # lambda_
        lambda_ = 1.0 / ( (((X .- mu).^2)' * e_etas_) / n)
        # nu
        function f(nu)
            return digamma(nu[1]/2.0) - log(nu[1]/2.0) - (1.0 + sum(e_log_etas_)/n - sum(e_etas_)/n)
        end
        for j in 1:1000
            tmp = nu - f([nu])/gradient(f, [nu])[1]
            if !isnan(tmp)
                nu = max(tmp, 1e-6)
                if abs(f([nu])) < 1e-6
                    break
                end
            else
                break
            end
        end
    end
    
    return mu, lambda_, nu
end

d = Normal(0.0, 1.0)
noise_d = Normal(100.0, 1.0)
data = rand(d, 10)
outlier = rand(noise_d, 1)
X = vcat(data, outlier)
mu, lambda_, nu = fit_t(X)

scatter(data, zeros(length(data)), label="data")
scatter!(outlier, zeros(length(outlier)), label="outlier")
plot!(fit_mle(Normal, X), label="Normal")

Xs = Array(range(-100, 120, step=0.1))
function t_pdf(x, mu, lamda_, nu)
    return gamma((nu+1)/2)./gamma(nu/2).*sqrt(lambda_/(pi*nu)).*(1 .+ lambda_*(x .- mu).^2/nu).^(-(nu+1)/2)
end
Ys = t_pdf(Xs, mu, lambda_, nu)
plot!(Xs, Ys.*0.035, label="Student t Distribution", title="MLE of the Student t distribution using EM Algorithm")

f:id:suzuzusu:20200804220738p:plain

上記の結果のように正規分布は、外れ値に引っ張られた推定をしていることが分かる。一方、t分布は外れ値に対してロバストに推定可能であることが分かる。

notebook

Jupyter Notebook Viewer

参考

2020-03-04

AV1で採用されているChroma from Luma Prediction (CfL)を使ってイントラ予測

技術系機械学習系

Chroma from Luma Prediction (CfL) はAlliance for Open Mediaが開発したオープンかつロイヤリティフリーな動画圧縮コーデックAV1などで採用されているイントラ予測手法です。イントラ予測というのは動画圧縮で他のフレームを参照しないフレーム符号化のことを言います。今回は、Juliaを用いてその効果を検証していきます。

Chroma（彩度）とLuma（輝度）の関係

CfLはLumaからChromaを推定する手法です。これにより保持すべきデータ容量が削減できます。ここでは、レナの画像を使ってChromaとLumaの関係を可視化します。

まず必要なパッケージを使ってレナの画像を準備します。

using Images, ImageView, TestImages
using Plots

lena = testimage("lena_color_256")
plot(lena)

f:id:suzuzusu:20200304055448p:plain

次にYCbCr色空間に変換します。

lena = YCbCr.(lena)
lena_arr = channelview(lena)

p_y = heatmap(reverse(lena_arr[1,:,:], dims=1), title="Y(Luma)", color=:grays)
p_cb = heatmap(reverse(lena_arr[2,:,:], dims=1), title="Cb", color=:grays)
p_cr = heatmap(reverse(lena_arr[3,:,:], dims=1), title="Cr", color=:grays)
plot(p_y, p_cb, p_cr, layout=(1, 3), size=(1200, 300), fmt=:png)

f:id:suzuzusu:20200304055512p:plain

左上の32x32のタイル上における、ビットごとのそれぞれChromaとLumaの関係を見ます。横軸にLuma、縦軸にChromaの散布図を作成します。

bs = 32 # block size
s = 1
e = s + bs - 1

plot(lena[s:e,s:e], fmt=:png)



x_min = minimum(lena_arr[1,s:e,s:e])
x_max = maximum(lena_arr[1,s:e,s:e])

# Cb
x = vec(lena_arr[1,s:e,s:e])
y = vec(lena_arr[2,s:e,s:e])
scatter(x, y, label="Cb", xlabel="Luma", ylabel="Chroma")

cb_a = (bs*bs*sum(x.*y) - sum(x)*sum(y)) / (sum(bs*bs*(x.^2)) - sum(x)^2)
cb_b = (sum(y) - cb_a * sum(x)) / (bs*bs)
x = range(x_min, x_max, length=1000)
y = cb_a .* x .+ cb_b
plot!(x, y, label="Cb prediction")

# Cr
x = vec(lena_arr[1,s:e,s:e])
y = vec(lena_arr[3,s:e,s:e])
scatter!(x, y, label="Cr")

cr_a = (bs*bs*sum(x.*y) - sum(x)*sum(y)) / (sum(bs*bs*(x.^2)) - sum(x)^2)
cr_b = (sum(y) - cr_a * sum(x)) / (bs*bs)
x = range(x_min, x_max, length=1000)
y = cr_a .* x .+ cr_b
plot!(x, y, label="Cr prediction", fmt=:png)

f:id:suzuzusu:20200304055625p:plain

f:id:suzuzusu:20200304055638p:plain

上記の結果から、ChromaのそれぞれのCr, CbはLumaに対してのある程度線形な関係にあることが分かりました。つまりCr, Cbはそれぞれ傾きと切片が分かっているならLumaからある程度推定可能ということです。この関係を用いてCfLを実装します。

Chroma from Luma Prediction (CfL)

先程の例でChromaとLumaの関係が分かっているのであとは定式化してみましょう。 Predicting Chroma from Luma in AV1 から式を引用するとCfLは以下のようなパラメータ $\alpha, \beta$ の線形モデルを使って推定します。 $C^{p}, L^{r}$ はChroma, Lumaです。

f:id:suzuzusu:20200304061247p:plain

$\alpha, \beta$ パラメータは次のように最小二乗法で求めます。

f:id:suzuzusu:20200304061319p:plain

最後に、8, 16, 32のブロックサイズごとにCfLをした画像とオリジナル画像を比較してみましょう。

function cfl(img_arr, bs=8)
    nc, h, w = size(img_arr)
    img_arr_est = zeros((nc, h, w))
    
    for i in 1:Int(h/bs)
        si = (i-1)*bs + 1
        ei = i*bs
        for j in 1:Int(w/bs)
            sj = (j-1)*bs + 1
            ej = j*bs
            img_arr_tmp = img_arr[:, si:ei, sj:ej]

            # Cb
            x = vec(img_arr_tmp[1, :, :])
            y = vec(img_arr_tmp[2, :, :])
            cb_a = (bs*bs*sum(x.*y) - sum(x)*sum(y)) / (sum(bs*bs*(x.^2)) - sum(x)^2)
            cb_b = (sum(y) - cb_a * sum(x)) / (bs*bs)

            # Cr
            x = vec(img_arr_tmp[1, :, :])
            y = vec(img_arr_tmp[3, :, :])
            cr_a = (bs*bs*sum(x.*y) - sum(x)*sum(y)) / (sum(bs*bs*(x.^2)) - sum(x)^2)
            cr_b = (sum(y) - cr_a * sum(x)) / (bs*bs)

            x = img_arr_tmp[1, :, :]
            cb_est = cb_a .* x .+ cb_b
            cr_est = cr_a .* x .+ cr_b
            img_arr_est_tmp = zeros((3, bs, bs))
            img_arr_est_tmp[1,:,:] .= x
            img_arr_est_tmp[2,:,:] .= cb_est
            img_arr_est_tmp[3,:,:] .= cr_est

            img_arr_est[:, si:ei, sj:ej] .= img_arr_est_tmp
        end
    end
    img_arr_est
end

p_org = plot(colorview(YCbCr, lena_arr), title="origin")
p8 = plot(colorview(YCbCr, cfl(lena_arr, 8)), title="CfL(block size 8)")
p16 = plot(colorview(YCbCr, cfl(lena_arr, 16)), title="CfL(block size 16)")
p32 = plot(colorview(YCbCr, cfl(lena_arr, 32)), title="CfL(block size 32)")
plot(p_org, p8, p16, p32, layout = (1, 4), size = (1200, 300), fmt=:png)

f:id:suzuzusu:20200304055719p:plain

上記の結果からオリジナル画像に対して、CfLでかなりきれいに推定できることが分かりました。 AV1では、上記のような基本的なCfLが用いられているわけではないのですが原理は同じです。さらに今後の取り組みとして非線形モデルを使って推定するなども考えられているらしいです。

f:id:suzuzusu:20200304062355p:plain

notebook

nbviewer.jupyter.org

参考

arxiv.org
people.xiph.org
Chroma from Luma Intra Prediction for AV1 from Luc Trudeau
www.slideshare.net

2019-12-15

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

機械学習系技術系

fableとは？

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 — 結果（拡大）

参考

2019-12-13

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

参考

2019-12-04

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

機械学習系技術系

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

線形分離不可能なデータ

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

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

参考

2019-11-24

幾何ブラウン運動を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"))

参考

2019-11-17

尤度関数におけるガウス分布とスチューデントの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);
}