2017-11-15

異常検知ライブラリを作ってみた

今回の記事は、前職消防士でゼロからプログラミングを始めた超未熟者の私が、異常検知ライブラリを作った話です。

なぜ作ったか

マインド的背景

消防士を辞めてエンジニアに転職してから1年、いろんな技術に触れました。TensorFlow、scikit-learn、Dockerなどなど、必死になって使い方を覚えました。しかしだんだん、「これ、コマンド覚えてるだけで自分に何も技術身についてなくない？」という疑問や焦りが湧いてきて、自分は便利なツールを使いこなせるようになりたいんじゃなくて、いつの日かみんなが使って便利なツールを作る側になりたいんだ、ということに気づきました。そのような思いから今回初めてライブラリを作り、公開してみました。

データサイエンス的背景

　世の中は時系列データで溢れています。ビジネスの場において、データの何かしらの変化の兆候を捉えて、いち早く意思決定をしたいという場面がよくあります。例えば、観測データが商品の受注データであれば、急に売れなくなった（もしくは売れ出した）変化をいち早く検知し、次の一手を打てると機会損失を防げるかもしれません。私も実際に業務で時系列データの変化を捉えたい場面がありましたが、変化の兆候というのはデータによって様々であり、明らかなものもあれば、そうでないものもありました。様々なタイプの変化に対して柔軟に検知するためには、閾値やルール設定といった力技では難しく、コンピュータに「いつもと違う状態」を自動に検知させることが求められています。これがいわゆる異常検知（ここでは特に変化点検知）の技術です。

作ったもの

ソースコード

Banpei（番兵）
github.com

使い方

インストールは下記のコマンドで行います。（Pythonパッケージ化はしていますが、PyPIには登録していないのでgit cloneしたのちにpip installする必要があります。）

git clone https://github.com/tsurubee/banpei.git
cd banpei
pip install .

一番簡単な使い方は、異常を検知したいデータをPythonのリスト型やNumPy arrayなどの配列にして、detect関数に渡してやると、データの変化度を計算して配列を返します。
下記のコードのように3行のみで異常検知できます。まず、Banpeiをimportし、SSTという特異スペクトル変換法（異常検知アルゴリズムの一つ）を使うためのインスタンスを作成し、detect関数にインプットデータを渡すだけです。インスタンス作成時の引数wはスライド窓の幅であり、チューニングが必要なパラメータです。（詳しくは以前の記事を参照ください）

import banpei
model   = banpei.SST(w=50)
results = model.detect(data)

resultsには変化度が入っており、インプットであるdataと、アウトプットであるresultsを可視化すると、
（data: 緑、results: 赤）※見やすいように縦軸は調整しています。

f:id:hirotsuru314:20171115011620p:plain

このようにデータの兆候が変化する点で、変化度である赤線が立ち上がっているのが見て取れます。この赤の立ち上がりをみることで変化点検知ができそうです。

リアルタイム異常監視への挑戦

もともと私がBanpeiをつくったときは、リアルタイム監視の環境を作りたいというモチベーションがありました。（売り上げのリアルタイム監視やサーバリソースのリアルタイム監視などなど）そこで今回、ブラウザを通してリアルタイム異常監視ができる簡単な仕組みを構築しました。
可視化のツールとしてBokehというPython製の可視化ライブラリを採用しました。簡単なデモ動画を載せておきます。動画の10秒あたりで変化点検知をしています。

特異スペクトル変換による傾きの変化点検知

このようなリアルタイム監視環境を構築したければ、pip install bokehでBokehをインストールして、下のようなスクリプトで可視化できます。（私の環境ではbokeh==0.12.10）

import banpei
import numpy as np
import pandas as pd
from bokeh.io import curdoc
from bokeh.models import ColumnDataSource, DatetimeTickFormatter
from bokeh.plotting import figure
from datetime import datetime
from math import radians
from pytz import timezone

data = []
results = []

def get_new_data():
    global data
    #ここでデータ取得（例：new_data = psutil.cpu_percent()）
    data.append(new_data)

def update_data():
    global results
    get_new_data()
    ret = model.stream_detect(data)
    results.append(ret)
    now = datetime.now(tz=timezone("Asia/Tokyo"))
    new_data = dict(x=[now], y=[data[-1]], ret=[results[-1]])
    source.stream(new_data, rollover=500)

# Create Data Source
source = ColumnDataSource(dict(x=[], y=[], ret=[]))
# Create Banpei instance
model = banpei.SST(w=30)
# Draw a graph
fig = figure(x_axis_type="datetime",
             x_axis_label="Datetime",
             plot_width=950,
             plot_height=650)
