ローカルレベルモデルとカルマンフィルタ

本to実装

ちょっとだけ話題（？）になっていたカルマンフィルタを参考文献の本で学んでいます。話題になっているから学んでいるわけではないですが、時系列分析の基本だと思うので。

今回は第2章のローカルレベルモデルに対するカルマンフィルタを実装します。ローカルレベルモデルとは次のような最もシンプルな状態空間モデルのことです。
観測値を $y_t$ 、観測できない潜在変数を $\alpha_t$ として
f:id:mt19:20190519101423j:plain
という仮定を置いたモデルです。ただし、これだと初期値について何も言及がないので
f:id:mt19:20190519101430j:plain
のように初期分布を置いてあげます。従って未知パラメータは $\sigma^2_{\epsilon}, \sigma^2_{\eta}, a_1, P_1$ の4つです。未知パラメータのうち初期分布に関わる $a_1, P_1$ については、何らかの事前知識があれば良いのですが今回は無いのでいわゆる無情報事前分布にしてしまいます。つまり、分散 $P_1$ を大きくしてしまいます。このとき $a_1$ はもはやどこでも変わらないので、常識的な範囲で適当に置いてもらって構いません。他のパラメータについては最尤推定で決める方法があるようですが、今回は既知として自分で設定しました。

使うデータは参考文献の本から提供されている体重の推移データです。データを使う場合は出版社のページから飛んでください。
カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― / 野村俊一著 | 共立出版

ローカルレベルモデルはシンプルであり、stanを使えばあっという間に実装できそうな気がしますが、今回はあくまでカルマンフィルタのお勉強です。ということで、アルゴリズムを導出して、その式を使って実装しています。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

#load data and plot
weight_data = pd.read_csv("./Weight.dat", header = None)

plt.figure(figsize = (8, 4))
plt.plot(weight_data)
plt.xlabel("day")
plt.ylabel("weight(kg)")
plt.grid(True)
plt.show()

f:id:mt19:20190519103005p:plain

ローカルレベルモデルのカルマンフィルタを実装です。

