2017-11-05

Pythonによる「べき乗法（Power Method）」の実装 - 最大特異値・特異ベクトルを求める

Python

はじめに

　理工学の分野において、行列の固有値と固有ベクトルを求める問題に直面することが多々あります。一般に、固有値や固有ベクトルは固有値分解という行列分解で求められますが、これは正方行列に対してのみ定義されています。特異値分解は、固有値分解を長方形行列に拡張したものであり、より適用範囲が広く、画像圧縮や自然言語処理における潜在意味解析などにも応用されています。特異値分解は比較的計算コストが高いため、例えば絶対値最大の特異値やその特異ベクトルを求めたいだけの場合は、他の手法を使った方が計算効率が高いです。
　今回はべき乗法（Power method）と呼ばれる数値計算アルゴリズムをPythonで実装し、行列の最大特異値および特異ベクトルを求める方法について書いていきます。

Pythonによる実装

　べき乗法の理論については、下記リンクの資料（P5〜）を参考にしました。

http://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf

import numpy as np

def power_method(A, iter_num=1):
    """
    Calculate the first singular vector/value of a target matrix based on the power method.
    Parameters
    ----------
    A : numpy array
        Target matrix
    iter_num : int
               Number of iterations

    Returns
    -------
    u : numpy array
        first left singular vector of A
    s : float
        first singular value of A
    v : numpy array
        first right singular vector of A
    """
    # set initial vector q
    q = np.random.normal(size=A.shape[1])
    q = q / np.linalg.norm(q)

    for i in range(iter_num):
        q = np.dot(np.dot(A.T, A), q)

    v = q / np.linalg.norm(q)
    Av = np.dot(A, v)
    s = np.linalg.norm(Av)
    u = Av / s

    return u, s, v

特異値分解との比較

　特異値分解の結果と一致しているか確かめます。Pythonで特異値分解を行うにはNumPyのlinalg.svgを使うとよいでしょう。
以下のように4×3行列を作成します。

A = np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])

power_method関数を使うと、

u, s, v = power_method(A, iter_num=1)

u
#-> array([-0.14091245, -0.34396479, -0.54701712, -0.75006945])
s
#-> 25.462398134435947
v
#-> array([-0.50388184, -0.57446659, -0.64505134])

NumPyで特異値分解を行うと

U, S, V = np.linalg.svd(A, full_matrices=False)

U
#-> array([[-0.14087668,  0.82471435,  0.54771037],
#         [-0.34394629,  0.42626394, -0.72755711],
#         [-0.54701591,  0.02781353, -0.1880169 ],
#         [-0.75008553, -0.37063688,  0.36786364]])
S
#-> array([  2.54624074e+01,   1.29066168e+00,   1.97614008e-15])
V
#-> array([[-0.50453315, -0.5745157 , -0.64449826],
#         [-0.76077568, -0.05714052,  0.64649464],
#         [-0.40824829,  0.81649658, -0.40824829]])

両者を比較すると、左特異ベクトル、特異値、右特異ベクトルがほぼ一致していることがわかります。

2017-10-11

特異スペクトル変換法による時系列データの異常検知（Python）

Python データサイエンス異常検知

はじめに

　今回は、特異スペクトル変換法というアルゴリズムをPythonで実装します。このアルゴリズムは時系列データの異常検知に対して非常に強い力を発揮します。また、ハイパーパラメータ（人が調整する必要のあるパラメータ）が少なく、比較的チューニングが容易であることも特徴の一つです。数学の理論については深追いはせず、アルゴリズムの概要と実装まで書いていきたいと思います。

【目次】

はじめに
時系列データについて
時系列データの異常と変化点検知　
特異スペクトル変換法の概要
Pythonによる実装
特異スペクトル変換法の課題
おわりに

時系列データについて

　時系列データとは、時間の推移ととともに観測されるデータのことです。昨今、様々な企業がデータ活用を推進していますが、世の中の実務の現場に貯まっていく多くのデータは時系列のデータです。
　データ分析において、時系列データの取り扱いは、他のデータ（クロスセクションデータなどと言われる）に比べて、理論的にも実銭的にも難しいです。例えば、分析をする上で、ヒストグラムを描くということをよくやると思いますが、ヒストグラムを描くということは、暗にデータの相関を無視し、互いに独立であることを仮定しています。時系列データでは、隣り合う時刻の観測値同士には明らかに相関関係があり、「各観測値が互いに独立である」という仮定は成り立ちません。このようなデータをきちんと取り扱えることは非常に有用です。

【参考】実務の現場に多い時系列データ分析の際に注意しておきたい点を列挙してみる - 六本木で働くデータサイエンティストのブログ

時系列データの異常と変化点検知　

　一口に異常と言ってもいろいろありますが、明らかに他の大多数のデータから外れた値を持っている場合は非常に検知が容易です。このような異常を検知することは外れ値検知と呼ばれ、正規分布を用いた単純な外れ値検知の話を以前記事に書きました。

blog.tsurubee.tech

一方、時系列データの異常パターンの特徴としては、各観測データの時刻（横軸）を適当にシャッフルすると検知できなくなるということです。実際に例を見てみましょう。

f:id:hirotsuru314:20171011215303p:plain

　データの真ん中あたりで波が激しくなって、また戻っているのがわかると思います。変化点という観点で言うと、波の周波数が激しくなっている点（横軸が30のあたり）と、またゆるやかに戻っている点（横軸が60のあたり）が変化点です。このように人間の目でデータを見ると、どこが異常かだいたいわかりますが、コンピュータにこの異常を検知させるのはけっこう難しいです。これを自動で検知してやろうというのが今回の目的です。今回は、1変数の時系列データにおける異常検知の話です。

特異スペクトル変換法の概要

　時系列データの異常検知の手法はいろいろありますが、今回は特異スペクトル変換というアルゴリズムを実装してみます。特異スペクトル変換法の概略は下の図に集約されます。