fig.title.text = "Realtime monitoring with Banpei"
fig.line(source=source, x='x', y='y', line_width=2, alpha=.85, color='blue', legend='observed data')
fig.line(source=source, x='x', y='ret', line_width=2, alpha=.85, color='red', legend='change-point score')
fig.circle(source=source, x='x', y='y', line_width=2, line_color='blue', color='blue')
fig.legend.location = "top_left"
# Configuration of the axis
format = "%Y-%m-%d-%H-%M-%S"
fig.xaxis.formatter = DatetimeTickFormatter(
    seconds=[format],
    minsec =[format],
    minutes=[format],
    hourmin=[format],
    hours  =[format],
    days   =[format],
    months =[format],
    years  =[format]
)
fig.xaxis.major_label_orientation=radians(90)
# Configuration of the callback
curdoc().add_root(fig)
curdoc().add_periodic_callback(update_data, 1000) #ms

get_new_data関数の部分をカスタマイズして、独自のデータ取得の処理を書けば、bokeh serve file名.pyでlocalhostが立ち上がりブラウザ上で上記の動画のようなものが見れます。Banpei側ではstream_detectという関数を使っていますが、これは渡されたインプットデータのうち、最新の1点の変化度を計算して、数値を返す関数です。

Banpeiの異常検知アルゴリズム

　Banpeiには異常検知アルゴリズムとして特異スペクトル変換を採用しています。概要は以前記事に書きました。

blog.tsurubee.tech

特異スペクトル変換の特徴としては、特定の確率分布を仮定していないため、①多様な変化に頑強に対応できる、②パラメータチューニングが容易、などが挙げられます。パラメータに関しては、スライド窓のサイズwの実質1個のみのチューニングで済みます。（他の必要なパラメータはwから経験的に算出）特異スペクトル変換は愚直に線形計算しているだけ、というところがすごく素敵です。
一方、内部の計算過程で特異分解と呼ばれる行列計算を用いており、これがかなり計算コストが高いです。したがって、検知を行いたい時間間隔や、パラメータwの最適値によっては、リアルタイム計算が厳しくなるかもしれません。
今後、特異スペクトル変換の高速化のために、クリロフ部分空間を用いた下の論文を実装したいと思ってます。

http://epubs.siam.org/doi/pdf/10.1137/1.9781611972771.54

さいごに

　作った感想としては、やっぱりモノ作り楽しいな、です。一旦、画面上で動くものまで作れたのですごく達成感はありました。ここからどのように実際の現場で使えるように持っていくかは真剣に考えています。例えば、前述のアルゴリズムの高速化でよりリアルタイム性を高めたり、メールやSlackに異常通知を送れる機能を実装したりと・・まだまだやりたいことが山積みなので、着実に開発を進めていきたいと思っています。

参考資料

異常検知って？

時系列データって？

ライブラリ（パッケージ）開発するには？

Bokehって？

特異スペクトル変換って？

2017-11-11

Bokehを用いたビットコイン価格のリアルタイム可視化（Python）

Python データサイエンス

　今回はBokehというライブラリを使って、ビットコイン価格のリアルタイム可視化を行う方法について書いていきます。Bokehを使うと、Pythonオンリーで可視化までできるので、非常に便利です。

Bokehとは

　BokehとはPython製の対話的可視化ライブラリです。対話的とはどういうことかというと、作成したグラフのズームアップ・ズームアウト・軸の調整といった様々な操作を手元で簡単に行えるということです。

f:id:hirotsuru314:20171111151354p:plain

これは、Bokeh内部で、HTML・CSS・JavaScriptを自動で生成する仕組みがあることで達成されており、ユーザは基本的にはPythonコードしか書く必要がありません。
Bokehの公式のギャラリーを見るとすごくファンシーな可視化がたくさんあるので見るだけでも楽しいです！（実際に対話的な操作もできます）

Gallery — Bokeh 0.13.0 documentation

※以前の記事に対話的なグラフの作成方法を書きましたので、併せてどうぞ

Bokehの特徴

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

対話的な可視化ができる
ストリーミングデータのリアルタイム可視化ができる
グラフを作るためにHTML・CSS・JavaScriptを書かなくてもよい

といった感じです。
今回は『ストリーミングデータのリアルタイム可視化ができる』に着目した記事になります。
対話的というだけでもすごいのですが、動的に流れてくるデータに対しても、Bokehにデータを流すことでリアルタイム可視化ができ、これが本当に便利です！

用いるデータについて

　ストリーミングデータの例として『ビットコイン価格』を取り上げて見ることにします。（私自身はシステムトレードなどは全くやっていません。笑）
　下記リンクの中のLast Tradeという部分を取得します。この部分は、数秒ごとにリロードすると値が変わるので、定期的にスクレイピングすることで、時系列データを得ることができます。

Bitcoincharts | bitFlyer - JPY Summary

f:id:hirotsuru314:20171111151600p:plain

ビットコイン価格のリアルタイム可視化

開発環境
Python 3.6.3
bokeh==0.12.10
beautifulsoup4==4.6.0
requests==2.13.0

作成プログラム