class LocalLevelModel():
    """
    this class gives LocalLevelModel using KalmanFilter.
    model accepts only univariate pandas data.
    if missing data is described np.NaN, model deals with it automatically
    """
    def __init__(self, data, a_ini = 0.0, P_ini = 1000, ep_ini = 1.0, eta_ini = 0.1):
        self.data = data.values
        self.n = len(data)
        self.a_ini = a_ini
        self.P_ini = P_ini
        self.ep_ini = ep_ini
        self.eta_ini = eta_ini
        self.missing_ind = np.where(np.isnan(self.data) == True)[0]
        self.dropna_data = self.data[~np.isnan(self.data)]
    
    
    #calculate filtered estimator
    def fit(self):
        v, F, a_con, a = np.empty(self.n), np.empty(self.n), np.empty(self.n), np.empty(self.n)
        P, P_con, K, L = np.empty(self.n), np.empty(self.n), np.empty(self.n), np.empty(self.n)

        #初期値代入
        a[0] = self.a_ini
        P[0] = self.P_ini
        
        
        #逐次計算
        for t in range(self.n - 1):
            F[t] = P[t] + self.ep_ini ** 2
            K[t] = P[t] / F[t]
            L[t] = 1 - K[t]
            
            if t in self.missing_ind:
                a_con[t] = a[t]
                P_con[t] = P[t]
                a[t + 1] = a_con[t]
                P[t + 1] = P_con[t] + self.eta_ini ** 2
            else:
                v[t] = self.data[t] - a[t]
                a_con[t] = a[t] + K[t] * v[t]
                P_con[t] = P[t] * L[t]
                a[t + 1] = a_con[t]
                P[t + 1] = P_con[t] + self.eta_ini ** 2
        
        t = t + 1
        F[t] = P[t] + self.ep_ini ** 2
        K[t] = P[t] / F[t]
        L[t] = 1 - K[t]
        if t in self.missing_ind:
            a_con[t] = a[t]
            P_con[t] = P[t]
        else:
            v[t] = self.data[t] - a[t]
            a_con[t] = a[t] + K[t] * v[t]
            P_con[t] = P[t] * L[t]
        
        return a_con, P_con, F, v, L
    
    #state smoothing
    def smoothing_fit(self):
        #smoothed state
        hat_alpha = np.empty(self.n)
        #smoothed state variance
        V = np.empty(self.n)
        
        r = np.empty(self.n)
        N = np.empty(self.n)
        
        #initialize
        r[self.n - 1] = 0
        N[self.n - 1] = 0
        
        a_con, P_con, F, v, L = self.fit()
        
        #逐次計算
        for t in range(self.n - 1):
            if t in self.missing_ind:
                r[self.n - 2 - t] = r[self.n - 1 - t]
                hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
                
                N[self.n - 2 - t] = N[self.n - 1 - t]
                V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
            else:
                r[self.n - 2 - t] = v[self.n - 1 - t] / F[self.n - 1 - t] + L[self.n - 1 - t] * r[self.n - 1 - t]
                hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
            
                N[self.n - 2 - t] = 1 / F[self.n - 1 - t] + (L[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
                V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
            
        t = t + 1
        hat_alpha[self.n - 1 - t] = a_con[self.n - 1 - t] + P_con[self.n - 1 - t] * r[self.n - 1 - t]
        V[self.n - 1 - t] = P_con[self.n - 1 - t] - (P_con[self.n - 1 - t] ** 2) * N[self.n - 1 - t]
        
        return hat_alpha, V

変数名は参考文献に寄っているので、本が無いと分かりにくいかもしれません。重要な変数については説明を置いているつもりです。欠測値（np.NaN）にも自動で対応するようにしました。

fit()ではフィルタ化推定量とその分散、smoothing_fit()では平滑化状態とその分散を返してくれます。

まずは体重の推移データで、欠測値が一切無いときをプロットしてみます。
f:id:mt19:20190526120150p:plain
f:id:mt19:20190526120158p:plain

平滑化状態は全ての観測データを利用した予測なので、フィルタ化推定量に比べて幅が狭くなっているのが分かります。

次に欠測値を混ぜてみます。20日から40日の間は体重を測定しなかったと仮定して、データをモデルに投げます。
f:id:mt19:20190526120227p:plain
f:id:mt19:20190526120234p:plain
欠測値のあたりだけ不自然に膨らんでいます。データが無いので幅が広くなっていることが分かります。

上の欠測値のあるときのフィルタ化推定量というのは実質n期先の予測に他なりません。ローカルレベルモデルでは直前の状態を参考にするしかないため、そうすると「直前の値をひたすらまっすぐ取り続けるであろう」というなんとも情けない予測になってしまいます。ここが良くないですね。

おまけ１（2019/05/25）

ためしに外れ値を混ぜてみたらどういう結果になるか調べてみる。
f:id:mt19:20190526120315p:plain
f:id:mt19:20190526120323p:plain
上がフィルタ化推定量、下が平滑化状態の推定量です。思いのほか引っ張られていない？よく分からぬ。

おまけ２（2019/05/25）

TFPでフィルタ化推定量を求めてみた。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import tensorflow_probability as tfp
import tensorflow as tf
import tensorflow.contrib.eager as tfe
tfe.enable_eager_execution()
    
tfl = tf.linalg
tfd = tfp.distributions


weight_data = pd.read_csv("./Weight.dat", header = None)
weight_data = weight_data.values

num_dims = 1
num_timesteps = weight_data.shape[0]
step_mu = 0.0
step_std = 1.0
noise_mu = 0.0
noise_std = 1.0
ini_loc = tf.convert_to_tensor(np.array([0], dtype = np.float32))
ini_scale = tf.convert_to_tensor(np.array([1000]), dtype = np.float32)
scale = tf.convert_to_tensor(np.array([0.1]), dtype = np.float32)
model = tfd.LinearGaussianStateSpaceModel(
    num_timesteps = num_timesteps,
    transition_matrix = tfl.LinearOperatorIdentity(num_dims) ,
    transition_noise = tfd.MultivariateNormalDiag(loc = tf.zeros([num_dims]), scale_diag = scale),
    observation_matrix = tfl.LinearOperatorIdentity(num_dims),
    observation_noise = tfd.MultivariateNormalDiag(loc = tf.zeros([num_dims]), scale_diag = tf.ones([num_dims])),
    initial_state_prior = tfd.MultivariateNormalDiag(loc = ini_loc, scale_diag = ini_scale)
)

input_data = tf.convert_to_tensor(weight_data.astype(np.float32))
output_data = model.forward_filter(input_data)
prediction = np.array(output_data[1]).flatten()
std = np.sqrt(np.array(output_data[2]).flatten())
x = range(len(weight_data))
z = stats.norm.ppf(q = 0.975)
plt.figure(figsize = (8, 4))
plt.scatter(x, weight_data)
plt.plot(weight_data)
plt.plot(prediction)
plt.fill_between(x, prediction - z * std, prediction + z * std, alpha = 0.2)
plt.grid(True)
plt.show()

f:id:mt19:20190526121640p:plain
たぶん同じように動いてくれているはず・・・。
もしコピーしてjupyterなどで自分で動かしてみようという物好きな方がいれば注意してください。kernelを再起動してまっさらな状態でeager_executionを実行しないとエラーが出ます。
python - tf.enable_eager_execution must be called at program startup ONLY in SPYDER IDE - Stack Overflow

参考文献

カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)

作者: 野村俊一
出版社/メーカー: 共立出版
発売日: 2016/09/08
メディア: 単行本
この商品を含むブログを見る

この記事については第2章を参考にしました。

2019-04-14

ジェフリーズ事前分布(Jeffreys prior distribution)

数理統計

たまにはブログ更新しないとと思ったので。
間違いがあればtwitterにでも気軽にリプくださると嬉しいです。

定義

ジェフリーズ事前分布という名前のついた事前分布があります。
en.wikipedia.org
ベイズ推論において尤度のパラメータ（多変量でも可）を $\theta$ としたとき次のような事前分布のことを言います。
\begin{equation}
p(\theta) \propto |J(\theta)|^{1/2}
\end{equation}
ただし $J(\theta)$ はフィッシャー情報量です。フィッシャー情報量って何ぞという方はwikipediaでチェックしましょう。
ja.wikipedia.org

これは何？

この定義の何が嬉しいかというと変数変換に対して分布が不変なのです。つまり、ジェフリーズ事前分布を変数変換するとそれもまたジェフリーズ事前分布になります（ただし適当な条件下で）。実際に書き下してみましょう。変数変換を $\phi = h(\theta)$ としたとき、
\begin{align}
p(\phi) &= p(\theta) |d\theta / d\phi| = p(\theta) |h'(\theta)|^{-1}\\
J(\phi) &= -E\left [\frac{d^2\log{p(y | \phi)})}{d\phi^2}\right ] \\
&= -E\left [ \left(\frac{d^2\log{p(y | \theta = h^{-1}(\phi))}}{d\theta^2} \right) \left (\frac{d\theta}{d\phi} \right)^2 \right] = J(\theta) \left (\frac{d\theta}{d\phi} \right)^2
\end{align}
となります。従って $p(\theta)$ がジェフリーズ事前分布ならば、
\begin{equation}
p(\phi) \propto |J(\theta)|^{1/2} |d\theta / d\phi| = |J(\phi)|^{1/2}
\end{equation}
となり、確かに変数変換してもジェフリーズ事前分布となります。