f:id:hirotsuru314:20171011214314p:plain

（出典：上の図は井手剛氏の著書「入門機械学習による異常検知―Rによる実践ガイド」のP200 図7.4を元に作成しました。）

履歴行列とテスト行列

　一次元時系列データの分析をする際に、ある任意の幅 $w$ を設定して、 $w$ 個の隣接した観測データをまとめて $w$ 次元のベクトルとして取り出すということをよくやります。さらに $w$ を時間の経過とともにずらしながら取り出すことで、時系列片をつくっていきます。この取り出す領域のことをスライド窓と呼び、スライド窓により取り出したベクトルを部分時系列と呼びます。
　上の概略図からわかるように、特異スペクトル変換法ではこの部分時系列をさらにまとめて行列としてデータを取り出します。この行列自体の特徴が、ある時刻の時系列データの特徴となるわけです。ある時刻 $t$ のまわりデータを用いて生成した行列をテスト行列と呼び、時刻 $t$ より以前のデータを用いて生成した行列を履歴行列と呼びます。特異スペクトル変換法では、この二つの行列の食い違いの大きさを変化度として異常検知を行います。

特異値分解

　特異スペクトル変換法の理論の理解には線形代数の知識がある程度必要なので、場合によっては「特異スペクトル変換法は現在と過去の行列を取り出して、その食い違いの大きさをもとに変化度を計算しているのだな」くらいの理解でよいかもしれません。（後述のPythonによる実装まで読み飛ばしていただいてかまいません。）
　特異値分解は、固有値分解を長方形行列に拡張したものと見なすことができます。詳細な説明については下記の参考リンクがわかりやすいと思います。固有値分解ってなに？という方にはプログラミングのための線形代数という本が超絶おすすめです。

特異値分解の定義，性質，具体例 | 高校数学の美しい物語

http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-kisosuri/singularvalue050129.pdf

　数式を追うことも大事ですが、今回やっていることをざっと書くと、前述の履歴行列とテスト行列に特異値分解を行い、特異値上位の左特異ベクトルを取り出し行列としてまとめてます。（行列のランクより少ない数を取り出す）これはマイナーな方向を無視し、ノイズ除去または次元削除をしたことに相当すると考えることができます。また、取り出したそれぞれの列ベクトルで張られる空間を主部分空間などと呼びます。
　後に出てくるコードを見ればわかりますが、いかにも難しそうな特異値分解もNumPyの力を借りれば、関数を1個呼び出すだけで計算結果が返ってきます。

numpy.linalg.svd — NumPy v1.13 Manual

変化度の定義

　時刻 $t$ における過去と現在の行列の主部分空間から、両者の食い違いを定量化することで、時刻 $t$ の変化度を定義します。これは部分空間同士の距離を求める問題になるのですが、空間同士の距離を求めるというのはすごく理論的にも難しいです。今回は、行列2ノルムが行列の最大特異値に等しいということを利用して、異常度を算出しています。
　まとめると、特異スペクトル変換法は、まずスライド窓を使って、現在と過去のデータを行列として取り出し、特異値分解に基づき行列の特徴パターンを求め、そこから両者の変化度を算出するアルゴリズムです。

Pythonによる実装

　作りたいものは、入力として1次元データのリスト（Pythonの組み込みリストやNumPyのarray）を渡すと、出力として特異スペクトル変換法に基づき変化度に変換した1次元データリストを返すというものです。
※今回処理の流れを掴みやすくするため、あえて2つの関数に詰め込んで実装しています。

import numpy as np

def extract_matrix(data, start, end, w):
    row = w
    column = end - start + 1
    matrix = np.empty((row, column))
    i = 0
    for t in range(start, end+1):
        matrix[:, i] = data[t-1:t-1+row]
        i += 1

    return matrix

def sst(data, w, m=2, k=None, L=None):
    """
    Parameters
    ----------
    data : array_like
           Input array or object that can be converted to an array.
    w    : int
           Window size
    m    : int
           Number of basis vectors
    k    : int
           Number of columns for the trajectory and test matrices
    L    : int
           Lag time

    Returns
    -------
    Numpy array contains the degree of change.
    """
    # Set variables
    if not isinstance(data, np.ndarray):
        data = np.array(data)
    if k is None:
        k = w // 2
    if L is None:
        L = k // 2
    T = len(data)

    # Calculation range
    start_cal = k + w
    end_cal = T - L + 1

    # Calculate the degree of change
    change_scores = np.zeros(len(data))
    for t in range(start_cal, end_cal + 1):
        # Trajectory matrix
        start_tra = t - w - k + 1
        end_tra = t - w
        tra_matrix = extract_matrix(data, start_tra, end_tra, w)

        # Test matrix
        start_test = start_tra + L
        end_test = end_tra + L
        test_matrix = extract_matrix(data, start_test, end_test, w)

        # Singular value decomposition(SVD)
        U_tra, _, _  = np.linalg.svd(tra_matrix, full_matrices=False)
        U_test, _, _ = np.linalg.svd(test_matrix, full_matrices=False)
        U_tra_m  = U_tra[:, :m]
        U_test_m = U_test[:, :m]
        s = np.linalg.svd(np.dot(U_tra_m.T, U_test_m), full_matrices=False, compute_uv=False)
        change_scores[t] = 1 - s[0]

    return change_scores

実際にこのコードを使って、先に挙げた周波数異常のデータを変化度に変換し、描画した結果が下記になります。（縦軸は見やすいようにスケールを調整しています。）

change_scores = sst(data, 50)

f:id:hirotsuru314:20171011214524p:plain