from bokeh.io import curdoc
from bokeh.models import ColumnDataSource, DatetimeTickFormatter
from bokeh.plotting import figure
import requests
from bs4 import BeautifulSoup
from datetime import datetime
from math import radians
from pytz import timezone

URL = "https://bitcoincharts.com/markets/"
MARKET="bitflyerJPY"

def get_data():
    res = requests.get(URL + MARKET + ".html", headers={"User-Agent": "Mozilla/5.0"})
    content = res.content
    soup = BeautifulSoup(content, "html.parser")
    last_trade = int(soup.find_all("p")[0].span.text)
    return last_trade

def update_data():
    now = datetime.now(tz=timezone("Asia/Tokyo"))
    new_data = dict(x=[now], y=[get_data()])
    source.stream(new_data, rollover=200)

#データソースの作成
source = ColumnDataSource(dict(x=[], y=[]))

#グラフの作成
fig = figure(x_axis_type="datetime",
             x_axis_label="Datetime",
             y_axis_label="Last Trade",
             plot_width=800,
             plot_height=600)
fig.title.text = "Bitcoin Charts"
fig.line(source=source, x="x", y="y", line_width=2, alpha=.85, color="blue")
fig.circle(source=source, x="x", y="y", line_width=2, color="blue")

#軸の設定
format = "%Y-%m-%d-%H-%M-%S"
fig.xaxis.formatter = DatetimeTickFormatter(
    seconds=[format],
    minsec =[format],
    minutes=[format],
    hourmin=[format],
    hours  =[format],
    days   =[format],
    months =[format],
    years  =[format]
)
fig.xaxis.major_label_orientation=radians(90)

#コールバックの設定
curdoc().add_root(fig)
curdoc().add_periodic_callback(update_data, 3000) #ms単位

上記プログラムにbitcoin_streaming.pyというファイル名をつけ、コマンドラインから下のように実行すると、私の場合

$ bokeh serve bitcoin_streaming.py
2017-11-11 12:24:32,513 Starting Bokeh server version 0.12.2
2017-11-11 12:24:32,520 Starting Bokeh server on port 5006 with applications at paths ['/bitcoin_streaming']
2017-11-11 12:24:32,520 Starting Bokeh server with process id: 1995

のように表示され、http://localhost:5006/bitcoin_streaming にアクセスすると、下のようにグラフが表示されます。

f:id:hirotsuru314:20171111151837p:plain

実際はリアルタイムで動くグラフなので、動画にしてみました。

Bokehを用いたビットコイン価格のリアルタイム可視化

ちょっとわかりづらいですが、3秒に1回スクレイピングでビットコイン価格を取得して、リアルタイムに可視化しています。

プログラムについての補足

スクレイピング
 ビットコイン価格のスクレイピングにはBeautifulSoupというライブラリを使っています。BeautifulSoupはPython製のHTMLパーサーで、手軽にスクレイピングができるので超便利です！日本語のドキュメントも豊富なので学習コストも少なく済みます。

Beautiful Soupドキュメント — BeautifulSoup Document 3.0 ドキュメント

ColumnarDataSource
BokehにはColumnarDataSourceというデータ構造があります。ColumnarDataSourceは列指向でデータを保持します。
リアルタイム可視化の場合には、このデータ構造にstream関数を使って、データを流し込むとBokehの方でデータを保持から可視化までやってくれます。
ちなみに

def update_data():
    now = datetime.now(tz=timezone("Asia/Tokyo"))
    new_data = dict(x=[now], y=[get_data()])
    source.stream(new_data, rollover=200)

ここのstream関数の引数であるrolloverは、保持するデータ数の上限を決めるもので、その上限がそのまま可視化にも反映されるので、重要です。

コールバック
作成したプログラムの最後の行の

curdoc().add_periodic_callback(update_data, 3000)

部分で、update_data関数を3秒（3000ms）に1回呼び出す設定を書いていて、update_data関数が実行されると、ColumnarDataSourceのインスタンスに新しいデータが追加され、可視化される仕組みになっています。

参考

Hassle Free Data Science Apps with Bokeh Webinar
　Bokehの開発者のスライド
Data Visualization on the Browser with Python and Bokeh | Udemy
　Udemyというオンライン学習プラットフォームの中のBokehを取り上げているコース。かなり参考にさせて頂いています。超おすすめですが、英語版しかありません。

おわりに

　上記のソースコードはデータ取得のget_data関数の部分だけ変えるといろんなリアルタイム可視化ができます。例えば、CPU使用率を監視しようと思うと、こんな感じに入れ替えればできます

import psutil
def get_new_data():
    return psutil.cpu_percent()

　いろいろ楽しいことができそうですね！

ちなみに・・・

　私は最近個人で『Banpei』という異常検知ライブラリを作っており、Bokehと連携させ、リアルタイム異常監視などができるように開発を進めています。リアルタイム異常監視については近日中にまた記事にします！
　異常検知ライブラリBanpeiはGitHubで公開しているので見ていただけると嬉しいです！

github.com

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