具体例

分散が既知の正規分布 $N(\mu, \sigma_0^2)$ を考えます。この事前分布としてジェフリーズ事前分布 $p(\mu)$ を考えます。定義通りフィッシャー情報量を求めれば十分です。
f:id:mt19:20190414002607p:plain
※何故かブログのtexが反応しないので画像で。（追記：期待値の外側に√が抜けています。結果には問題ありません。）
このように定数になります。定数というのはよくある無情報事前分布です。しかしパラメータは $\mu \in (-\infty, \infty)$ を動くのでジェフリーズ事前分布を積分すると発散してしまい積分が1という確率の公理が満たされません。このような事前分布はimproperな事前分布と言われます。実際のところ、そんなに広いところまで $\mu$ が動くとは思えないので、ある程度の範囲で制限して閉区間にしたら良いと思います。そうすれば確率の公理を満たすproperな事前分布になります。

この具体例は変数変換に対して不変ということに何も触れていないので、具体例としてはイマイチです。肝心なこの定義の嬉しさを何も伝えられていません。正直なところ、どういうシチュエーションでこの事前分布がうまく働くのか私はあまり理解できていないので言及はしないでおきます。ちなみに参考文献2では少し言及されていて、どうやらこの事前分布を使える状況はかなり狭そうです。