変化度を示している赤のラインを見ると周波数変化が起こっている2点で、きわめて明瞭にピークをもつことがわかります。このように、波の周波数が激しくなっている点も、ゆるやかに戻っている点のどちらも検知できています。若干遅れて検知しているのは、パラメータの調整で改善できる部分でもありますが、ある程度限界もあります。
　今回の例は周波数の異常パターンでしたが、同アルゴリズムを使って、データの傾きの変化点なども検知することができます。

特異スペクトル変換法の課題

　特異スペクトル変換法の課題のひとつとして、計算コストが高いことが挙げられます。一連の処理の中でも特に特異値分解の部分の計算量がネックになります。計算が重く時間を要すると言うことは、データのリアルタイム処理には使えないということであり、時系列解析の幅を制限する要因にもなり得ます。
　この課題を解決する方法の一つとして、クリロフ部分空間法という手法と併用することで計算の高速化することができるようです。また、そもそもリアルタイム処理には別のアルゴリズムを使うという選択肢もあるでしょう。例えば、ChangeFinderは逐次的な計算アルゴリズムになっているため、リアルタイム処理には向いています。（一方、ハイパーパラメータの闇は深そうなので、要は適材適所を見極める必要があるということですね。）

おわりに

　現在、今回記事に書いた特異スペクトル変換法による異常検知も実装した『Banpei』という異常検知パッケージを絶賛開発中です。
　 github.com

　今後、様々なアルゴリズムを実装したり、ストリーミングデータの取り扱いができるように機能拡張したり、などいろいろ計画中です。
　私自身エンジニア歴がまだ1年にも満たない未熟者なので、書き方がおかしいところなどあればぜひご指導をよろしくお願いいたしますorz

2017-09-28

プログラマのための数学勉強会@福岡 #6に登壇しました

異常検知データサイエンス

先日「プログラマのための数学勉強会@福岡 #6」に登壇しました。

maths4pg-fuk.connpass.com

こちらがそのときの発表スライドになります。

speakerdeck.com

内容は、正規分布を使って1次元データの異常検知をする話で、理論・実装の詳細は下記の記事に書きました。

blog.tsurubee.tech

このイベントは、@tkengoさんが主催されており、@tkengoさんは最近機械学習の数学の本も出版されたとてもすごい方です！

やさしく学ぶ機械学習を理解するための数学のきほん ~アヤノ&ミオと一緒に学ぶ機械学習の理論と数学、実装まで~

作者: LINE Fukuoka株式会社立石賢吾
出版社/メーカー: マイナビ出版
発売日: 2017/09/21
メディア: 単行本（ソフトカバー）
この商品を含むブログ (2件) を見る

私のような末端エンジニアが@tkengoさんから登壇のお誘いをいただき大変恐縮でしたが、発表準備の過程で、私自身がかなり勉強させていただきました。また、発表のおかげで、今回のテーマである「異常検知」についてかなり興味を持ち、継続して勉強中なので、そのあたりはまたアウトプットできればと思います！

2017-09-21

1次元正規分布に基づく異常検知の理論とPythonによる実装

異常検知 Python データサイエンス

はじめに

　異常検知とは、大多数のデータとは振る舞いが異なるデータを検出する技術のことです。異常検知は、膨大なデータが収集可能となった現代におけるデータ活用のひとつとして脚光を浴びています。

統計的異常検知の考え方

　異常検知にもいろいろアプローチがあります。例えば、距離に基づくアプローチ、統計学に基づくアプローチ、機械学習に基づくアプローチなどが挙げられます。今回はその中でも統計学に基づくアプローチについて書きます。統計学に基づく異常検知のアプローチは統計的異常検知とも呼ばれます。
　統計的異常検知とは、データの生成機構が確率モデル（確率分布）で表現できると仮定した場合の異常検知の方法論です。これはデータから得られる「知識」を確率モデルという形で表現していると捉えることもできます。

統計的異常検知の手順

統計的異常検知は基本的に以下の3ステップからなります。

観測データからデータ生成の確率モデルを学習
さらに細かく分けると以下の２ステップあります。
① 未知パラメータを含む確率分布を仮定
② データから未知パラメータを推定
学習したモデルを基に、データの異常度合いをスコアリング（異常度の定義）
閾値の設定

今回は、統計的異常検知の中でも最も基本的な、1次元正規分布に基づく異常検知について、上記の３ステップを一つずつ解説していきます。

データの準備

用いるデータは、Rのデータセットにある「Davisデータ」というもので、200人の性別、測定した身長・体重、および自己申告の身長・体重が記録されています。 https://vincentarelbundock.github.io/Rdatasets/datasets.html

　前提として、今回扱うデータは測定した体重だけであり、１変数です。すなわち確率的に揺らぐ値（確率変数）が1種類だけ存在するということです。もちろん、これから示す方法論は、多変量にも拡張できますが、1変数から始めた方が数式が追いやすいです。
　横軸を標本番号、縦軸を体重として、データを可視化したものが以下になります。

f:id:hirotsuru314:20170919224410p:plain:w550

　また、このデータのヒストグラムは下のようになります。

f:id:hirotsuru314:20170919224425p:plain:w550

　どちらもパッと見で異常値が見つけられると思います。このデータの異常判定をしたいとき、「120kg以上は異常だ！」としてif文で検知するのは簡単ですが、その120kgという値には何の根拠もなく、説得力に欠けると思います。これを統計学を駆使して、より客観的な水準で判定してやろうというのが今回の話です。
観測データがN個あるとき、データをまとめて $D$ という記号で表すとします。

$\displaystyle D= \{x_1,x_2,...,x_N\}$

ここで、観測データが多次元の場合は、 $x$ が列ベクトルとして表されます。また、このデータの中には異常な観測データが含まれていないか、含まれていたとしてもその影響は無視できると仮定します。
これから前述した統計的異常検知の3ステップに沿って話を進めていきます。

1. 観測データからデータ生成の確率モデルを学習する

