Python geopandas + Bokeh で地理情報をプロットしたい
数日前、pandas を利用して地理情報をプロットするという非常によいエントリが翻訳されていた。
上のエントリ、前処理が手間に見えるが pd.read_html や .str アクセサを使えばもっと簡単に書けると思う、、、がそれは本題でない。
pandas で地理情報を扱う場合、geopandas という拡張パッケージを利用すると便利なため、その使い方を書きたい。また、処理を Python で完結させるため、QGIS ではなく Bokeh でプロットしたい。
geopandas のインストール
pip で。
$ pip install geopandas geopy
このエントリでは依存パッケージである shapely、geopy の機能も利用する。shapely は自動的にインストールされるはずだが、geopy については上のように別途インストールが必要。
地理情報の読み込み
以降の操作は IPython Notebook で行う。まずは必要なパッケージをインポートする。
%matplotlib inline import numpy as np import pandas as pd import geopandas as gpd
例として、東京都の地図を読み込んで利用したい。geopandas は、依存パッケージである fiona がサポートする地理情報ファイル形式を読み込むことができる。GitHub に 地球地図日本 のデータから作成した GeoJSON ファイルがあったため、こちらを利用させていただく。
ファイルの読み込みは、geopandas.read_file で行う。データは pandas.DataFrame を継承した geopandas.GeoDataFrame インスタンスとして読み込まれる。
df = gpd.read_file('tokyo.geojson') type(df) # geopandas.geodataframe.GeoDataFrame isinstance(df, pd.DataFrame) # True df.head() # area_en area_ja code geometry # 0 Tokubu 都区部 131211 POLYGON ((139.821051 35.815077, 139.821684 35.... # 1 Tokubu 都区部 131059 POLYGON ((139.760933 35.732206, 139.761002 35.... # 2 Tokubu 都区部 131016 POLYGON ((139.770135 35.705352, 139.770172 35.... # 3 Tokubu 都区部 131067 POLYGON ((139.809714 35.728135, 139.809705 35.... # 4 Tokubu 都区部 131091 (POLYGON ((139.719199 35.641847, 139.719346 35...
GeoDataFrame.plot を使うと、データを matplotlib を利用してプロットできる。
df.plot()
島嶼部があると細部がわからないため、都区部だけをフィルタしてプロットする。
df[df['area_en'] == 'Tokubu'].plot()
地理情報の処理
GeoDataFrame の各行はそれぞれひとつの領域 (ジオメトリオブジェクト) に対応している。各領域のジオメトリオブジェクトの実体は geometry キー中に shapely の Polygon もしくは MultiPolygon インスタンスとして保存されている。
Polygon には、領域各点の経度/緯度からなる tuple がリストとして保存されている。以下の処理で領域各点の経度/緯度をそれぞれ別のリストとして取り出すことができる。この処理は、続けて記載する Bokeh でのプロット時に必要になる。
# 1行目 (1領域目) を取得 s = df.iloc[0, :] s # area_en Tokubu # area_ja 都区部 # code 131211 # geometry POLYGON ((139.821051 35.815077, 139.821684 35.... # ward_en Adachi Ku # ward_ja 足立区 # Name: 0, dtype: object type(s) # <class 'pandas.core.series.Series'> # 領域のジオメトリオブジェクトの型を表示 type(s['geometry']) # <class 'shapely.geometry.polygon.Polygon'> # 領域の経度/緯度を取得 list(s['geometry'].exterior.coords[:5]) # [(139.821051, 35.815077), (139.821684, 35.814887), (139.822599, 35.814509), # (139.8229070000001, 35.81437), (139.822975, 35.814339)] # 経度、緯度別のリストに分割 lons, lats = zip(*list(s['geometry'].exterior.coords)) # 経度を表示 lons[:5] # (139.821051, 139.821684, 139.822599, 139.8229070000001, 139.822975, 139.822977)
地理情報のプロット
このデータを Bokeh でプロットしたい。Bokeh についてはこちらを。
Bokeh で地理情報をプロットするサンプルとして、ギャラリーにコロプレス図を描く例がある。Bokeh では 描画する各領域の x 座標 (経度)、y 座標 (緯度) をリストで指定する。各領域は複数座標からなるため、"各領域の x/y 座標のリスト" のリスト ( 例: [[境界1の点1のx座標, 境界1の点2のx座標, ...], [領域2の点1のx座標, ...], ...] ) を引数として用意する必要がある。
import shapely # "各領域の x 座標 (経度) のリスト" のリスト countries_xs = [] # "各領域の y 座標 (緯度) のリスト" のリスト countries_ys = [] for i, row in df.iterrows(): if row['area_en'] != 'Tokubu': # 都区部のみ描画 continue polygons = row['geometry'] if isinstance(polygons, shapely.geometry.multipolygon.MultiPolygon): # MultiPolygon に含まれる Polygon を取り出す polygons = [x for x in polygons] elif isinstance(polygons, shapely.geometry.polygon.Polygon): polygons = [polygons] else: raise ValueError for p in polygons: # 境界の各点を経度、緯度に分割 lons, lats = zip(*list(p.exterior.coords)) countries_xs.append(lons) countries_ys.append(lats)
作成したデータを Bokeh に渡す。出力先は IPython Notebook にした。
from bokeh import plotting # 出力先を指定 plotting.output_notebook() # figure を作成 p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700) # 各領域を描画 (塗り分けはしない) p.patches(countries_xs, countries_ys, fill_color='white', line_color="black", line_width=0.5) plotting.show(p)
Geocoding の利用
geopandas は geopy を利用して Geocoding (住所から緯度経度への変換) を行うための API geopandas.geocode.geocode を持っている。が、geopandas v0.1.1, geopy v1.10.0 では、これらのインターフェースに不整合があり そのままでは利用できない。詳細と回避策は以下 Stack Overflow を。
from geopandas.geocode import geocode # geopandas.geocode の API をそのまま使うとエラー geocode(['千代田区千代田1-1']) # AttributeError: 'module' object has no attribute 'MapQuest' # 回避策として geopy の以下プロパティを明示的に None にする import geopy geopy.geocoders.MapQuest = None geocode(['千代田区千代田1-1']) # address geometry # 0 1-1 Chiyoda, Chiyoda-ku, Tōkyō-to 100-0001, Japan POINT (139.7539454 35.6838012)
Geocoding を利用したデータの追加
上で作成した地図に、道路交通センサスにある各地点の交通量をプロットしたい。もっと細かいデータが見られるとよいのだが、公開情報ではこれくらいしかなさそうだ。以下のリンクから東京都の箇所別基本表 kasyo13.xls をダウンロードし、利用させていただく。
-
- 出典：「平成22年度 全国道路・街路交通情勢調査（道路交通センサス）一般交通量調査 集計表」（国土交通省）
pdf = pd.read_excel('kasyo13.xls', header=8, index_col=False) pdf.head() # (km) Unnamed: 1 Unnamed: 2 Unnamed: 3 Unnamed: 4 Unnamed: 5 # 0 13110100010 1 1010 東名高速道路 1 14110100090 # 1 13111010010 1 1101 中央自動車道富士吉田線 4 13200400170 # 2 13111010020 1 1101 中央自動車道富士吉田線 6 NaN # 3 13111010030 1 1101 中央自動車道富士吉田線 6 NaN # 4 13111010040 1 1101 中央自動車道富士吉田線 6 NaN # ヘッダがセル結合されており正しく読めないため、連番で上書き pdf.columns = range(len(pdf.columns)) # 計測地点の住所が入力されている列を表示 pdf.iloc[:, 28].head() # 0 東京～東名川崎 # 1 高井戸ＩＣ～世田谷区・東京都境間 # 2 高井戸ＩＣ～世田谷区・東京都境間 # 3 世田谷区・東京都境～調布ＩＣ間 # 4 高井戸ＩＣ～世田谷区・東京都境間 # Name: 28, dtype: object # NaN, 正しい住所でないデータをフィルタ pdf = pdf.dropna(subset=[28]) pdf = pdf[pdf[28].str.startswith('東京都')] # 計測地点住所と 24時間交通量の列のみにフィルタ pdf = pdf.iloc[:, [28, 36]] pdf.columns = ['計測地点', '交通量合計'] pdf.head() # 計測地点 交通量合計 # 22 東京都中央区八丁堀３丁目 107885 # 23 東京都中央区八丁堀１丁目 98482 # 24 東京都中央区新富１丁目 89089 # 25 東京都中央区銀座１丁目 94143 # 26 東京都中央区築地１丁目 93844
上のデータから、計測地点 カラムの住所をもとに、Geocoding で緯度経度を取得する。
codes = geocode(pdf['計測地点']) type(codes) # <class 'geopandas.geodataframe.GeoDataFrame'> codes.head() # address geometry # 22 3 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032... POINT (139.7756944 35.6755445) # 23 1 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032... POINT (139.7771192 35.6774165) # 24 1 Chome Shintomi, Chūō-ku, Tōkyō-to 104-0041, ... POINT (139.7725002 35.6729586) # 25 Ginza Itchome Station, 1 Chome-6 Ginza, Chūō-k... POINT (139.767045 35.67435) # 26 1 Chome Tsukiji, Chūō-ku, Tōkyō-to 104-0045, J... POINT (139.7716809 35.6702089)
取得したデータは上記のような GeoDataFrame となる。geometry カラムに 経度/緯度を示す Point インスタンスが含まれる。この座標と交通量を利用してバブルチャートが描きたい。Bokeh のギャラリーにバブルチャートをプロットする例があるため、それに倣うと、
pdf = pdf.join(codes) # 各点の x 座標 points_x = [] # 各点の y 座標 points_y = [] # 各点の円の半径 points_v = [] for i, row in pdf.iterrows(): points_x.append(row['geometry'].x) points_y.append(row['geometry'].y) points_v.append(row['交通量合計']) # 円の半径を適当に調整 points_v = np.sqrt(np.array(points_v)) / 50000 # Bokeh での描画 plotting.output_notebook() p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700) p.patches(countries_xs, countries_ys, fill_color='white', line_color="black", line_width=0.5) # バブルチャートを描画 p.scatter(points_x, points_y, radius=points_v, fill_alpha=0.1) plotting.show(p)
まとめ
geopandas を使って、地理情報に対する以下のような処理を行う方法を記載した。
- 各種地理情報フォーマットの読み込み/書き込み
- 各領域 (ジオメトリオブジェクト) からの座標 (緯度/経度) の取得
matplotlib,Bokehでのプロット- Geocoding
公式ドキュメント には、上以外に ジオメトリオブジェクトを操作 / 変換する例も記載されている。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (12件) を見る
Python pandas パフォーマンス維持のための 3 つの TIPS
pandas でそこそこ大きいデータを扱う場合、その処理速度が気になってくる。公式ドキュメントではパフォーマンス向上のために Cython や Numba を使う方法を記載している。
Enhancing Performance — pandas 0.16.2 documentation
が、軽く試したいだけなのに わざわざ Cythonや Numba を使うのは手間だし、かといってあまりに遅いのも嫌だ。そんなとき、pandas 本来のパフォーマンスをできるだけ維持するためのポイントを整理したい。
pandas に限らず、パフォーマンス改善の際にはボトルネックの箇所によってとるべき対策は異なる。pandas では速度向上/エッジケース処理のために データの型や条件によって内部で処理を細かく分けており、常にこうすれば速くなる！ という方法を出すのは難しい。以下はこの前提のうえで、内部実装からみて まあほとんどの場合あてはまるだろう、、、という内容をまとめる。
環境構築
環境は EC2 の c4.xlarge インスタンス上に作成する。パフォーマンスを重視する場合、環境構築の時点で以下を行っておく。
numpyを各種 数値計算ライブラリBLAS,LAPACK,ATLASとリンクさせる。numexpr,bottleneckをインストールする。
EC2 の Amazon Linux であれば以下のコマンドでできる。
$ sudo yum -y install blas-devel lapack-devel atlas-devel $ sudo python -m pip install numpy $ sudo python -m pip install numexpr bottleneck
環境構築が正しくできたかを確認する。以降の処理は IPython Notebook 上で行う。
import numpy as np import numpy.distutils.system_info as sysinfo sysinfo.get_info('blas') # {'libraries': ['blas'], 'library_dirs': ['/usr/lib64'], 'language': 'f77'} sysinfo.get_info('lapack') # {'libraries': ['lapack'], 'library_dirs': ['/usr/lib64'], 'language': 'f77'} sysinfo.get_info('atlas') # {'define_macros': [('ATLAS_INFO', '"\\"3.8.4\\""')], # 'include_dirs': ['/usr/include'], # 'language': 'f77', # 'libraries': ['lapack', 'f77blas', 'cblas', 'atlas'], # 'library_dirs': ['/usr/lib64/atlas']} import pandas as pd pd.show_versions() # ... # pandas: 0.16.2 # ... # bottleneck: 1.0.0 # ... # numexpr: 2.4.3 # ...
サンプルデータの作成
サンプルデータ作成には pandas がテスト用に用意しているメソッドを使う。もっとも、似たようなデータが作れれば方法はなんでもいい。
pd.util.testing.rands_array: 指定した文字数 / 要素数のランダムな文字列のnp.ndarrayを作成。pd.util.testing.choice:np.ndarrayからsize個をサンプリング。
# 2 文字 / 3 要素の array を作成 pd.util.testing.rands_array(2, 3) # array(['Xn', 'ZC', 'zj'], dtype=object) # array から 5 回サンプリング pd.util.testing.choice(np.array(['a', 'b']), size=5) # array(['a', 'b', 'b', 'b', 'a'], dtype='|S1')
サンプルデータとして 100 万レコード / 3 列のデータを用意した。
np.random.seed(0) # レコード数 N = 1000000 chars1 = pd.util.testing.rands_array(5, 100) chars2 = pd.util.testing.rands_array(5, 10000) df = pd.DataFrame({'x': np.random.randn(N), 'y': pd.util.testing.choice(chars1, size=N), 'z': pd.util.testing.choice(chars2, size=N)}) df.shape # (1000000, 3) df.head()
y 列は 100 通りの値を、z 列は 10000 通りの値をとる。
chars1[:10] # array(['SV1ad', '7dNjt', 'vYKxg', 'yym6b', 'MNxUy', 'rLzni', 'juZqZ', # 'fpVas', 'JyXZD', 'ttoNG'], dtype=object) len(np.unique(chars1)) # 100 len(np.unique(chars2)) # 10000
1. 行に対するループ / DataFrame.apply は 使わない
pandas.DataFrame は列ごとに異なる型を持つことができる。DataFrame は内部的に 同じ型の列をまとめて np.ndarray として保持している。列ごとに連続したデータを持つことになるため、そもそも行に対するループには向かない。また、DataFrame.iterrows でのループの際には 異なる型を持つ列の値を Series として取り出すため、そのインスタンス化にも時間がかかる。
また 行ごと / 列ごとに 関数を適用するメソッドに DataFrame.apply があるが、このメソッドでは Python の関数を繰り返し呼び出すためのコストがかかる。apply は利便性を重視したメソッドのため、パフォーマンスを気にする場合は避けたほうがよい。
参考 Python pandas データのイテレーションと関数適用、pipe - StatsFragments
上記ふたつの処理の組み合わせである、各行への関数適用 DataFrame.apply(axis=1) について処理時間を %timeit で計測する。まずは 単純に y 列と z 列の値を文字列結合する場合。
def f1(s): return s['y'] + s['z'] %timeit df.apply(f1, axis=1) # 1 loops, best of 3: 13.7 s per loop
この処理は、Series として列を取り出し ベクトルとして行うほうが格段に速い。
%timeit df['y'] + df['z'] # 10 loops, best of 3: 74.9 ms per loop
ただ、関数中に条件分岐を含む場合など、そのままベクトル化できない場合もある。x 列の値によって 結合の順序を変える例を考える。
def f2_1(s): if s['x'] > 0: return s['y'] + s['z'] else: return s['z'] + s['y'] %timeit df.apply(f2_1, axis=1) # 1 loops, best of 3: 16 s per loop
こういった場合は np.vectorize で関数をベクトル化し、引数として各列を Series (もしくは np.ndarray ) として渡したほうが apply よりは速い。
def f2_2(x, y, z): if x > 0: return y + z else: return z + y %timeit pd.Series(np.vectorize(f2_2)(df['x'], df['y'], df['z']), index=df.index) # 1 loops, best of 3: 334 ms per loop
補足 numpy のドキュメント に記載されているとおり、np.vectorize もパフォーマンスを最重視した方法ではない。さらに速くしたい場合は個別にベクトル化した関数を用意する。参考として、Cython で書いた場合の処理時間は以下となる。np.vectorize の 4-5 倍は速いようだ。
%load_ext cython
%%cython import numpy as np from numpy cimport ndarray def f2_3(ndarray[double, ndim=1] x, ndarray[object, ndim=1] y, ndarray[object, ndim=1] z): cdef: int i, length = len(x) double xval object yval, zval ndarray[object, ndim=1] result = np.empty(length, dtype=object) for i in range(length): xval = x[i] yval = y[i] zval = z[i] if xval > 0: result[i] = yval + zval else: result[i] = zval + yval return result
%timeit pd.Series(f2_3(df['x'].values, df['y'].values, df['z'].values), index=df.index) # 10 loops, best of 3: 67.5 ms per loop
補足 Cython については最近以下の書籍が出ていた。翻訳はわからないが、原著はわかりやすく要点がまとまっていたのでおすすめ。
- 作者: Kurt W. Smith,中田秀基,長尾高弘
- 出版社/メーカー: オライリージャパン
- 発売日: 2015/06/19
- メディア: 大型本
- この商品を含むブログ (3件) を見る
2. object 型は使わない
pandas には、グループ化 ( DataFrame.groupby ) や 結合 ( DataFrame.merge ) など、データの値をキーにして行われる処理がある。こういった操作を行う場合は object 型を避けたほうが速くなることが多い。
pandas では文字列を含む列は object 型になることに注意する。numpy の文字列型とならない理由は簡単にいうと欠損値 ( NaN ) の処理のため。例えば文字列で保存された商品コードをキーにしてグループ化 / 結合する場合など、そのままでは object 型として処理されてしまう。
df.dtypes # x float64 # y object # z object # dtype: object
object 型の列をキーとしてグループ化 / 集約したときの処理時間は以下。
%timeit df.groupby('y').mean() # 10 loops, best of 3: 52.3 ms per loop
処理速度をあげるためには、キーとなる値を カテゴリ型 ( pd.Categorical ) に変換しておくとよい。これは R でいう factor にあたる型。
df['y'] = df['y'].astype('category') %timeit df.groupby('y').mean() # 100 loops, best of 3: 6.89 ms per loop
カテゴリ型では、カテゴリのラベル (pd.Categorical.categories ) と ラベルの位置 (pd.Categorical.codes ) を分けて扱う。内部処理は int 型で保存されたラベルの位置について行われるため、object に対する処理と比較すると速くなる。
# カテゴリ型の内部データを表示 df['y'].values # [ewJ6t, JaFfE, ttoNG, F2Sy9, OopuJ, ..., dS9oG, juZqZ, 5gvFn, 9brWJ, 9fxRG] # Length: 1000000 # Categories (100, object): [0dmK0, 1DdJN, 1IZE1, 1m5Qu, ..., x7c5I, ydsVd, yym6b, zYGBj] # カテゴリのラベル df['y'].cat.categories # Index([u'0dmK0', u'1DdJN', u'1IZE1', u'1m5Qu', u'219tH', u'2aMtU', u'403al', # ... # u'yym6b', u'zYGBj'], # dtype='object') # カテゴリのラベルの位置 df['y'].cat.codes.head() # 0 66 # 1 30 # 2 90 # 3 23 # 4 38 # dtype: int8
列がとりうる値の数が増えた場合も、カテゴリ型にしたほうが速い。
# object 型 %timeit df.groupby('z').mean() # 10 loops, best of 3: 78.4 ms per loop # カテゴリ型 df['z'] = df['z'].astype('category') %timeit df.groupby('z').mean() # 100 loops, best of 3: 17.9 ms per loop
列の型 dtype 以外に処理に影響を与えうるのは、欠損値 ( NaN ) の有無 (無いほうが速い)、object 型内での型の混在 (文字列と数値が混ざっている) など。
3. ユニークでない / ソートされていない index は使わない
最後。index についても、上のとおり object 型は避ける。また、特に理由がない場合 値を ユニークにし、かつソートしておく。上記のサンプルデータ作成時点では、特に index を指定していないため int 型の index が昇順で振られている。
df.index # Int64Index([ 0, 1, 2, 3, 4, 5, 6, 7, # 8, 9, # ... # 999990, 999991, 999992, 999993, 999994, 999995, 999996, 999997, # 999998, 999999], # dtype='int64', length=1000000)
この index で 結合 ( DataFrame.join ) した場合の処理時間は以下 (自身との結合のため意味のない処理だが)。
%timeit df.join(df, rsuffix='right_') # 10 loops, best of 3: 70.6 ms per loop
このとき、 DataFrameに対して前処理 (スライス / サンプリングなど) を行うと、処理に応じて index も並びかわる。例えば サンプリング ( DataFrame.sample ) を行った場合、
df2 = df.sample(n=len(df)) df2.shape # (1000000, 3) df2.index # Int64Index([372149, 955220, 320961, 244572, 254656, 74192, 143279, 246307, # 764579, 96091, # ... # 827030, 681861, 492455, 894210, 758153, 327280, 245717, 952466, # 26440, 532620], # dtype='int64', length=1000000)
このような DataFrame で結合を行うと、上と比較して処理時間が長くなっていることがわかる。
%timeit df2.join(df, rsuffix='right_') # 1 loops, best of 3: 185 ms per loop
index の値がユニークかどうか / ソートされているかは 以下のプロパティで確認できる。このプロパティによって一部の内部処理が分かれる。
# ユニークかどうか df.index.is_unique # True # 昇順でソートされているか df.index.is_monotonic_increasing # True df2.index.is_unique # True df2.index.is_monotonic_increasing # False
まとめ
pandas にその本来のパフォーマンスを発揮させるためには、以下 3 点の処理を避けるとよい。
- 行に対するループ /
DataFrame.applyは 使わない object型は使わない- ユニークでない / ソートされていない
indexは使わない
上記を行った上で処理速度に不満がある場合は Cython や Numba で高速化する。もしくは、Dask を使って並列化を考える。
9/24 追記 Dask についてはこちらを。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (10件) を見る
Chainer で Deep Learning: Bokeh で Live Monitoring したい
概要
Deep Learning の学習には時間がかかるため、進捗が都度 確認できるとうれしい。その際、テキストのログ出力では味気ないので、リアルタイムでプロットを眺めたい。
いくつかの Deep Learning パッケージではそういった機能 (Live Monitoring) が提供されている。
- Undefined Intelligence: Monitoring Experiments in Pylearn2
- Live plotting — Blocks 0.0.1 documentation
- fchollet/hualos · GitHub
同じことを Chainer でやりたい。自分は EC2 を使うことが多いので、リモート環境でも利用できるものがいい。そのため、ここでは Bokeh を使うことにした。
Bokeh とは
Bokeh とは、D3.js を利用したブラウザベースのインタラクティブな可視化を実現するパッケージ。どんなものかは公式の Gallery が充実しているのでそちらを。
補足 R 用のパッケージ {rbokeh} もある。
インストール
pip で。
$ pip install bokeh
準備
環境は EC2 上に作成する。Bokeh からは、ファイル、IPython、Bokeh 組み込みのWebサーバ ( bokeh-server ) 上の画面 の 3種類に対して出力することができるが、今回は Ipython を使うことにする。
補足 もっとも、常に IPython を使うわけではないため、 bokeh-server 上でも描画できるようにしたい。少し試したが、EC2 上に Bokeh を置くと bokeh-server の画面は開くが個々のプロットが表示できなかったため諦めた (ローカルでは問題ない)。できた方いたらやり方教えてください。
bokeh-server はプロットの描画画面だけでなく、動的なデータ更新のためのAPIも提供している。この例では 学習の進捗をリアルタイムで描画するために bokeh-server の機能を利用する。bokeh-server はシェルから以下のコマンドで起動できる。このとき、外部からの接続を受け入れるには自身のホスト名 / IP を引数 ip として指定する。
$ bokeh-server --ip=ec2-xxx-xxx-xxx-xxx.ap-northeast-1.compute.amazonaws.com
補足 EC2 では (特に設定していなければ) グローバル/プライベートの IP が異なるため、IP 指定では正しくデータが更新されないようだ。そのため、パブリックDNSを指定した。
Live Monitoring
以降は IPython Notebook から。Live Monitoring のためのプロットを行うクラスを定義する。
import numpy as np import chainer import bokeh.plotting as plotting from bokeh.models import GlyphRenderer class LiveMonitor(object): def __init__(self, server='chainer', url='http://localhost:5006/', **kwargs): # 出力先に IPython Notebook を指定 plotting.output_notebook(url=url) # トレーニング/テストデータの loss, accuracy を描画する figure を定義 self.train_loss = self._initialize_figure(title='Train loss', color='#FF0000', **kwargs) self.train_acc = self._initialize_figure(title='Train accuracy', color='#0000FF', **kwargs) self.test_loss = self._initialize_figure(title='Test loss', color='#FF0000', **kwargs) self.test_acc = self._initialize_figure(title='Test accuracy', color='#0000FF', **kwargs) # figure を 2 x 2 のグリッド (サブプロット) として配置 self.grid = plotting.gridplot([[self.train_loss, self.test_loss], [self.train_acc, self.test_acc]]) # グリッドを描画 plotting.show(self.grid) def _initialize_figure(self, color=None, line_width=2, title=None, title_text_font_size='9pt', plot_width=380, plot_height=280): """ figure の初期化用のメソッド""" figure = plotting.figure(title=title, title_text_font_size=title_text_font_size, plot_width=plot_width, plot_height=plot_height) # 空のデータで折れ線グラフを作成 x = np.array([]) y = np.array([]) figure.line(x, y, color=color, line_width=line_width) return figure def update(self, train_loss=None, train_accuracy=None, test_loss=None, test_accuracy=None): """ プロットを更新するためのメソッド 指定したキーワード引数に対応する figure が更新される """ self._maybe_update(self.train_loss, train_loss) self._maybe_update(self.train_acc, train_accuracy) self._maybe_update(self.test_loss, test_loss) self._maybe_update(self.test_acc, test_accuracy) def _maybe_update(self, figure, value): """ figure の値を更新するメソッド""" if value is not None: # Variable から np.array に戻す if isinstance(value, chainer.Variable): value = chainer.cuda.to_cpu(value.data) # figure が利用している data_source を取得 renderer = figure.select(dict(type=GlyphRenderer)) ds = renderer[0].data_source # data_source 中の値を更新 y = np.append(ds.data['y'], value) ds.data['y'] = y ds.data['x'] = np.arange(len(y)) # session へ返す (とプロットが更新される) plotting.cursession().store_objects(ds)
このクラスをインスタンスにする。URL としては、bokeh-server の起動時に指定したものを 以下 URL の形式で指定する。
monitor = LiveMonitor(url="http://ec2-xxx-xxx-xxx-xxx.ap-northeast-1.compute.amazonaws.com:5006/")
あとは 学習 / テスト中に、monitor.update に適当な引数を渡せばよい。chainer/examples/mnist/train_mnist.py を例にすると、
# ... 略 optimizer.zero_grads() loss, acc = forward(x_batch, y_batch) loss.backward() optimizer.update() monitor.update(train_loss=loss, train_accuracy=acc) # ... 略 loss, acc = forward(x_batch, y_batch, train=False) monitor.update(test_loss=loss, test_accuracy=acc)
IPython 上 ( LiveMonitor インスタンスを作成した直後) のセルに、以下のようなプロットが表示される。プロットは update が呼ばれるたびに更新される。
また、各プロットはインタラクティブにズーム/パンといった操作ができる (左上)。
まとめ
これでリアルタイムに進捗を眺めて楽しむことができる。
- 作者: 岡谷貴之
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
Chainer で Deep Learning: model zoo で R-CNN やりたい
ニューラルネットワークを使ったオブジェクト検出の手法に R-CNN (Regions with CNN) というものがある。簡単にいうと、R-CNN は以下のような処理を行う。
- 入力画像中からオブジェクトらしい領域を検出し切り出す。
- 各領域を CNN (畳み込みニューラルネットワーク) にかける。
- 2での特徴量を用いて オブジェクトかどうかをSVMで判別する。
R-CNN については 論文著者の方が Caffe (Matlab) での実装 (やその改良版) を公開している。
が、自分は Matlab のライセンスを持っていないので Python でやりたい。Python でやるなら 今 流行りの Chainer を使ってみたい。その試行の記録。
準備
とりあえず論文の再現は目指さず、R-CNN っぽい処理 (オブジェクト検出 -> CNN) を Python から回せるようにしたい。
まずは オブジェクト検出について調べてみると、上述の論文 & 実装で利用されている Selective Search という手法は ( Python の wrapper はあるが) Matlab なしでは利用できないようだ。Python から使える代替手法として GOP (Geodesic Object Proposals) という方法を見つけた。著者がパッケージも作っているのだが、どうもオブジェクト検出の方法がよくわからない、、、が矩形を切り取ることはできそうだ。
ということで今回は以下の処理が流せるようにしたい。モデルの修正やファインチューニングはせず、 caffenet をそのまま使う(済まぬ、、)。環境は EC2 の GPU インスタンス上に作成する。
- オブジェクトっぽい矩形をクロップ (Geodesic K-means)
- クロップした領域を CNN で判別 (
Chainermodel zoo で caffenet を利用) - 判別した結果を描画
都合上、 2 -> 1 -> 3 の順で記載する。
Chainer model zoo の利用
本日時点で PyPI にリリースされている v1.0.1 標準にはない機能のため、GitHub からダウンロードする。自分が利用したのは 本日時点 5222fe572b のリビジョン。
model zoo のダウンロード
chainer/examples/modelzoo から以下を実行してモデル と mean ファイルをダウンロードする。
$ python download_model.py caffenet $ python download_mean_file.py
また、データセットの ID / ラベルを含む synset_words.txt を別途ダウンロードしておく。これは以下のファイルに含まれている。
$ wget http://dl.caffe.berkeleyvision.org/caffe_ilsvrc12.tar.gz $ mkdir ilsvrc $ tar zxvf caffe_ilsvrc12.tar.gz -C ilsvrc
model zoo の読み込み
以降は IPython Notebook で。最初に必要なパッケージをロードする。
%matplotlib inline import os import sys import numpy as np import chainer from chainer import cuda import chainer.functions as F from chainer.functions import caffe import matplotlib.pyplot as plt
次に、先ほどダウンロードした caffenet のパスを指定し、Chainer のモデルとして読み込む。
また、synset_words.txt を np.array として読み込む。
base = os.path.join('chainer', 'examples', 'modelzoo') model = os.path.join(base, 'bvlc_reference_caffenet.caffemodel') func = caffe.CaffeFunction(model) cuda.init(0) func.to_gpu() synset_words = np.array([line[:-1] for line in open(os.path.join(base, 'ilsvrc', 'synset_words.txt'), 'r')]) synset_words[:10] # array(['n01440764 tench, Tinca tinca', # 'n01443537 goldfish, Carassius auratus', # 'n01484850 great white shark, white shark, man-eater, man-eating shark, Carcharodon carcharias', # 'n01491361 tiger shark, Galeocerdo cuvieri', # 'n01494475 hammerhead, hammerhead shark', # 'n01496331 electric ray, crampfish, numbfish, torpedo', # 'n01498041 stingray', 'n01514668 cock', 'n01514859 hen', # 'n01518878 ostrich, Struthio camelus'], # dtype='|S131')
続けて、画像に対する caffenet の判別結果を返す関数を準備する。画像のリストを入力とし、一度に入力された画像群に対してはバッチで処理するようにした。
in_size = 224 cropwidth = 256 - in_size start = cropwidth // 2 stop = start + in_size # caffenet meannpy = os.path.join(base, 'ilsvrc_2012_mean.npy') mean_image = np.load(meannpy) mean_image = mean_image[:, start:stop, start:stop].copy() def predict(images): """画像のリストに対する判別を行う""" global mean_image, in_size, cropwidth, start, stop def swap(x): x = np.array(x)[:, :, ::-1] x = np.swapaxes(x, 0, 2) x = np.swapaxes(x, 1, 2) return x if not isinstance(images, list): images = [images] batch_size = len(images) x_data = np.ndarray((batch_size, 3, in_size, in_size), dtype=np.float32) for i, image in enumerate(images): image = swap(image) image = image[:, start:stop, start:stop].copy().astype(np.float32) image -= mean_image x_data[i] = image x_data = cuda.to_gpu(x_data) x = chainer.Variable(x_data, volatile=True) y, = func(inputs={'data': x}, outputs=['fc8'], train=False) y.data = cuda.to_cpu(y.data) indexer = y.data.argmax(axis=1) return synset_words[indexer], y.data.max(axis=1)
動作確認
自分は ImageNet へのアクセス権を持っていないため、wikimedia で適当な画像をみつくろって試す。そのため、URL を指定して画像をダウンロードする関数、画像を caffenet へ入力するためにリサイズする関数を作成する。
def get_image(url): """URLから画像をダウンロード""" import urllib import StringIO import PIL.Image as Image return Image.open(StringIO.StringIO(urllib.urlopen(url).read())).convert("RGB") def resize(image): """画像をリサイズ""" import PIL.Image as Image return image.resize((256, 256), Image.ANTIALIAS)
いくつか結果を添付する。
url1 = "https://upload.wikimedia.org/wikipedia/commons/0/0e/BIRD_PARK_8_0189.jpg"
img = get_image(url1)
plt.imshow(img)
img = resize(img) plt.imshow(img)
predict(img) # (array(['n01806143 peacock'], # dtype='|S131'), array([ 26.1235218], dtype=float32))
別の画像。
url2 = "https://upload.wikimedia.org/wikipedia/commons/thumb/f/f4/Ostrich_Ngorongoro_05.jpg/640px-Ostrich_Ngorongoro_05.jpg"
img = get_image(url2)
plt.imshow(img)
img = resize(img) plt.imshow(img)
predict(img) # (array(['n01518878 ostrich, Struthio camelus'], # dtype='|S131'), array([ 22.44941521], dtype=float32))
定量的に測った訳ではないが、魚類は正解しにくい気がする。
url3 = "https://upload.wikimedia.org/wikipedia/commons/thumb/f/ff/Peixe010eue.jpg/640px-Peixe010eue.jpg"
img = get_image(url3)
plt.imshow(img)
img = resize(img) plt.imshow(img)
predict(img) # (array(['n01443537 goldfish, Carassius auratus'], # dtype='|S131'), array([ 24.29219818], dtype=float32))
複数画像を処理する場合は、画像のリストを渡す。
images = [resize(get_image(url1)), resize(get_image(url2)), resize(get_image(url3))] labels, scores = predict(images) labels # array(['n01806143 peacock', 'n01518878 ostrich, Struthio camelus', # 'n01443537 goldfish, Carassius auratus'], # dtype='|S131') scores # array([ 26.1235218 , 22.44941521, 24.29219818], dtype=float32)
矩形のクロップ (Geodesic K-means)
Geodesic Object Proposals 著者のサイトから Code, Data をダウンロードする。README が少しわかりにくいが、以下の方法でインストールできた。
# eigen のダウンロード $ wget http://bitbucket.org/eigen/eigen/get/3.2.5.zip $ unzip 3.2.5.zip # GOP のダウンロード $ wget http://googledrive.com/host/0B6qziMs8hVGieFg0UzE0WmZaOW8/code/gop_1.3.zip $ unzip gop_1.3.zip # GOP の解凍ディレクトリ/external/eigen に Eigen ディレクトリをコピー $ mkdir gop_1.3/external/eigen $ cp -r eigen-eigen-bdd17ee3b1b3/Eigen/ gop_1.3/external/eigen/ $ ls gop_1.3/external/eigen # Eigen $ cd gop_1.3/build # -DUSE_PYTHON で Python 2.x とのリンクを指定 $ cmake .. -DCMAKE_BUILD_TYPE=Release -DUSE_PYTHON=2 $ sudo make install
また、ダウンロードした Data は gop_1.3/data ディレクトリの中に入れておく。
$ unzip gop_data.zip
GOP を IPython Notebook から使ってみる。GOP は site-packages にインストールされないため、ロードするには sys に対して path を追加してやる必要がある。矩形を選択する方法のうち、今のところ動かせたのは Geodesic K-means のみなのでそれを使う。GOP では Geodesic K-means の後にオブジェクト検出を行うようなのだが、その結果を利用して矩形選択する方法はまだわかってない (そのため、オブジェクト検出はまだできていないという理解だが Proposal の処理を追っていないので自信はない)。
sys.path.append('./gop_1.3/src/') from gop import * from util import * prop = proposals.Proposal(setupLearned(200, 4, 0.8)) detector = contour.MultiScaleStructuredForest() detector.load("./gop_1.3/data/sf.dat") def get_boundaries(image): # 画像データを直接入力とする方法が不明のため、一時ファイルに保存 img.save('tmp.png') s = segmentation.geodesicKMeans(imgproc.imread('tmp.png'), detector, 100) b = prop.propose(s) boxes = s.maskToBox(b) return boxes
以下の画像を使って確認してみる。
url = "https://upload.wikimedia.org/wikipedia/commons/8/8f/Dogs_playing_on_the_beach_in_the_sand.jpg"
img = get_image(url)
plt.imshow(img)
fig, ax = plt.subplots(1, 1) ax.imshow(img) boxes = get_boundaries(img) def plot_boxes(ax, boxes, labels=None): """ax への矩形とラベルの追加""" if labels is None: labels = [None] * len(boxes) history = [] from matplotlib.patches import FancyBboxPatch for box, label in zip(boxes, labels): coords = (box[0], box[1]) b = FancyBboxPatch(coords, box[2]-box[0], box[3]-box[1], boxstyle="square,pad=0.", ec="b", fc="none", lw=0.5) mindist = 100000 if len(history) > 0: mindist = min([sum((box - h) ** 2) for h in history]) # ほぼ重なる矩形は描画しない if mindist > 30000: if label is not None: ax.text(coords[0], coords[1], label, color='b') ax.add_patch(b) history.append(box) plot_boxes(ax, boxes)
クロップされた画像はこんな感じ。
cropped = img.crop(boxes[13])
plt.imshow(cropped)
リサイズして caffenet にかける。犬という意味では近いが、品種は当たってない。
predict(resize(cropped)) # (array(['n02105412 kelpie'], # dtype='|S131'), array([ 10.38668251], dtype=float32))
各領域への CNN の適用
Geodesic K-means で抽出された各領域をクロップした画像のリストを作り、caffenet に渡す。
images = [resize(img.crop(box)) for box in boxes] labels = [] scores = [] nunit = len(images) / 5 unit = len(images) / nunit for i in range(nunit+1): l, s = predict(images[i*unit:min(i*unit+unit, len(images))]) labels.extend(l.tolist()) scores.extend(s.tolist()) import pandas as pd df = pd.DataFrame({'Labels': labels, 'Scores': scores}) df = df.sort('Scores', ascending=False) df.head()
判別結果のうちスコアがよい部分をラベル付きで描画する。やはりうまくオブジェクトが切り出せていない感じがする。
fig, ax = plt.subplots(1, 1) ax.imshow(img) heads = df.head(n=50) labels = [l.split(' ', 1)[1] for l in heads['Labels'].tolist()] plot_boxes(ax, boxes[heads.index], labels=labels)
まとめ
今回、Python から以下の処理をまわすことはできた。
- オブジェクトっぽい矩形をクロップ
- クロップした領域を CNN で判別
- 判別した結果を描画
が、判別結果自体は微妙な感じなので、以下2つについてはもう少し調べたい。
同日追記
パラメータ調整したら矩形抽出は多少ましになったような、、(画像末端まで選択されてしまうことが減っている)。元論文では候補を 2000 個抽出しているため、同程度以上にしておけばよさそう。
- 作者: 岡谷貴之
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
Python pandas データのイテレーションと関数適用、pipe
pandas ではデータを 列 や 表形式のデータ構造として扱うが、これらのデータから順番に値を取得 (イテレーション) して何か操作をしたい / また 何らかの関数を適用したい、ということがよくある。このエントリでは以下の 3 つについて整理したい。
- イテレーション
- 関数適用
pipe(0.16.2 で追加)
それぞれ、Series、DataFrame、GroupBy (DataFrame.groupbyしたデータ) で可能な操作が異なるため、順に記載する。
まずは必要なパッケージを import する。
import numpy as np import pandas as pd
イテレーション
Series
Series は以下 2つのイテレーション用メソッドを持つ。各メソッドの挙動は以下のようになる。
図で表すとこんな感じ。矢印が処理の方向、枠内が 1 処理単位。
s = pd.Series([1, 2, 3], index=['a', 'b', 'c']) for v in s: print(v) # 1 # 2 # 3 for i, v in s.iteritems(): print(i) print(v) print('') # a # 1 # # b # 2 # # c # 3
DataFrame
DataFrame は以下 3つのイテレーション用メソッドを持つ。同様に挙動を示す。
__iter__:DataFrameの列名 (columns) のみをイテレーションiteritems:DataFrameの列名と 列の値 (Series) からなるtupleをイテレーションiterrows:DataFrameの行名と 行の値 (Series) からなるtupleをイテレーション
df = pd.DataFrame({'A': [1, 2, 3], 'B': [4, 5, 6]}, index=['a', 'b', 'c'])
df
# A B
# a 1 4
# b 2 5
# c 3 6
for col in df:
print(col)
# A
# B
for key, column in df.iteritems():
print(key)
print(column)
print('')
# A
# a 1
# b 2
# c 3
# Name: A, dtype: int64
#
# B
# a 4
# b 5
# c 6
# Name: B, dtype: int64
for key, row in df.iterrows():
print(key)
print(row)
print('')
# a
# A 1
# B 4
# Name: a, dtype: int64
#
# b
# A 2
# B 5
# Name: b, dtype: int64
#
# c
# A 3
# B 6
# Name: c, dtype: int64
GroupBy
__iter__:GroupByのグループ名と グループ (DataFrameもしくはSeries) からなるtupleをイテレーション
df = pd.DataFrame({'group': ['g1', 'g2', 'g1', 'g2'],
'A': [1, 2, 3, 4], 'B': [5, 6, 7, 8]},
columns=['group', 'A', 'B'])
df
# group A B
# 0 g1 1 5
# 1 g2 2 6
# 2 g1 3 7
# 3 g2 4 8
grouped = df.groupby('group')
for name, group in grouped:
print(name)
print(group)
print('')
# g1
# group A B
# 0 g1 1 5
# 2 g1 3 7
#
# g2
# group A B
# 1 g2 2 6
# 3 g2 4 8
関数適用
Series
Series の各値に対して 関数を適用する方法は以下の 2 つ。挙動はほぼ一緒だが、関数適用する場合は apply を使ったほうが意図が明確だと思う
Series.apply:Seriesの各値に対して関数を適用。Series.map:Seriesの各値を、引数を用いてマッピングする。引数として、dictやSeriesも取れる。
s = pd.Series([1, 2, 3], index=['a', 'b', 'c']) s.apply(lambda x: x * 2) # a 2 # b 4 # c 6 # dtype: int64 # apply の引数には、Series の値そのものが渡っている s.apply(type) # a <type 'numpy.int64'> # b <type 'numpy.int64'> # c <type 'numpy.int64'> # dtype: object # 関数が複数の返り値を返す場合、結果は tuple の Series になる s.apply(lambda x: (x, x * 2)) # a (1, 2) # b (2, 4) # c (3, 6) # dtype: object # 結果を DataFrame にしたい場合は、返り値を Series として返す s.apply(lambda x: pd.Series([x, x * 2], index=['col1', 'col2'])) # col1 col2 # a 1 2 # b 2 4 # c 3 6 # map の挙動は apply とほぼ同じ (map では結果を DataFrame にすることはできない) s.map(lambda x: x * 2) # a 2 # b 4 # c 6 # dtype: int64 s.map(type) # a <type 'numpy.int64'> # b <type 'numpy.int64'> # c <type 'numpy.int64'> # dtype: object # map は 関数以外に、 mapping 用の dict や Series を引数として取れる s.map(pd.Series(['x', 'y', 'z'], index=[1, 2, 3])) # a x # b y # c z # dtype: object
DataFrame
DataFrame に対して 関数を適用する方法は以下の 2 つ。
DataFrame.apply:DataFrameの各列もしくは各行に対して関数を適用。行 / 列の指定はaxisキーワードで行う。DataFrame.applymap:DataFrameの各値に対して関数を適用。
df = pd.DataFrame({'A': [1, 2, 3], 'B': [4, 5, 6]}, index=['a', 'b', 'c'])
df
# A B
# a 1 4
# b 2 5
# c 3 6
# 各列に対して関数適用
df.apply(lambda x: np.sum(x))
A 6
B 15
dtype: int64
# 各行に対して関数適用
df.apply(lambda x: np.sum(x), axis=1)
a 5
b 7
c 9
dtype: int64
# 各値に対して関数適用
df.applymap(lambda x: x * 2)
# A B
# a 2 8
# b 4 10
# c 6 12
# apply で適用される関数には、各列もしくは各行が Series として渡される
df.apply(type)
# A <class 'pandas.core.series.Series'>
# B <class 'pandas.core.series.Series'>
# dtype: object
# applymap で適用される関数には、値そのものが引数として渡される
df.applymap(type)
# A B
# a <type 'numpy.int64'> <type 'numpy.int64'>
# b <type 'numpy.int64'> <type 'numpy.int64'>
# c <type 'numpy.int64'> <type 'numpy.int64'>
GroupBy
GroupBy については、GroupBy.apply で各グループに関数を適用できる。
df = pd.DataFrame({'group': ['g1', 'g2', 'g1', 'g2'],
'A': [1, 2, 3, 4], 'B': [5, 6, 7, 8]},
columns=['group', 'A', 'B'])
df
# group A B
# 0 g1 1 5
# 1 g2 2 6
# 2 g1 3 7
# 3 g2 4 8
grouped = df.groupby('group')
grouped.apply(np.mean)
# A B
# group
# g1 2 6
# g2 3 7
補足 処理最適化のため、対象となるグループの数 == 関数適用の実行回数とはならないので注意。関数中で破壊的な処理を行うと意図しない結果になりうる。
# 適用される関数 def f(x): print('called') return x # グループ数は 2 grouped.ngroups # 2 # f の実行は 3 回 grouped.apply(f) # called # called # called
pipe
先日 リリースされた v0.16.2 にて pipe メソッドが追加された。これは R の {magrittr} というパッケージからインスパイアされたもので、データへの連続した操作を メソッドチェイン (複数のメソッドの連続した呼び出し) で記述することを可能にする。
Series.pipe、DataFrame.pipe それぞれ、x.pipe(f, *args, **kwargs) は f(x, *args, **kwargs) と同じ。つまり、データ全体に対する関数適用になる。
補足 GroupBy.pipe は v0.16.2 時点では存在しない。
# 渡される型は呼び出し元のインスタンス s.pipe(type) # pandas.core.series.Series df.pipe(type) # pandas.core.frame.DataFrame np.random.seed(1) df = pd.DataFrame(np.random.randn(10, 10)) # DataFrame を引数として heatmap を描く関数を定義 def heatmap(df): import matplotlib.pyplot as plt fig, ax = plt.subplots() return ax.pcolor(df.values, cmap='Greens') # heatmap(df) と同じ。 df.pipe(heatmap)
まとめ
イテレーション、関数適用、pipe について整理した。特に関数適用は データの前処理時に頻出するため、パターンを覚えておくと便利。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas + folium + Jupyter でリーフレット / コロプレス図を描きたい
引き続き、 R の可視化を Python に持ってくるシリーズ。R には以下のようなパッケージがあり、地図上へのリーフレット配置やコロプレス図の描画がカンタンにできる。それぞれの概要はリンク先を。
{leaflet}: リーフレット配置{choroplethr}: コロプレス図の描画
これを Python でやりたい。調べてみると folium というパッケージが上記 両方をサポートしているようなので使ってみる。
インストール
pip で。
pip install folium
準備
以降の操作は Jupyter Notebook から行う。まずはパッケージをロードする。
import numpy as np import pandas as pd import folium
folium は プロット結果を html としてエクスポートすることを想定して作成されているようだ。そのため、結果を Jupyter 上に埋め込みたい場合は 以下のような関数を定義する必要がある。
from IPython.display import HTML def inline_map(m): # 中間生成される json が必要なプロットがあるため、一度 html として書き出し m.create_map(path='tmp.html') iframe = '<iframe srcdoc=\"{srcdoc}\" style=\"width: 100%; height: 400px; border: none\"></iframe>' return HTML(iframe.format(srcdoc=m.HTML.replace('\"', '"')))
リーフレット
手順は以下のようになる。
folium.Mapで地図を描画する範囲を緯度経度/ズームにより指定Map.simple_markerで緯度経度を指定してリーフレットを配置 (複数配置する場合は繰り返し)Map.create_mapで地図を含むhtmlを生成 (Jupyter上に描画する場合は上で定義した関数inline_mapを呼ぶ )
m = folium.Map(location=[33.763, -84.392], zoom_start=17) m.simple_marker([33.763006, -84.392912], popup='World of Coca-Cola') inline_map(m)
補足 Jupyter 上ではスクロール / 拡大縮小できる。
既定では OpenStreetMap が利用されるが、Mapbox や maps.stamen を使うこともできる。また、リーフレットのマーカーとしては以下のものが利用できる。
- Simple Markers: シンプルなマーカー (上のもの)
- Circle Markers: 円形のマーカー
- Polygon Markers: 多角形のマーカー
- Lat/Lng Popups: 緯度経度を表示するマーカー
- Click-for-Marker: クリックで配置可能なマーカー
- Vincent Popups: Vincent によるプロットを埋め込めるマーカー
それぞれの描画サンプルは README で確認することができる。
コロプレス図
コロプレス図を描くには以下２つのデータソースが必要である。
GeoJSONもしくはTopoJSON形式のファイル- コロプレス図を色分けするための値を含む
pandasのDataFrame
サンプルとして 国別のマクドナルドの店舗数 をプロットする。
GeoJSON ファイルの準備
国別にプロットするため、国別の GeoJSON ファイルがほしい。以下リポジトリの countries.geo.json をローカルに保存して使う。
中身は 以下のように ISO 3166-1 alpha-3 の国コードを id としたデータとなっている。
{"type":"FeatureCollection",
"features":[{"type":"Feature","id":"AFG",
"properties":{"name":"Afghanistan"},
"geometry":{"type":"Polygon","coordinates":[[[61.210817,35.650072],[62.230651,35.270664],...
DataFrame の準備
DataFrame は 上で準備した GeoJSON と紐づけるためのキー (一般には GeoJSON の id ) と値の 2 列をもつ必要がある。 国別のマクドナルドの店舗数 データを pd.read_html で読み込む。
url = "https://en.wikipedia.org/wiki/List_of_countries_with_McDonald%27s_restaurants" df = pd.read_html(url, header=0, index_col=0)[0] # 列名を変更 df.columns = ['Country', 'Date', 'First outlet location', 'Number', 'Source', 'Note', 'CEO'] df[['Country', 'Number']].head() # Country Number # # # 1 United States 14267 # 2 Canada 1427 # 3 Puerto Rico 108 # 4 U.S. Virgin Islands 6 # 5 Costa Rica 54
こちらのデータには国名が入っているため、GeoJSON と紐づけるためには ISO 国コードに変換する必要がある。変換には pycountry を使う。インストールしていない方は pip で。
import pycountry pycountry.countries.get(name='Japan').alpha3 # 'JPN' # データ中の国名が pycountry のものと違う場合の mapper countries = {'United Kingdom': 'United Kingdom', 'Russia': 'Russian Federation', 'South Korea': 'Korea, Republic of', 'Taiwan': 'Taiwan, Province of China', 'Vietnam': 'Viet Nam', } def f(x): for k, v in pd.compat.iteritems(countries): if k in x: x = v try: return pycountry.countries.get(name=x).alpha3 except KeyError: return np.nan # 国コードへの変換 df.loc[:, 'Code'] = df['Country'].apply(f) # 国コードに変換できなかったデータは除外 df = df.dropna(subset=['Code']) # 欠損値を 0 でパディング df.loc[:, 'Number'] = df['Number'].fillna('0') # 数値に変換できない文字列があるため、余計な文字を削除 df.loc[:, 'Number'] = df['Number'].str.replace('[+,]', '') # 数値 (float) 型に変換 df.loc[:, 'Number'] = df['Number'].astype(float) df[['Code', 'Country', 'Number']].sort('Number', ascending=False).head() # Code Country Number # # # 1 USA United States 14267 # 8 JPN Japan 3164 # 49 CHN China 2000 # 11 DEU Germany 1468 # 2 CAN Canada 1427
これで、国コード / プロットする値をもつ DataFrame が準備できた。アメリカすごいな、、、。
コロプレス図の描画
Map.geo_json で コロプレス図を描画するレイヤーを Map に追加できる。ここで使っている引数の意味は以下。
geo_path:GeoJSONファイルのパスdata: 色分けのための値をもつDataFramecolumns:data中GeoJSONと紐づけるキー / 値 を含む列名key_on:GeoJSON側で紐付けに使うキーthreshold_scale: 色分けをする際の閾値fill_color色分けに使う color-brewer の名前reset: 既存のレイヤがある場合に削除する
m = folium.Map(location=[10, 35], zoom_start=1.5) geojson = r'countries.geo.json' m.geo_json(geo_path=geojson, data=df, columns=['Code', 'Number'], key_on='feature.id', threshold_scale=[1, 100, 500, 1000, 2000, 4000], fill_color='BuPu', reset=True) inline_map(m)
まとめ
folium を使えば リーフレット / コロプレス図の描画がカンタンにできる。Jupyter 上で Javascript を使用してのデータ可視化、結構使えるのでは。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas のデータを Highcharts/Highstock + Jupyter でプロットしたい
R を使っている方はご存知だと思うが、R には {htmlwidgets} というパッケージがあり、R 上のデータを任意の Javascript ライブラリを使ってプロットすることが比較的カンタンにできる。{htmlwidgets} って何？という方には こちらの説明がわかりやすい。
同じことを Python + pandas を使ってやりたい。サンプルとして利用する Javascript ライブラリは 上の資料と同じく Highcharts、Highstock にする。
補足 pandas-highcharts という Python パッケージもあるが、このエントリでは任意の Javascript ライブラリで使えるであろう方法を記載する。
Highcharts でのプロット
以降の操作は Jupyter Notebook 上で行う。まずは必要パッケージをロードする。
import numpy as np import pandas as pd from IPython.display import HTML
続けて、Highcharts、Highstock を読み込む。これは %%html magic を使うのが楽。
%%html
<script src="http://code.highcharts.com/highcharts.js"></script>
<script src="http://code.highcharts.com/stock/highstock.js"></script>
<script src="http://code.highcharts.com/modules/exporting.js"></script>
Highcharts でプロットする準備ができたため、%%html magic を使って サンプル Your first chart に記載のサンプルをプロットしてみる。これでプロット時のスクリプト / データ構造が確認できる。
%%html
<div id="container" style="width:100%; height:400px;"></div>
<script>
plot = function () {
$('#container').highcharts({
chart: {
type: 'bar'
},
title: {
text: 'Fruit Consumption'
},
xAxis: {
categories: ['Apples', 'Bananas', 'Oranges']
},
yAxis: {
title: {
text: 'Fruit eaten'
}
},
series: [{
name: 'Jane',
data: [1, 0, 4]
}, {
name: 'John',
data: [5, 7, 3]
}]
});
};
plot();
</script>
補足 Jupyter から実行すればアニメーションする。
pandas のデータをプロットしたい場合は、上のスクリプトと同じように .highcharts() の引数にあわせた形式でデータを渡し、HTML としてレンダリングしてやればうまくいきそうだ。pandas で元データとなる DataFrame を定義して、上のスクリプトの形式に変換していく。
df = pd.DataFrame({'Jane': [1, 0, 4], 'John': [5, 7, 3]},
index=['Apples', 'Bananas', 'Oranges'])
df
# Jane John
# Apples 1 5
# Bananas 0 7
# Oranges 4 3
スクリプトは文字列結合で作ってもよいが、引数となる json 形式に対応する辞書型のデータをつくってからjson.dumps したほうが簡単だろう。DataFrame.to_json でデータを直接 変換できると楽なのだが、フォーマットが違うため無理そうだ。
# NG! df.to_json() # '{"Jane":{"Apples":1,"Bananas":0,"Oranges":4}, # "John":{"Apples":5,"Bananas":7,"Oranges":3}}'
そのため、個々の要素ごとに変換を考えていく。うまいこと辞書ができたら json.dumps する。
[{'name': c, 'data': col.tolist()} for c, col in df.iteritems()]
# [{'data': [1, 0, 4], 'name': 'Jane'}, {'data': [5, 7, 3], 'name': 'John'}]
chartdict = {'chart': {'type': 'bar'},
'title': {'text': 'Fruit Consumption'},
'xAxis': {'categories': df.index.tolist()},
'yAxis': {'title': {'text': 'Fruit eaten'}},
'series': [{'name': c, 'data': col.tolist()} for c, col in df.iteritems()]
}
import json
json.dumps(chartdict)
# '{"series": [{"data": [1, 0, 4], "name": "Jane"}, {"data": [5, 7, 3], "name": "John"}],
# "yAxis": {"title": {"text": "Fruit eaten"}}, "chart": {"type": "bar"},
# "xAxis": {"categories": ["Apples", "Bananas", "Oranges"]},
# "title": {"text": "Fruit Consumption"}}'
あとは 必要な HTML / Javascript のテンプレートを文字列として作って format すればよい。
template = """ <script src="http://code.highcharts.com/highcharts.js"></script> <script src="http://code.highcharts.com/modules/exporting.js"></script> <div id="{chart}" style="width:100%; height:400px;"></div> <script type="text/javascript"> plot = function () {{ $("#{chart}").highcharts({data}); }}; plot(); </script> """ HTML(template.format(chart='container2', data=json.dumps(chartdict))) # 略
Highstock でのプロット
同様に、Highstock へプロットすることもできる。サンプル Single line series を pandas のデータを使ってプロットしてみる。
補足 株価の取得には以下のパッケージを使う。
import japandas as jpd toyota = jpd.DataReader(7203, 'yahoojp', start='2015-01-01') toyota.head() # 始値 高値 安値 終値 出来高 調整後終値* # 日付 # 2015-01-05 7565 7575 7416 7507 9515300 7507 # 2015-01-06 7322 7391 7300 7300 12387900 7300 # 2015-01-07 7256 7485 7255 7407 11465400 7407 # 2015-01-08 7500 7556 7495 7554 10054500 7554 # 2015-01-09 7630 7666 7561 7609 10425400 7609
Highstock に渡すデータの形式は以下のサイトがわかりやすい。引数としては [UNIX時間(ミリ秒), 始値, 高値, 安値, 終値] を入れ子のリストにして渡せばよいようだ。
時刻は UNIX時間で渡す必要があるので、少し操作が必要だ。上で取得した DataFrame は日時型の Index である DatetimeIndex を持っている。DatetimeIndex はUNIXエポックを基準とした現在時刻をナノ秒で保存しているため、int 型に変換して 1000000 で割ればUNIX時間(ミリ秒)となる。したがって、Highstock へ渡すデータは以下のようにして作れる。
toyota.index.astype(int) / 1000000 # array([1420416000000, 1420502400000, 1420588800000, 1420675200000, # 1420761600000, 1421107200000, 1421193600000, 1421280000000, # .... # 1433721600000, 1433808000000, 1433894400000, 1433980800000, # 1434067200000]) toyota['time'] = toyota.index.astype(int) / 1000000 toyota[['time', u'始値', u'高値', u'安値', u'終値']].values.tolist() # [[1420416000000, 7565, 7575, 7416, 7507], # [1420502400000, 7322, 7391, 7300, 7300], # ... # [1433980800000, 8250, 8326, 8241, 8322], # [1434067200000, 8387, 8394, 8329, 8394]]
描画する。アニメーションも含め、うまく動いているようだ。
chartdict = {'rangeSelector': {'selected': 1},
'title': {'text' : 'Stock Price'},
'series': [{'name' : u'トヨタ',
'data': toyota[['time', u'始値', u'高値', u'安値', u'終値']].values.tolist(),
'tooltip': {'valueDecimals': 2}}]}
template = """
<div id="{chart}" style="width:100%; height:400px;"></div>
<script type="text/javascript">
plot = function () {{
$('#{chart}').highcharts('StockChart', {data});
}};
plot();
</script>
HTML(template.format(chart='container3', data=json.dumps(chartdict)))
まとめ
pandas のデータを任意の Javascript ライブラリでプロットする際には、
- 当該の Javascript ライブラリが利用するデータ構造を確認する
- そのデータ構造にあうように
pandasのデータを変換し、辞書型を作る json.dumpsでスクリプトに渡す
とやればだいたいできると思う。
補足 ここで今 話題の spyre 上でもプロットしたいな？と思ったのだが、spyre で普通に HTML としてレンダリングするとうまく動かない。spyre が内部で使っている cherrypy も DOM を書き換えているのが原因かと思っているのだが、自分にはそのあたりの知識がまったくないのでよくわからない。できた方がいれば教えてください。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python spyre によるデータ分析結果のWebアプリ化
R を使っている方はご存知だと思うが、R には {Shiny} というパッケージがあり、データ分析の結果を インタラクティブな Web アプリとして共有することができる。{Shiny} って何？という方には こちらの説明がわかりやすい。
Python でも {Shiny} のようなお手軽可視化フレームワークがあるといいよね、とたびたび言われていたのだが、spyre という なんかそれっぽいパッケージがあったので触ってみたい。
インストール
pip で。
pip install dataspyre
使い方
現時点で ドキュメンテーションはない ので、README と examples ディレクトリを見る。サンプルとして株価を取得してプロットするWebアプリを作ってみたい。spyre で Webアプリを作る手順は以下の3つ。
spyre.server.Appを継承したクラスを作る。- 描画をコントロール/指示するクラス変数を指定する。
- 描画を行うメソッド
getData,getPlotを書く。メソッド名は描画内容 (output_type) ごとに固定。
最初は examples のものを書き換えながら作るのが楽。各クラス変数/メソッドは相互に関連するため、プログラム全体を示した上で必要と思われる箇所にコメントを入れた。
補足 株価の取得には以下のパッケージを使う。
#!/usr/bin/env python # -*- coding: utf-8 -*- from spyre import server import pandas as pd pd.options.display.mpl_style = 'default' # あらかじめデータを取得しておく # 終値のみを取得し、一つのDataFrameに結合 import japandas as jpd toyota = jpd.DataReader(7203, 'yahoojp', start='2015-01-01')[[u'終値']] toyota.columns = [u'トヨタ'] honda = jpd.DataReader(7267, 'yahoojp', start='2015-01-01')[[u'終値']] honda.columns = [u'ホンダ'] df = toyota.join(honda) # spyre.server.App を継承したクラスを作る class StockExample(server.App): title = u"株価のプロット" # 左側のペインに表示する UI 要素を辞書のリストで指定 # ここではドロップダウン一つだけを表示 inputs = [{"input_type":'dropdown', # ドロップダウン自体の表示ラベル "label": 'Frequency', # ドロップダウンの選択項目を指定 # label はドロップダウン項目の表示ラベル # value は各項目が選択された時にプログラム中で利用される値 "options" : [ {"label": "月次", "value":"M"}, {"label": "週次", "value":"W"}, {"label": "日次", "value":"B"}], # 各 UI 要素の入力は各描画メソッド (getData, getPlot) に # 辞書型の引数 params として渡される # その辞書から値を取り出す際のキー名 "variable_name": 'freq', "action_id": "update_data" }] # 画面を更新する設定 controls = [{"control_type" : "hidden", "label" : "update", "control_id" : "update_data"}] # 描画するタブの表示ラベルを文字列のリストで指定 tabs = [u"トヨタ", u"ホンダ", u"データ"] # tabs で指定したそれぞれのタブに描画する内容を辞書のリストで指定 outputs = [{"output_type" : "plot", # matplotlib のプロットを描画する "output_id" : "toyota", # 描画するタブに固有の id "control_id" : "update_data", "tab" : u"トヨタ", # 描画するタブの表示ラベル (tabs に含まれるもの) "on_page_load" : True }, {"output_type" : "plot", "output_id" : "honda", "control_id" : "update_data", "tab" : u"ホンダ", "on_page_load" : True }, {"output_type" : "table", # DataFrameを描画する "output_id" : "table_id", "control_id" : "update_data", "tab" : u"データ", "on_page_load" : True }] def getData(self, params): """ output_type="table" を指定したタブを描画する際に呼ばれるメソッド DataFrameを返すこと params は UI 要素の入力 + いくつかのメタデータを含む辞書 UI 要素の入力は inputs で指定した variable_name をキーとして行う """ # ドロップダウンの値を取得 # 値にはユーザの選択によって、options -> value で指定された M, W, B いずれかが入る freq = params['freq'] # freq でグループ化し平均をとる tmp = df.groupby(pd.TimeGrouper(freq)).mean() return tmp def getPlot(self, params): """ output_type="plot" を指定したタブを描画する際に呼ばれるメソッド matplotlib.Figureを返すこと """ tmp = self.getData(params) # 同じ output_type で複数のタブを描画したい場合は、 params に含まれる # output_id で分岐させる # output_id は タブの表示ラベルではなく、outputs 中で指定した output_id if params['output_id'] == 'toyota': ax = tmp[[u'トヨタ']].plot(legend=False) return ax.get_figure() elif params['output_id'] == 'honda': ax = tmp[[u'ホンダ']].plot(legend=False) return ax.get_figure() else: raise ValueError app = StockExample() # port 9093 で Webサーバ + アプリを起動 app.launch(port=9093)
実行後、ブラウザで http://127.0.0.1:9093 を開くと以下のような画面が表示される。inputs で指定した項目が左側のメニューに、 tabs で指定した項目がタブとして表示されている。
ドロップダウンやタブの選択を切り替えると、動的にグラフやデータが更新されることがわかる。
DataFrame の描画はちょっとおかしく、index である日時自体が描画されずにその名前 ("日付") だけが表示されている。これは後日 issue あげて直そうかと思う。
06/13追記 index を描画しないのは spyre の仕様で、そのとき名前が残ってしまうのは pandas のバグです。v0.17.0で修正予定。
まとめ
お手軽可視化フレームワーク spyre を試してみた。現行の {Shiny} と比べると 機能/UI とも十分ではないが、いくつかのデータ / プロットをインタラクティブに共有する用途であれば十分使えそうだ。
使いごこちの印象として、はじめて {Shiny} を触ったあの日の気持ちを思い出したような気がする、、。
R で 状態空間モデル: {dlm} の最尤推定を可視化する
{dlm} において、状態空間モデルが最尤推定される過程がみたい。以下内容の補足的なエントリ。
「状態空間時系列分析入門」では 引き続き 第8章 に相当。
- 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
- 出版社/メーカー: シーエーピー出版
- 発売日: 2008/09
- メディア: 単行本
- 購入: 2人 クリック: 4回
- この商品を含むブログを見る
データの準備
データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。
library(dlm) # install_github('sinhrks/ggfortify') library(ggfortify) ukdrivers <- read.table('UKdriversKSI.txt', skip = 1) ukdrivers <- ts(ukdrivers[[1]], start = c(1969, 1), frequency = 12) ukdrivers <- log(ukdrivers) y <- ukdrivers
モデルの作成と状態推定値のプロット
サンプルとして、第3章3節 確率的レベルと確定的傾きをもつローカル線形トレンドモデル を作成する。
parm <- log(c(var(y), 0.001)) buildFun <- function(parm){ dlmModPoly(order = 2, dV = exp(parm[1]), dW = c(exp(parm[2]), 0)) } dlmFit <- dlmMLE(y, parm = parm, build = buildFun) # カルマンフィルタ dlmFiltered <- dlmFilter(y, buildFun(dlmFit$par)) # 平滑化 dlmSmoothed <- dlmSmooth(y, buildFun(dlmFit$par))
第8章ではモデルの推定値として 3 つの値が出てくるため、その差異がわかるようにプロットする。用語は状態空間時系列分析入門のものを利用。
p <- autoplot(y) p <- autoplot(dlmFiltered$a[, 1], colour = 'red', linetype = 'dashed', p = p) p <- autoplot(dlmFiltered$m[, 1], colour = 'red', p = p) p <- autoplot(dlmSmoothed, colour = 'blue', p = p) p + ylim(c(6.8, 8))
- 黒実線 (
y): 観測時系列 - 赤破線 (
dlmFiltered$a): カルマンフィルタの1期先予測 (予測ステップでの予測値) - 赤実線 (
dlmFiltered$m): カルマンフィルタのフィルタ化状態 (更新ステップでの予測値) - 青実線 (
dlmSmoothed): 平滑化状態
一部期間を拡大するとそれぞれの差異がわかる。これらのうち、最尤推定時の尤度計算式中で直接 利用されるのは "1期先予測" 時の予測誤差 (赤破線) とその標準偏差。
p + xlim(as.Date('1975', '%Y'), as.Date('1977', '%Y')) + ylim(c(6.8, 8))
補足 上のモデルでは FF = [1, 0] のため、観測値の推定値 == 状態推定値。
dlmFiltered$mod$FF # [,1] [,2] # [1,] 1 0 all((dlmFiltered$a %*% t(dlmFiltered$mod$FF)) == dlmFiltered$a[, 1]) # [1] TRUE
補足 カルマンフィルタの動きについては以前こんなエントリを書いた。2つ目のエントリにカルマンフィルタの予測 / フィルタ化状態がどういった感じで動いているか記載がある (平滑化は記載ない)。
予測誤差と予測誤差分散
最尤推定の尤度計算式は先のエントリに記載。計算に必要な (1期先) 予測誤差と予測誤差標準偏差は dlm:::residuals.dlmFiltered で求めることができる。
# 既定では標準化されるため、 type = 'raw' を指定して実値をとる resid <- residuals(dlmFiltered, type = 'raw') # 予測誤差 resid$res[1:12] # [1] 7.430707083 -3.827523556 0.111504168 -0.021616238 0.225665269 -0.044983932 # [7] 0.045295381 0.061233329 -0.021088564 0.049603994 0.270853650 0.007138991 # 予測誤差標準偏差 resid$sd[1:12] # [1] 4472.1359566 2236.0679824 0.1922632 0.1596384 0.1484746 0.1429560 # [7] 0.1396797 0.1375109 0.1359691 0.1348166 0.1339224 0.1332084
ふつうに計算して求めるなら、dlm::dlmFiltered インスタンスの各プロパティを利用して、
# 予測誤差 (y - (matrix(as.numeric(dlmFiltered$a), 192, 2) %*% t(dlmFiltered$mod$FF)))[1:12] # [1] 7.430707083 -3.827523556 0.111504168 -0.021616238 0.225665269 -0.044983932 # [7] 0.045295381 0.061233329 -0.021088564 0.049603994 0.270853650 0.007138991 # もしくは (y - dlmFiltered$f)[1:12] # [1] 7.430707083 -3.827523556 0.111504168 -0.021616238 0.225665269 -0.044983932 # [7] 0.045295381 0.061233329 -0.021088564 0.049603994 0.270853650 0.007138991 # 予測誤差標準偏差 f <- function(i) { FF <- dlmFiltered$mod$FF V <- dlmFiltered$mod$V diag(crossprod(dlmFiltered$D.R[i, ] * t(FF %*% dlmFiltered$U.R[[i]])) + V) } drop(t(sqrt(sapply(seq(along = dlmFiltered$U.R), f))))[1:12] # [1] 4472.1359566 2236.0679824 0.1922632 0.1596384 0.1484746 0.1429560 # [7] 0.1396797 0.1375109 0.1359691 0.1348166 0.1339224 0.1332084
最尤推定の可視化
最尤推定時のパラメータの変化を記録するため、dlm::dlmMLE を書き換えた関数を定義する。LL_ckbookcompat はテキスト記載の対数尤度 (に近いもの) を計算するための関数 ( 詳細はこちら )。
# optim 1 ステップごとに作成されたモデルを保存 models <<- list() dlmMLE_wrapper <- function (y, parm, build, method = "L-BFGS-B", ..., debug = FALSE) { logLik <- function(parm, ...) { mod <- build(parm, ...) # グローバル変数に追加 models[[length(models) + 1]] <<- mod return(dlmLL(y = y, mod = mod, debug = debug)) } return(optim(parm, logLik, method = method, ...)) } # 推定するパラメータ数が 2 のため、総和をとる範囲を変更 LL_ckbookcompat <- function(y, mod) { n <- length(y) filtered <- dlmFilter(y, mod) resid <- residuals_dlmcompat(filtered) F <- (as.numeric(resid$sd)^2)[3:n] v <- as.numeric(resid$res)[3:n] return(-0.5 * (n * log(2 * pi) + sum(log(F) + (v^2/F)))) }
続けて、アニメーション用の処理を定義して実行。
library(animation) models <<- list() dlmMLE_wrapper(y, parm = parm, build = buildFun) saveGIF({ for (mod in models[1:65]) { filtered <- dlmFilter(y, mod) p <- autoplot(filtered, fitted.colour = 'red', fitted.linetype = 'dashed') p <- autoplot(dlmSmooth(y, mod), colour = 'blue', p = p) label <- paste('CKbook LL:', round(LL_ckbookcompat(y, mod), 2), 'dlmLL:', round(dlmLL(y, mod), 2)) p <- p + annotate('text', label = label, x = as.Date('1980', '%Y'), y = 7.8, size = 6) print(p) Sys.sleep(0.01) } }, interval = 0.3, movie.name = "dlm.gif", ani.width = 500, ani.height = 300)
最尤推定時に各状態推定値がどのように変化していくかをみる。左側 "CKbookLL" はテキスト記載の対数尤度計算式 LL_ckbookcompat の値。右側 "dlmLL" は dlm::dlmLL の値 ( optim で最適化される値 )。
テキストに記載された数値より、収束時には左側の値が 0.6247935 x 観測時系列の長さ (192) = 119.9604 程度になるはずだ。
補足 {dlm} は既定では optim の際に method = 'L-BFGS-B' を利用する。可視化目的であれば Nelder-Mead 法なんかを使ったほうが 素人目には "最適化らしい" 動きになる気がする。
models <<- list() dlmMLE_wrapper(y, parm = parm, build = buildFun, method = 'Nelder-Mead') # 以降略
R で 状態空間モデル: 状態空間時系列分析入門を {rstan} で再現したい
前の記事でもリンクさせていただいているが、サイト 「状態空間時系列分析入門」をRで再現する では以下のテキストを {dlm}, {KFAS} で再現されており非常にありがたい。これらのパッケージの使い方については リンク先を読めば困らない感じだ。
自分も勉強のために似たことやりたい、、でも同じことやるのもなあ、、と考えた結果 同テキストの内容 {rstan} を使ってやってみた。
補足 Stan には状態空間表現用の関数 gaussian_dlm_obs ( 利用例 ) があるのだが、自分は使ったことがない。7章までのモデルは全て漸化式で表現されているため、それらを Stan のモデルとして記述した。
- 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
- 出版社/メーカー: シーエーピー出版
- 発売日: 2008/09
- メディア: 単行本
- 購入: 2人 クリック: 4回
- この商品を含むブログを見る
github リポジトリ
各スクリプト / モデルは GitHub にあげた。現在 7 章までのモデルと kick する R スクリプト, markdown ファイルを置いている。
モデルの一覧
各章のモデルについて、その実行結果を RPubs にあげている。各モデルの概要とリンクを以下に記載する。自分のモデリング能力不足のため いくつか課題があると考えており、特に怪しいものには * をつけた。
- 第1章 はじめに
- 第2章 ローカル・レベル・モデル
- 第3章 ローカル線形トレンド・モデル
- 第4章 季節要素のあるローカル・レベル・モデル
- 第5章 説明変数のあるローカル・レベル・モデル
- 第6章 干渉変数のあるローカル・レベル・モデル
- 第7章 英国シートベルト法とインフレーション・モデル
かんたんな説明
必要パッケージ
まず、実行には以下のCRAN非公開パッケージが必要である。後者はどうでもよいが {pforeach} はいろいろ捗るのでインストールしていない方はただちに入れるべきだ。
hoxo-m/pforeach:Stan並列化に利用。sinhrks/ggfortify: 結果のプロットに利用。
データの読み込み
データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。各データは適宜 ts インスタンス化しておく。
Stan の実行
Stan の呼び出しは全て以下のように並列化して行う。
model_file <- 'モデルを記述したファイルのパス' # モデルのコンパイルを実行 stan_fit <- stan(file = model_file, chains = 0) # Stan を並列実行する。実行結果に sflist2stanfit を適用し、stanfit インスタンスを得る fit <- pforeach(i = 1:4, .final = sflist2stanfit)({ stan(fit = stan_fit, data = standata, chains = 1, seed = i) })
また、モデルが収束したかどうかは 全ての Rhat が 1.1 未満となっているかどうかで確認する。そのままでは収束しない / しにくいモデルには、初期値 / 事前分布 / 制約を入れて収束させた。
is.converged <- function(stanfit) { summarized <- summary(stanfit) all(summarized$summary[, 'Rhat'] < 1.1) } stopifnot(is.converged(fit))
また、テキストに最尤推定の結果が記載されている場合は、以下のようにして近い値になっているかを確認した。値がずれる場合は stopifnot を外し、差異を出力する。
# is.almost.fitted の定義は各リンク先を参照 stopifnot(is.almost.fitted(mu[[1]], 7.4157))
結果のプロット
テキストのグラフを再現する形で行う。使っているパッケージの説明はこちら。
- RPubs - Plotting Time Series with ggplot2 and ggfortify
- RPubs - Plotting State Space Time Series with ggplot2 and ggfortify
課題
以下の2つについてはもう少し Stan に習熟できたら直したい。
複数の確率的状態をもつモデルをうまく収束させたい
このようなモデルは 収束しにくく、また 最尤推定した状態推定値 ( テキスト、{dlm}、{KFAS}の結果 ) とずれやすい。特に 現状 制約を入れて収束させているモデルについて、事前分布 / 初期値 / 尤度関数の範囲でなんとかしたい。
モデルの一覧で * をつけたものがこれにあたる。
テキストの尤度関数を再現したい
テキストに記載の式 (概要は以下記事参照) に揃えようとすると、自分の記述では収束しにくくなってしまう。