何か分かればまた追記しますね。

参考文献

1.Andrew Gelman『Bayesian Data Analysis 3rd Edition』Chapman and Hall/CRC
2.8Noninformative prior distributionsに説明があります。

2. 渡辺澄夫『ベイズ統計の理論と方法』コロナ社
定義はp169にあります。

2018-12-17

統計検定1級に合格しました

雑記

合格したので感想を。合格した勉強法の1例として参考になれば。
2018年11月下旬に統計検定1級を受験して統計数理と統計応用（理工学を選択）の両方に合格しました。以下、その受験記です。
*1

背景

自分のことを少しだけ書きます。
普段は数学を専攻している大学生です。とは言っても大学で統計学を学んだことはなく、2018年3月ごろから趣味で勉強を始めました。2018年5月には統計検定2級に合格しています。良ければその受験記も見て下さい。
統計検定2級に合格しました - lazy diary
さて、「2級に合格したなら次は1級でしょ」ということで受験しました。準1級を飛ばしたのは単純に時期が合わなかったからです。

統計検定1級の概要

年1回11月下旬開催の試験です。試験範囲をチェックしましょう。
www.toukei-kentei.jp
試験内容は統計数理と統計応用の2つに分かれています。
この両方に合格すると統計検定1級合格となります。

・統計数理

試験日の午前に行なわれる90分の筆記試験です。60/100点で合格です。
5つの問題から3問を選んで解きます。確率分布に関わる様々な計算がメインです。
確率密度関数の計算、確率母関数を用いた計算、検定論が頻出です。

・統計応用

試験日の午後に行なわれる90分の筆記試験です。60/100点で合格です。
統計応用は「人文科学」「社会科学」「理工学」「医薬生物学」の中から1つの分野を試験申し込み時に選んで解くことになります。
事前に選んだ分野の5つの問題から3問を選んで解きます。

勉強したこと

勉強期間は3ヵ月くらいです。2ヵ月ほど本で勉強、残りは演習といった感じですかね。

さて、まず参考書選びです。1級に合格したというブログを調べるとちらほら出てくるのでそこにオススメしてある参考書と、あとは公式サイトにある合格者の声にリストアップされている参考書を書店でチェックしました。
www.toukei-kentei.jp
その中で私が選んだのがこちらの本。

・現代数理統計学の基礎 (共立講座数学の魅力)

現代数理統計学の基礎 (共立講座数学の魅力)

作者: 久保川達也
出版社/メーカー: 共立出版
発売日: 2017/04/07
メディア: 単行本
この商品を含むブログ (1件) を見る

この本が無ければ合格は無かったです。それくらい良い本でした。（※ゴリゴリの数式に慣れていない方は少々難しいかもしれません。）
統計検定1級統計数理の範囲となっているが本には掲載されていない事項の一部を著者がWebページに補足資料として載せています。それくらい統計検定に最適化された本です。この本で統計数理のほとんどの試験範囲を抑えることができます。以下、『現代数理統計学の基礎』の目次です。

第1章　確率
第2章　確率分布と期待値
第3章　代表的な確率分布
第4章　多次元確率変数の分布
第5章　標本分布とその近似
第6章　統計的推定
第7章　統計的仮説検定
第8章　統計的区間推定
第9章　線形回帰モデル
第10章　リスク最適性の理論
第11章　計算統計学の方法
第12章　発展的トピック：確率過程
引用：https://www.kyoritsu-pub.co.jp/bookdetail/9784320111660

第1章から第9章までは十分に理解しておくべきでしょう。試験問題の半分くらいは適切に確率密度関数や期待値、分散などの計算ができるかどうかです。この計算がスムーズにできるかどうかが鍵だと思います。
第10章「リスク最適性の理論」は飛ばしても問題ないと思います。
第11章、第12章についても理解しておくべきでしょう。ただし、マルチンゲールといった試験に出ないであろう内容もあるので試験範囲と照らし合わせて必要に応じていくつか飛ばしても良いと思います。