① 未知パラメータを含む確率分布を仮定
　上記のヒストグラムによると、体重データの分布は若干左右非対称ではありますが、おおむね山の形になっています。そこで、それぞれの観測データ $x$ が平均 $\mu$ 、分散 $\sigma^2$ の正規分布に従うと仮定します。

$\displaystyle N(x|\mu, \sigma^2)=\frac{1}{\sqrt{2\pi \sigma^2}}\exp{\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}}$

② データから未知パラメータを推定
　正規分布に含まれているパラメータは平均 $\mu$ 、分散 $\sigma^2$ の二つです。言い換えると、この二つのパラメータが決定すると、分布の形が一意に決まるということです。これらのパラメータを観測データから推定します。パラメータ推定には最尤推定を用います。最尤推定は、「観測データが得られる確率」をパラメータの関数とみなした尤度関数を最大化するようにパラメータを決定する推定法です。

$\displaystyle \begin{eqnarray} p(D|\mu, \sigma^2) & = & \prod_{n=1}^N N(x_n|\mu, \sigma^2) \\ & = & \prod_{n=1}^N\frac{1}{\sqrt{2\pi \sigma^2}}\exp{\left\{-\frac{(x_n-\mu)^2}{2\sigma^2}\right\}} \\ & = & (2\pi \sigma^2)^{-\frac{N}{2}}\exp{\left\{-\frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n-\mu)^2 \right\}} \end{eqnarray}$

ここで、計算の都合上、尤度関数の自然対数をとった対数尤度関数を最大化します。対数は単調増加であるため、尤度関数の最大化は対数尤度関数の最大化と同値です。
$\beta=\frac{1}{\sigma^2}$ とした場合、対数尤度尤度関数は、

$\displaystyle \ln p(D|\mu, \sigma^2)=\frac{N}{2}\ln \beta - \frac{N}{2}\ln 2\pi - \frac{\beta}{2}\sum_{n=1}^{N}(x_n-\mu)^2$

これを最大化するパラメータを求めるため、 $\mu$ と $\beta$ でそれぞれ偏微分してゼロと等置します。

$\displaystyle \frac{\partial}{\partial \mu}\ln p(D|\mu, \sigma^2)=0 \tag{1}$

$\displaystyle \frac{\partial}{\partial \beta}\ln p(D|\mu, \sigma^2)=0 \tag{2}$

$(1)$ より

$\displaystyle \begin{eqnarray} \frac{\partial}{\partial \mu}\ln p(D|\mu, \sigma^2) & = & \beta\sum_{n=1}^{N}(x_n-\mu) \\ & = & \beta \biggl(\sum_{n=1}^{N}x_n-N\mu\biggl) \\ \end{eqnarray}$

$\displaystyle ∴　\mu=\frac{1}{N}\sum_{n=1}^{N}x_n$

$(2)$ より

$\displaystyle \frac{\partial}{\partial \beta}\ln p(D|\mu, \sigma^2) = \frac{N}{2\beta}-\frac{1}{2}\sum_{n=1}^{N}(x_n-\mu)^2$

$\displaystyle ∴　\sigma^2=\frac{1}{N}\sum_{n=1}^{N}(x_n-\mu)^2$

　平均・分散の推定値は私たちがよく知っている平均・分散の定義そのものだということがわかると思います。これらの結果は、平均・分散の推定値として、観測データの標本平均・標本分散を採用することを意味しています。また、推定した平均・分散にはデータからの推定値であることを明示するため「^（ハット）」をつけます。
したがって、学習モデル（予測分布） $p(x|\hat{\mu}, \hat{\sigma^2})$ が得られたことになります。
ここで、 $\hat{\mu}$ と $\hat{\sigma^2}$ はデータセットが与えられると一意に決まります。（普通にデータの平均・分散を出すだけです。）

2. 学習したモデルを基に、データの異常度合いをスコアリング（異常度の定義）

　一般に、異常度の定義として、負の対数尤度を採用されることが多いです。（これは情報理論におけるシャノン情報量と呼ばれるものです）
新たな観測値 $x'$ に対する異常度 $a(x')$ は、