この本を読むのに1ヵ月と少しかかりました。さっそく過去問に取り組みます。2017年の統計数理の問題を解きました。自己採点をしたら合格ライン上の点数だったと思います。それより時間が足りません。明らかに演習不足なのを実感しましたが、同時に解説を読めば理解できるということも分かったので統計数理の勉強はひとまず終了にしました。

次に統計応用の勉強をします。私は「理工学」分野を選択しました。理由は統計数理とさほど範囲が変わらないからです。過去問を見てもらえば分かると思います。「理工学」以外の分野に明るい場合はそれを選べば良いと思います。特にそういうのがないですという場合は、「理工学」分野を選ぶのが無難だと私は思います。
参考書は次の2冊です。

・多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

作者: 小西貞則
出版社/メーカー: 岩波書店
発売日: 2010/01/27
メディア: 単行本（ソフトカバー）
購入: 14人クリック: 347回
この商品を含むブログ (9件) を見る

１　はじめに
２　線形回帰モデル
３　非線形回帰モデル
４　ロジスティック回帰モデル
５　モデル評価基準
６　判別分析
７　ベイズ判別
８　サポートベクターマシーン
９　主成分分析
10　クラスター分析
引用：https://www.iwanami.co.jp/book/b265441.html

この本はかなり丁寧に記述されています。統計検定云々は関係なく、多変量解析を学びたい人にはオススメしたい1冊でした。
（※もし「試験日まで時間がない」とか「とにかく近道して受かりたい」といった場合はこの本は飛ばしても問題ないでしょう。これら多変量解析の内容は準1級および1級試験の試験範囲に含まれています。従って試験に出題されてもおかしくないです。実際に過去問を見てみると準1級試験ではこれらの出題事例があります。例えば、2018年準1級試験ではマーク問題(not記述式)として線形判別やSVMが出題されています。一方で1級試験では私の知る限り出題事例が無く、せいぜい線形重回帰くらいだったと思います。）

・日本統計学会公式認定統計検定1級対応統計学

日本統計学会公式認定統計検定1級対応統計学

作者: 二宮嘉行,大西俊郎,小林景,椎名洋,笛田薫,田中研太郎,岡田謙介,大屋幸輔,廣瀬英雄,折笠秀樹,日本統計学会,竹村彰通,岩崎学
出版社/メーカー: 東京図書
発売日: 2013/04/08
メディア: 単行本
この商品を含むブログ (5件) を見る

第I部　統計数理
第1章　確率と確率変数
第2章　種々の確率分布
第3章　統計的推定
第4章　仮説検定
第5章　主なデータ解析手法
第II部統計応用
第6章　統計応用共通手法
第7章　人文科学分野キーワード
第8章　社会科学分野キーワード
第9章　理工学分野キーワード
第10章　医学生物学分野キーワード
引用：http://www.tokyo-tosho.co.jp/books/ISBN978-4-489-02150-3.html

次は公式の参考書です。上2冊で学習できていない分散分析や直交表などを補いつつ、これまで学んだ事を一気に復習します。この本は某レビューで酷評されているようですが、私は使い方の問題だと思っていて、この本でイチから学ぼうとすると確かに難しい気がしますが、復習に使うならコンパクトに試験範囲がまとまっているので良い本だと思います。

さて、これで統計応用の範囲も全てチェックしました。
私は時間にゆとりがあったのでこの間にもいくつか別の本をやっていますが、これだけで十分だと思います。後はひたすら演習です。時間をきっちり計って過去問を解きまくりましょう。
試験日まで時間がある人はこのあたりで気づくと思うのですが、年1回しか試験がないので過去問が少ないです。もっと演習をしたいという人にはこれが良いかもしれません。
Exam past papers
王立統計学会（RSS）が実施している試験の過去問です。これも統計検定公式サイトに載っている合格者の声で知りました。このうち「Graduate Diploma」の過去問を解いてみるのが良いと思います。難易度は1級よりやや簡単ですが、問題に慣れるのには十分でしょう。実際に似たような問題が1級試験でも出題されています。

受験当日

受験票に写真を貼ることと、電卓、時計を忘れないようにしましょう。
会場に行きます。最初にアンケートがあるので回答しましょう。
さっそく統計数理の試験開始です。どこかで見たことがある問題を見つけ1問サクッと解答。残り2問は最後がダメでしたがそれ以外は一通り解けたと思います。そしてお昼休憩を挟んで、午後から統計応用の試験です。
計算量が多そうな条件付き正規分布の計算問題はスルーして、とにかく分かる問題を落とさずしっかり得点するという気持ちで解いていきました。

12月半ばにWeb合格発表が行われ無事合格しました。

最後に

統計検定1級を通して、いわゆる古典的な統計は一通り学べたと思っています。今はDeep Learningが流行っていますが、ベイジアンになりたいなら確率関数の操作だったりはそれらを理解するための基本なので、その理解度をチェックする上で統計検定を受験してみてはいかがでしょうか。

*1:統計検定の公式キャラクターらしいです http://www.toukei-kentei.jp/info/toukeiakari/

2018-10-22

非線形SVMによる分類

本to実装

参考書はこちら

多変量解析入門――線形から非線形へ

作者: 小西貞則
出版社/メーカー: 岩波書店
発売日: 2010/01/27
メディア: 単行本（ソフトカバー）
購入: 14人クリック: 347回
この商品を含むブログ (9件) を見る

ハードマージンSVM、ソフトマージンSVMと続いて今回でSVMは最後です。
今までは全て線形SVMでしたが、今回は非線形のSVMとなります。非線形 SVMでは特徴空間次元の内積計算が多数必要になり、それを避けるためにカーネル法という手法が用いられます。カーネル法にも種類は様々あるのですが今回は最も一般的なガウスカーネルを用います。データはソフトマージンSVMのときと同じデータです。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class NonLinearSoftMarginSVM():
  
  def __init__(self, X, y, lam, sigma):
    self.X = X
    self.y = y
    self.lam = lam
    self.sigma = sigma
    
  #ガウスカーネル
  def kernel(self, x, y):
    return np.exp(np.dot(x - y, x - y) / (-2 * self.sigma))
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * self.kernel(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      for j in range(num_of_data):
        b = b + alpha[j] * label.iloc[j] * self.kernel(self.X.iloc[j], self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index, alpha]

前回のソフトマージンSVMと比べて変化しているのは、内積計算が全てカーネルによる計算に置き換えられているという点だけです。続いて、ここにデータを入れてみます。

lam = 1.0
sigma = 1.0
X = data[["x", "y"]]
y = data["c"]
label = data["c"] * 2 - 1
svm = NonLinearSoftMarginSVM(X , y, lam, sigma)
[w, b, Sup_index, alpha] = svm.fit()

ここまでは実はさほど計算時間はかかりません。決定境界をプロットしてみましょう。

def f(x, alpha, Sup_index, b, X, label):
  result = b
  for i in range(len(Sup_index)):
    j = Sup_index[i]
    result = result + alpha[j] * label[j]  * svm.kernel(X.iloc[j], x)
  return result

#色の設定
col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")

#メッシュグリッドを作成して決定境界をプロット
seq = np.linspace(0, 10, 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
for i in range(X1.size):
  Z.append(f(np.array([x1[i][0], x2[i][0]]), alpha, Sup_index, b, X, label))

plt.title("Non Lnear Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.0])
plt.scatter(data["x"], data["y"], c = col)

f:id:mt19:20181022155837p:plain
このように非線形の決定境界となります。SVM自体の話はここまでで、あとはどうやってモデルを良いものにしていくかが難しいところです。例えば今回はソフトマージンなので正則化係数、さらにはガウスカーネルなので分散 $\sigma$ を事前に設定する必要があります。他には多項式カーネルであれば2つのハイパーパラメータを設定しなければなりません。クロスバリデーションなどを適用しながらこれらを適切に設定することが求められます。

参考

非線形の決定境界を書く際に参考になりました。
python で等高線を描くなら meshgrid して contour | コード７区

2018-10-22

ソフトマージンSVMによる分類

本to実装

参考書はこちら

多変量解析入門――線形から非線形へ

作者: 小西貞則
出版社/メーカー: 岩波書店
発売日: 2010/01/27
メディア: 単行本（ソフトカバー）
購入: 14人クリック: 347回
この商品を含むブログ (9件) を見る

前回はハードマージンSVMを実装しました。今回はそれの拡張のようなもので、ソフトマージンSVMを行ないます。ハードマージンSVMではデータは完全に線形分離できるものにしてくださいと書きました。
ハードマージンSVMによる分類 - lazy diary
しかし現実問題そんな都合の良いデータばかりではないでしょう。なんらかの理由でラベルの異なるグループに入ってしまうデータもあるはずです。ソフトマージンSVMでは、そのようなデータが入ることを許すよう条件を緩和します。しかし、どの程度まで許すかということが問題になります。そこで条件を緩和しすぎず、なおかつマージンを大きくすることを目指したのがソフトマージンSVMです。
実際に解析するにあたってその程度を決める指標として正則化係数を私たちは与える必要があります。これについてはクロスバリデーションで決める方法がおそらく最もシンプルだと思われます。

コード自体はハードマージンSVMとほとんど変わりません。二次計画問題で、条件が増える点と、サポートベクトル集合にも条件がさらに加わることだけ注意します。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class SoftMarginSVM():
  
  def __init__(self, X, y, lam):
    self.X = X
    self.y = y
    self.lam = lam
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

データを作ります。ハードマージンSVMと違い多少は入り乱れたデータでも平気です。

#create data
x_1 = 4.2
y_1 = 7
x_2 = 7
y_2 = 4.2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

これをSVMにかけてみましょう。

lam = 10.0
svm = SoftMarginSVM(data[["x", "y"]], data["c"], lam)
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.title("Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022142912p:plain
このように、ハードマージンSVMの厳しすぎる条件を緩和することでソフトマージンSVMでは分類を達成します。

2018-10-22

ハードマージンSVMによる分類

本to実装

参考書はこちら

多変量解析入門――線形から非線形へ

作者: 小西貞則
出版社/メーカー: 岩波書店
発売日: 2010/01/27
メディア: 単行本（ソフトカバー）
購入: 14人クリック: 347回
この商品を含むブログ (9件) を見る

今回はハードマージンSVMというものを実装します。
ソフトでもない上に非線形でもないという最も基本的なSVMです。ハードマージンSVMさえ理解してしまえば他のSVMはその拡張みたいなものです。

必要なものをimportしておきます。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

このcvxoptというのは二次計画問題のソルバーです。詳しくはこの記事の最後にあるブログを参考にしてください。google colabにもデフォルトではインストールされていないようなので、

pip install cvxopt

で導入しておきましょう。

データは相変わらずいい加減に作ったものです。とは言いつつも今回はハードマージンSVMですので、線形分離できないデータは困ります。かならず超平面で分離できるデータをチョイスしてください。2次元なら直線で完全に分離できるものです。

#create data
x_1 = 3
y_1 = 7
x_2 = 7
y_2 = 2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

ここからが本編です。流れとしては
1. データから二次計画問題を解く
2. パラメータwを求める
3. サポートベクトル集合を求める
4. サポートベクトル集合から切片bを求める
となっています。

class HardMarginSVM():
  
  def __init__(self, X, y):
    self.X = X
    self.y = y
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    G = matrix((-1) * np.identity(num_of_data))
    h = matrix(np.zeros((num_of_data, 1)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon:
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

これを上のデータに当てはめてみましょう。

svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022105408p:plain
黄色になっている部分はサポートベクトルです。SVMの特徴として実質的にサポートベクトルにしか依存しないので、次のような外れ値が入っても影響がありません。

data.loc[100] = [0, 20, 1]
svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022110158p:plain

参考

今回、二次計画問題を解くにあたって参考にしました。SMOアルゴリズムについてはまだ理解できていないため、今回はcvxoptというソルバーに頼りました。
二次計画法(Quadratic Programming)の概要とPythonソルバーcvxoptの使い方 - MyEnigma

2018-10-11

XORのようなデータを非線形ロジスティック回帰で分類する

本to実装

参考書は相変わらずこちら

多変量解析入門――線形から非線形へ

作者: 小西貞則
出版社/メーカー: 岩波書店
発売日: 2010/01/27
メディア: 単行本（ソフトカバー）
購入: 14人クリック: 347回
この商品を含むブログ (9件) を見る

今回は、基底関数展開を利用した非線形のロジスティック回帰を行ないます。基底関数を求めて、データを変換、最後にロジスティック回帰を行なうだけです。今回は基底関数としてガウスカーネルを使いました。
非線形の恩恵を確かめるために、通常の線形ロジスティック回帰では分類できないXORのようなデータを今回は使います。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.cluster import KMeans
from sklearn import linear_model

from scipy.stats import multivariate_normal

cov = 0.5 * np.identity(2)
data1 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, -2.0), cov = cov, size = 25))
data2 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, -2.0), cov = cov, size = 25))
data3 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, 2.0), cov = cov, size = 25))
data4 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, 2.0), cov = cov, size = 25))