$\displaystyle \begin{eqnarray} a(x') & = & -\ln p(x'|\hat{\mu}, \hat{\sigma^2}) \\ & = & -\ln \biggl[\frac{1}{\sqrt{2\pi \hat{\sigma}^2}}\exp{\left\{-\frac{(x'-\hat{\mu})^2}{2\hat{\sigma}^2}\right\}} \biggr] \\ & = & \frac{1}{2}\ln (2\pi \hat{\sigma}^2)+\frac{1}{2\hat{\sigma}^2}(x'-\hat{\mu})^2 \end{eqnarray}$

第1項は観測データ $x'$ に依存しないので無視します。全体に2をかけると、異常度は、

$\displaystyle a(x')=\biggl(\frac{x'-\hat{\mu}}{\hat{\sigma}}\biggl)^2$

という形で表されます。
分子 $(x'-\hat{\mu})$ は観測データの標本平均からのずれの大きさを表しており、ずれが大きくなると異常度が大きくなります。
分母は $\hat{\sigma}$ であり、標準偏差で割ることで、もともとばらつきの大きいものは多少のずれでも多めに見るようになっており、異常度の表現として非常に直感に近い形になっていることがわかります。

3. 閾値の設定

　異常度が決まると、それに閾値を設定することで異常判定をすることができます。閾値は経験に基づき主観的に決める方法もありますが、できるだけ客観的基準に基づいて決めるのが望ましいです。
　ホテリング理論は、観測データが正規分布に従うことを仮定した下で異常検知を行う非常に有名な古典理論であり、広く応用されています。ホテリング理論が有効な理由の一つは、異常度が従う確率分布を、明示的に導くことができる点です。異常度の確率分布が既知であれば、それに基づいて閾値を設定するのは容易です。例えば、「正常には1パーセント未満でしか起こらないくらい稀な値だから、きっと正常ではないだろう」という論理になります。つまり、パーセント値により異常判定を行うことができるようになります。
　ホテリング理論によると異常度 $a(x’)$ はデータ数Nが十分に大きい時、自由度1のカイ二乗分布に従うということが数学的に証明できます。
　下のグラフは自由度1のカイ二乗分布の確率密度関数です。

f:id:hirotsuru314:20170919224519p:plain:w550

　ここで、カイ二乗分布とは、 $Z_1,Z_2,...Z_n$ を独立な標準正規分布 $N(0, 1)$ に従う確率変数としたとき、

$\displaystyle \chi^2=Z_1^2+Z_2^2+...+Z_n^2$

とすると、確率変数 $\chi^2$ が従う確率分布を自由度 $n$ のカイ二乗分布といいます。
（異常度がカイ二乗分布に従うことの数学的な証明は難解なため避けますが、Nが十分に大きいとき、標本平均・標本分散は母平均・母分散に近づき、 $a(x’)$ が「確率変数から母平均を引いて、母標準偏差で割ったもの」すなわち標準正規分布 $N(0, 1)$ の2乗に従うようになるからだと考えてよいのだと思います）

Pythonによる実践

　これまで説明した1次元正規分布に基づく異常検知をPythonで実装したいと思います。
　ここで、入力として、データと異常度の閾値(%)を与えると、出力として、異常値とそのインデックス番号を返す関数を作りたいと思います。

import numpy as np
from scipy import stat

def hotelling_1d(data, threshold):
    """
    Parameters
    ----------
    data : Numpy array
    threshold : float

    Returns
    -------
    List of tuples where each tuple contains index number and anomalous value.
    """
    #Covert raw data into the degree of abnormality
    avg = np.average(data)
    var = np.var(data)
    data_abn = [(x - avg)**2 / var for x in data]

    #Set the threshold of abnormality
    abn_th = stats.chi2.interval(1-threshold, 1)[1]

    #Abnormality determination
    result = []
    for (index, x) in enumerate(data_abn):
        if x > abn_th:
            result.append((index, data[index]))
    return result

　この関数に引数のdataとして、200人の体重データを渡すと以下のように、

hotelling_1d(data, 0.01)
#-> [(11, 166), (20, 119)]

　異常値として2点返ってきます。
　結果を可視化すると下のようになります。縦軸は異常度です。グラフ内の赤の破線は異常度の閾値です。（プログラム中の「abn_th」の値に対応）

f:id:hirotsuru314:20170919224604p:plain:w550

　プログラムを見てもわかりますが、結局ほぼデータの平均と分散を計算して、ちょろっと異常度に変換するとできてしまいます。カイ二乗分布の確率密度関数は非常に複雑ですが、SciPyを用いるとその積分値は簡単に出せます。

ホテリング理論の課題

　ホテリング理論では、観測データ同士が互いに独立で、それぞれが単一の正規分布に従っていると仮定していますが、それは非常に強い制約と言えます。したがって、下記のようなデータには適用は困難です。

１. 複数のモードがあるような複雑な系への適用
　例えば、複数のピークを持ち、単一の正規分布で仮定できないような場合です。

２. 値が動的に変化する時系列データへの適用が困難
　正規分布を仮定するということは、データの分布がひと山でだいたい安定しているとみなしいますが、平均の値が揺らぐようなデータもあります。いわゆる時系列データと呼ばれる類のものですが、そのようなデータにホテリング理論は適用が困難です。
　時系列データの異常検知についてはまた別の機会に書きたいと思います。

参考文献

入門機械学習による異常検知―Rによる実践ガイド

作者: 井手剛
出版社/メーカー: コロナ社
発売日: 2015/02/19
メディア: 単行本
この商品を含むブログ (4件) を見る

現在日本語で読める異常検知の本ではおそらく一番やさしく、読みやすい。時系列データの異常検知についても載っている。

2017-09-02

Python製の対話的可視化ライブラリBokehを使ってみる

Python データサイエンス

　Pythonには様々なデータ可視化ライブラリがありますが、私は最近Bokehというライブラリを知って、その便利さにハマってます！今回はBokehの簡単なチュートリアル的な内容を書きたいと思います。

Bokehとは？

　Bokehって何？の答えを知るには下記の公式ページを見るのが一番良いと思います。ただ、日本語版はないようです。

bokeh.pydata.org

　Bokehの特徴をいくつか挙げると、

対話的な可視化ができる（対話的って何？というのは後で書きます）
斬新奇抜なグラフが描ける
ビッグデータをリアルタイムで処理できる
グラフを作るためにHTML・CSS・JavaScriptを書く必要がない

といった感じです。
ちなみにBokehの主な開発者は、有名なPythonパッケージであるAnacondaを開発している人たちのようです。
実際の開発者によるBokehについてのプレゼン資料がslideshareにあり、非常に面白いです。

Hassle Free Data Science Apps with Bokeh Webinar

Bokehを触って見る

　それでは実際にBokehを触ってみようと思います。

#Importing libraries
from bokeh.plotting import figure
from bokeh.io import output_file, show

#prepare data
x = [1,2,3,4,5]
y = [6,7,8,9,10]

#prepare the output file
output_file("Line.html")

#create figure
fig = figure()
fig.line(x, y)
show(fig)

こんな短いコードで対話的なグラフが描けます。今回データについてはコードに直書きしています。こちらを実行すると（コマンドラインでpython file名を叩くと）、ブラウザ上にページが立ち上がり、下のようにグラフが描画されています。

f:id:hirotsuru314:20170902100536p:plain

これを右上にドラッグすると、

f:id:hirotsuru314:20170902100545p:plain

スクロールすると、

f:id:hirotsuru314:20170902100559p:plain

対話的とはこういうことで、ユーザーの操作に対して即座に動きが返ってくる仕組みになっています。ぜひ実際にやってみてグラフをグリグリできる楽しさを味わってみてください！（動画を載せようと思ったのですが、パッと調べた感じYouTube等を使う方法しかなく、画像にしました。わかりづらくてすみません、いい方法あれば教えてください^^;）

上記のコードの中の

output_file("Line.html")

の部分で、Bokehが生成してくれたHTMLファイルも出力するようにしているので、出力されたHTMLの中をみてみましょう。

<!DOCTYPE html>
<html lang="en">
    <head>
        <meta charset="utf-8">
        <title>Bokeh Plot</title>

<link rel="stylesheet" href="https://cdn.pydata.org/bokeh/release/bokeh-0.12.2.min.css" type="text/css" />

<script type="text/javascript" src="https://cdn.pydata.org/bokeh/release/bokeh-0.12.2.min.js"></script>
<script type="text/javascript">
    Bokeh.set_log_level("info");
</script>
        <style>
          html {
            width: 100%;
            height: 100%;
          }
          body {
            width: 90%;
            height: 100%;
            margin: auto;
          }
        </style>
    </head>
    <body>

        <div class="bk-root">
            <div class="plotdiv" id="a7161ecc-bd75-43de-9089-247eb36da918"></div>
        </div>

        <script type="text/javascript">
            Bokeh.$(function() {
            Bokeh.safely(function() {
                var docs_json = {"ceb05a6f-303a-42ec-8fe1-717a3c3ba468":{"roots":{"references":[{"attributes":{"overlay":{"id":"c8df6c4b-20d9-4d00-ab91-8696851c6517","type":"BoxAnnotation"},"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"fe131afa-6cb9-4027-8e5e-81f50a7f7c2f","type":"BoxZoomTool"},{"attributes":{"formatter":{"id":"e9a2431f-b252-4bc6-86ca-98dd06d792ca","type":"BasicTickFormatter"},"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"},"ticker":{"id":"1e5e43c3-495b-471d-a871-7f085567f5fe","type":"BasicTicker"}},"id":"272fbeb6-e972-4283-94f3-c5c01022aaf3","type":"LinearAxis"},{"attributes":{"callback":null},"id":"50aa3059-7d61-402f-96d9-e51c4d7b4ab4","type":"DataRange1d"},{"attributes":{"bottom_units":"screen","fill_alpha":{"value":0.5},"fill_color":{"value":"lightgrey"},"left_units":"screen","level":"overlay","line_alpha":{"value":1.0},"line_color":{"value":"black"},"line_dash":[4,4],"line_width":{"value":2},"plot":null,"render_mode":"css","right_units":"screen","top_units":"screen"},"id":"c8df6c4b-20d9-4d00-ab91-8696851c6517","type":"BoxAnnotation"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"},"ticker":{"id":"6ccb7a55-c184-446a-b680-3135cd27775d","type":"BasicTicker"}},"id":"1844c54b-118b-452b-83bd-ed647412592f","type":"Grid"},{"attributes":{"data_source":{"id":"a13993f5-d7b0-4174-88bb-b90f2e6c1412","type":"ColumnDataSource"},"glyph":{"id":"bafa0919-8be4-4742-8c66-2805ce1f9da4","type":"Line"},"hover_glyph":null,"nonselection_glyph":{"id":"29ac1919-026e-4ed7-a866-b13a07c96bfe","type":"Line"},"selection_glyph":null},"id":"a7d00cc9-6e84-413a-a3e9-e1eb94796522","type":"GlyphRenderer"},{"attributes":{"formatter":{"id":"f521c266-6f21-4f57-8c1d-5d8196ce5fa2","type":"BasicTickFormatter"},"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"},"ticker":{"id":"6ccb7a55-c184-446a-b680-3135cd27775d","type":"BasicTicker"}},"id":"449bd4d4-7483-4290-88a0-f1593170ee1d","type":"LinearAxis"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"d64e9c84-783f-4d2c-9254-46a5b09d79e8","type":"PanTool"},{"attributes":{"callback":null},"id":"963839b2-34d1-414c-be56-c84bdb1ef5cc","type":"DataRange1d"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"9a747cea-226e-4006-b86a-6d6b186765f1","type":"WheelZoomTool"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"8109f175-be7e-4f0f-b4b3-3d0bc4762177","type":"HelpTool"},{"attributes":{"line_alpha":{"value":0.1},"line_color":{"value":"#1f77b4"},"x":{"field":"x"},"y":{"field":"y"}},"id":"29ac1919-026e-4ed7-a866-b13a07c96bfe","type":"Line"},{"attributes":{},"id":"e9a2431f-b252-4bc6-86ca-98dd06d792ca","type":"BasicTickFormatter"},{"attributes":{"dimension":1,"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"},"ticker":{"id":"1e5e43c3-495b-471d-a871-7f085567f5fe","type":"BasicTicker"}},"id":"386d29c4-a962-48fe-815f-18666d39826c","type":"Grid"},{"attributes":{},"id":"f9cb966f-a936-4979-9465-afee14f8f21d","type":"ToolEvents"},{"attributes":{},"id":"6ccb7a55-c184-446a-b680-3135cd27775d","type":"BasicTicker"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"2de10c54-d536-4ec1-9904-d39d9c57d1cb","type":"SaveTool"},{"attributes":{"callback":null,"column_names":["x","y"],"data":{"x":[1,2,3,4,5],"y":[6,7,8,9,10]}},"id":"a13993f5-d7b0-4174-88bb-b90f2e6c1412","type":"ColumnDataSource"},{"attributes":{"active_drag":"auto","active_scroll":"auto","active_tap":"auto","tools":[{"id":"d64e9c84-783f-4d2c-9254-46a5b09d79e8","type":"PanTool"},{"id":"9a747cea-226e-4006-b86a-6d6b186765f1","type":"WheelZoomTool"},{"id":"fe131afa-6cb9-4027-8e5e-81f50a7f7c2f","type":"BoxZoomTool"},{"id":"2de10c54-d536-4ec1-9904-d39d9c57d1cb","type":"SaveTool"},{"id":"5dc6dca4-c75b-428e-a030-60e291cd1162","type":"ResetTool"},{"id":"8109f175-be7e-4f0f-b4b3-3d0bc4762177","type":"HelpTool"}]},"id":"c1e5b658-fb25-422c-947d-a9e05d1babb9","type":"Toolbar"},{"attributes":{"line_color":{"value":"#1f77b4"},"x":{"field":"x"},"y":{"field":"y"}},"id":"bafa0919-8be4-4742-8c66-2805ce1f9da4","type":"Line"},{"attributes":{},"id":"1e5e43c3-495b-471d-a871-7f085567f5fe","type":"BasicTicker"},{"attributes":{},"id":"f521c266-6f21-4f57-8c1d-5d8196ce5fa2","type":"BasicTickFormatter"},{"attributes":{"below":[{"id":"449bd4d4-7483-4290-88a0-f1593170ee1d","type":"LinearAxis"}],"left":[{"id":"272fbeb6-e972-4283-94f3-c5c01022aaf3","type":"LinearAxis"}],"renderers":[{"id":"449bd4d4-7483-4290-88a0-f1593170ee1d","type":"LinearAxis"},{"id":"1844c54b-118b-452b-83bd-ed647412592f","type":"Grid"},{"id":"272fbeb6-e972-4283-94f3-c5c01022aaf3","type":"LinearAxis"},{"id":"386d29c4-a962-48fe-815f-18666d39826c","type":"Grid"},{"id":"c8df6c4b-20d9-4d00-ab91-8696851c6517","type":"BoxAnnotation"},{"id":"a7d00cc9-6e84-413a-a3e9-e1eb94796522","type":"GlyphRenderer"}],"title":{"id":"57403a00-d3e6-47c5-a4e0-0fb338855f24","type":"Title"},"tool_events":{"id":"f9cb966f-a936-4979-9465-afee14f8f21d","type":"ToolEvents"},"toolbar":{"id":"c1e5b658-fb25-422c-947d-a9e05d1babb9","type":"Toolbar"},"x_range":{"id":"50aa3059-7d61-402f-96d9-e51c4d7b4ab4","type":"DataRange1d"},"y_range":{"id":"963839b2-34d1-414c-be56-c84bdb1ef5cc","type":"DataRange1d"}},"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"},{"attributes":{"plot":{"id":"89f15627-1da2-4be9-8b49-7242d04f73e2","subtype":"Figure","type":"Plot"}},"id":"5dc6dca4-c75b-428e-a030-60e291cd1162","type":"ResetTool"},{"attributes":{"plot":null,"text":null},"id":"57403a00-d3e6-47c5-a4e0-0fb338855f24","type":"Title"}],"root_ids":["89f15627-1da2-4be9-8b49-7242d04f73e2"]},"title":"Bokeh Application","version":"0.12.2"}};
                var render_items = [{"docid":"ceb05a6f-303a-42ec-8fe1-717a3c3ba468","elementid":"a7161ecc-bd75-43de-9089-247eb36da918","modelid":"89f15627-1da2-4be9-8b49-7242d04f73e2"}];

                Bokeh.embed.embed_items(docs_json, render_items);
            });
        });
        </script>
    </body>
</html>

データはJSON形式に変換されて、JavaScriptに渡されています。このようにHTML・CSS・JavaScriptとして出力されるため、Webとの親和性も高く、FlaskやDjangoを使ったアプリへの組み込みも容易にできます。

さいごに

　Bokehの公式ページにあるギャラリーには超ファンシーなグラフたちがあり、実際に触れるようになっています。えっほんとにこんなことまでできるの？っていうものもあるので、ぜひ見てみてください！

https://bokeh.pydata.org/en/latest/docs/gallery.html

2017-08-13

Pythonで状態空間モデルを使う（StatsModels）

Python

今回は、代表的な時系列モデルである状態空間モデルをPythonで使う方法を書いていきます。先日、『時系列データ分析とPython』という記事を書きましたが、今回はその内容の実装部分にあたります。（状態空間モデルって？という方はぜひ前回の記事を見て下さい！）

blog.tsurubee.tech

用いるデータは、時系列のデータセットとして大変有名な『各年ごとのナイル川流量』のデータです。

f:id:hirotsuru314:20170702224936p:plain

このデータは今回用いる統計ライブラリであるStatsModelsの中のデータセットにも入っています。

statsmodels/nile.csv at master · statsmodels/statsmodels · GitHub

さっそくですが、このデータに状態空間モデルをあてはめてみます。

import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline
import statsmodels.api as sm

#日付形式で読み込む
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y') 
#CSV読み込み
data = pd.read_csv('nile.csv', index_col='year', date_parser=dateparse)

# ローカルレベルモデルによる状態推定
model = sm.tsa.UnobservedComponents(data, 'local level')
result = model.fit()
fig = result.plot_components()

f:id:hirotsuru314:20170813224510p:plain

　状態空間モデルの実質的なコードは最後のたった3行だけです。
　黒線のObserved（観測データ）に対して、青線が推定した状態（潜在変数）になります。
　今回あてはめたモデルは状態空間モデルの中でも一番単純なローカルレベルモデルと呼ばれるモデルです。このモデルはランダムウォークプラスノイズモデルとも呼ばれ、ランダムウォークで推移している潜在変数（真の状態）に、何かしらのノイズ（不規則要素）がのっかって、観測データが得られるといったシンプルなモデルです。
　このようにStatsModelsを使うと、とても短いコードで状態空間モデルを使うことができます。下記のリンクがStatsModelsの状態空間モデル部分のリファレンスです。

statsmodels.tsa.statespace.structural.UnobservedComponents — statsmodels 0.9.0 documentation

今回の例では

model = sm.tsa.UnobservedComponents(data, 'local level')

のように第2引数として'local level'を渡していますが、上記リンクのStatsModelsリファレンスに記載のとおり、ここが実際にどのようなモデルを用いるかを決めます。したがって、例えば引数に、‘local linear trend’を渡せばローカルレベルモデルに線形のトレンドを加えたモデルに拡張することができます。また、『seasonal=12』といった引数を渡せば月の季節性成分を入れることも容易にできます。
　実際のモデリングのプロセスでは、データをよく再現するモデルを見つけるため、様々な要素（トレンド成分、季節性成分など）を入れながら、モデルの試行錯誤を行います。その上で参考となる指標も、下記のように簡単にsummaryとして出すことができます。

print(result.summary())

f:id:hirotsuru314:20170813230712p:plain

　推定したパラメータだけでなくLjung-Box検定などの検定の結果もまとめて出すことができます。
　今回の状態推定はStatsModelsのUnobservedComponentsの中のfit関数を呼び出しただけです。公式リファレンスによると、『Fits the model by maximum likelihood via Kalman filter.』と書いてあり、裏側ではカルマンフィルタによる最尤推定を行なっているようです。

参考　

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

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

　本のタイトルは『カルマンフィルタ』とピンポイントですが、状態空間モデルの基礎からカルマンフィルタ、粒子フィルタによる非線形非ガウス状態空間モデルまで幅広く学べる良書です。

2017-07-02

時系列データ分析とPython

Python

　先日、『時系列データ分析とPython』というタイトルでLTをしたので、そのときのスライドをこちらに載せておきます。

時系列データ分析とPython from Hirofumi Tsuruta

www.slideshare.net

　LTで話したとは言っても、私自身、数ヶ月前まで時系列データなんてほとんど触ったことなくて、ここ最近興味を持ち、勉強を始めました。スライドには最小限のことしか載せてないので、こちらに内容の補足を書いていきます。

時系列データの取り扱いは難しい

　時系列データとは、時間の推移ととともに観測されるデータのことです。私たちの身の回りには時系列データが溢れています。例えば、気温・雨量といった気象データも時系列データですし、株価や為替といった金融データも時系列データです。
　これらのデータは時間の推移とともに観測されるというのはもちろんですが、多くの場合で時間依存性を持ちます。例えば、株価の場合、前日の値動きはどうか、上昇傾向なのか、下落傾向なのかが重要だというのはイメージしやすいかと思います。この時間依存性がデータ分析の際にやっかいな存在です。
　多くのデータ分析手法において、「それぞれのデータを独立と仮定する」という仮定を置いた上で分析を行うケースがよくあります。しかし、時間依存性を持っていると、この仮定が成り立たず、取り扱いが非常に難しくなります。
　下のグラフは、ナイル川流量の各年ごとのデータを表しており，有名な時系列データセットです。

f:id:hirotsuru314:20170702224936p:plain

Pythonの統計ライブラリStatsModels

　Pythonで統計モデリングしようと思うと、StatsModelsを使うことになるかと思います。

StatsModels: Statistics in Python — statsmodels 0.9.0 documentation

　StatsModelsは数多くの統計モデリング手法を提供してくれるライブラリです。一般化線形モデル(generalized linear model; GLM)に始まり、今回使った状態空間モデルもStatsModelsを使えば、簡単に構築できます。ちなみに状態空間モデルについてはこちらに書いてあります。

statsmodels.tsa.statespace.structural.UnobservedComponents — statsmodels 0.9.0 documentation

状態空間モデル

　状態空間モデルでは、データの変動の原因を「水準の変化」と「観測誤差」という二つに分けて考えます。イメージ図は下記のような感じです。 f:id:hirotsuru314:20170702225141p:plain

　状態空間モデルを使えば、欠測データや時変回帰係数、さらに多変量への拡張まで状態空間の枠組みの中で、柔軟にモデリングが行えます。状態空間モデルについては、下記リンクと、その関連ページを読むと、だいたい感覚はつかめると思います。

状態空間モデル | Logics of Blue

　今回、状態推定はStatsModelsのUnobservedComponentsの中のfit関数を呼び出しただけです。公式リファレンスによると、『Fits the model by maximum likelihood via Kalman filter.』と書いてあり、裏側ではカルマンフィルタによる最尤推定を行なっているようです。状態空間モデルにおけるパラメータ推定法は下記の記事がよくまとまっていて、参考になります。カルマンフィルタ・粒子フィルタ・マルコフ連鎖モンテカルロ法（Markov chain Monte Carlo methods; MCMC）について整理してあります。

状態空間モデルの推定方法の分類 | Logics of Blue

実装部分についても書きました（2017年8月13日更新）

blog.tsurubee.tech

参考書　

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

作者: 横内大介,青木義充
出版社/メーカー: 技術評論社
発売日: 2014/02/18
メディア: 単行本（ソフトカバー）
この商品を含むブログ (6件) を見る

　おそらく現時点で日本語で読める時系列解析の本ではこの本が一番やさしいと思います。ゼロから始めるには最適だと思います。時系列データとは？に始まり、自己相関・定常過程・単位根過程・ホワイトノイズ・ランダムウォークといった時系列解析でよく出てくる用語を理解するのに非常に役立ちます。実践の部分は、GARCHモデルがメインで、状態空間モデルについては取り扱っていません。サンプルコードはRです。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

　時系列うんぬんの前に統計モデリングやるなら結局みどり本を読まなきゃですね。私もいつもかなり勉強させていただいているTJOさんの記事にもみどり本程度の統計モデリングの知識は必要なスキルセットとしてあげられてました。

データサイエンティストもしくは機械学習エンジニアになるためのスキル要件とは（2017年夏版） - 六本木で働くデータサイエンティストのブログ