data1["c"] = 0
data2["c"] = 1
data3["c"] = 0
data4["c"] = 1

data = pd.concat([data1, data2, data3, data4], join = "inner")
data = data.reset_index(drop = True)
data = data.rename(columns = {0:"x", 1:"y", "c":"c"})

x = data["x"]
y = data["y"]
clusters = data["c"]

plt.scatter(x, y, c = clusters)

f:id:mt19:20181109161138p:plain
色付けをまともにしてないので見づらいですが、たしかにXORです。ここから前回と同じようにk-meansからの基底関数の構成およびデータの変換を行ないます。

#k-meansなどでクラスタリングする
dist = []
for i in range(1,30):
  km = KMeans(n_clusters = i)
  km.fit(data[["x", "y"]])
  dist.append(km.inertia_)
  
plt.plot(range(1,30), dist, marker = ".")

f:id:mt19:20181109161217p:plain
エルボー則に従って4くらいかなと考えられます。

from sklearn import linear_model
#クラスタ数を4とする

#ガウス型基底関数
def gauss_func(x, mu, sigma):
  return np.exp((np.sum(np.power(x - mu, 2))) / (-2 * sigma))

m = 4
km = KMeans(n_clusters = m)
km.fit(data[["x", "y"]])
cluster_center = km.cluster_centers_
dist = np.power(km.transform(data[["x", "y"]]), 2)
cluster_num = km.predict(data[["x", "y"]])
least_dist = np.zeros(shape = (m, 1))
for i in range(len(data)):
  for j in range(m):
    if cluster_num[i] == j:
      least_dist[j] += dist[i][j]
for i in range(m):
  least_dist[i] = least_dist[i] / np.sum(cluster_num == i)

transformed_X = pd.DataFrame(np.zeros(shape = (len(data), m)))
for i in range(len(data)):
  for j in range(m):
    transformed_X.iloc[i, j] = gauss_func(data[["x", "y"]].iloc[i,:], cluster_center[j], least_dist[j])

さて、カーネルパラメータの推定とカーネルを使ってデータの変換まで終わったのでロジスティック回帰を実行します。

#ロジスティック回帰の実行
log_reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga", C = 1000.0, max_iter = 1000).fit(transformed_X, clusters)
print(log_reg.score(transformed_X, clusters)) #->1.0

とりあえず分類できました。決定境界をプロットしてみます。

seq = np.linspace(data["x"].min(), data["x"].max(), 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
testdata = np.hstack((x1, x2))
transformed_testdata = pd.DataFrame(np.zeros(shape = (len(x1), m)))
for i in range(len(x1)):
  for j in range(m):
    transformed_testdata.iloc[i, j] = gauss_func(testdata[i,:], cluster_center[j], least_dist[j])
Z = log_reg.predict_proba(transformed_testdata)

plt.title("Non Linear Logistic Regression")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z[:,0]).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.5])
plt.scatter(data["x"], data["y"], c = clusters)

f:id:mt19:20181109161626p:plain
ガウスカーネルっぽい丸みを帯びた決定境界です。
基底関数になるものならばどれを使ったって構いません。（うまく分類できるかは別ですが）。例えば $x^n(n = 0, 1, 2,\dots)$ を用いても良いです。

おまけ１（2019/05/25）

おまけ２（2019/05/25）

参考文献

定義

これは何？

具体例

参考文献

背景

統計検定1級の概要

・統計数理

・統計応用

勉強したこと

・現代数理統計学の基礎 (共立講座 数学の魅力)

・多変量解析入門――線形から非線形へ

・日本統計学会公式認定 統計検定1級対応 統計学

受験当日

最後に

参考

参考

・現代数理統計学の基礎 (共立講座数学の魅力)

・日本統計学会公式認定統計検定1級対応統計学