メモ: Importance sampling(重点サンプリング)

概要

・サンプリングの難しい分布からのサンプリングを工夫することで期待値計算を効率的にする

Importance sampling

自分もよく分かっていない部分もあるので怪しいですが。ある期待値計算
$$
E \left( s \right) = \int s \cdot f \left( s \right) {\rm d}s
$$
をしたいけれど、\( f \left( \cdot \right) \)からサンプリングするのが難しいとする。
そんなときにサンプリングが簡単にできる\( g \left( \cdot \right) \)をもってくると
$$
E \left( s \right) = \int s \cdot f \left( s \right) {\rm d}s
= \int s \cdot \frac{f \left( s \right)}{g \left( s \right)} \cdot g \left( s \right) {\rm d}s
$$
と変形することができ、\( g \left( \cdot \right) \)からのサンプルで期待値を計算することができる。ここで\( \frac{f \left( s \right)}{g \left( s \right)} \)は重みと呼ばれる。

例題

参考サイト1に記載されている問題を考える。
「You want to simulate the mean of a standard normal distribution, truncated to the unit interval [0,1].」

[0, 1]だとImportance samplingの有り難みが分からないので、ここでは[2, 2.1]とする。
つまり、[2, 2.1]の範囲の標準正規分布の期待値を計算する。[2, 2.1]だとそもそもの範囲が狭いのと、標準正規分布で2以上の値を出すにはそれなりの試行回数が必要になると考えられるので、それなりにImportance samplingの有り難みが分かるはず

また、\( f \left( \cdot \right) \)の形は
$$
f \left( s \right) = \frac{\phi \left( x \right)}{\int_2^{2.1} \phi \left( x \right) {\rm d}z}
$$
となりまして、気合で標準正規乱数からサンプリングすることもできますが、\( f \left( \cdot \right) \)から直接サンプリングしたい場合には、はてどうやるんですかねという感じですね。\( \phi \left( \cdot \right) \)は標準正規分布の密度関数。

また、\( g \left( \cdot \right) \)は一様分布を用います。言わずもがな乱数を発生させるのが簡単なのと、[2, 2.1]の場合、\( g \left( x \right) = 10 \)となり、重みの計算も簡単になります。

Rコード

Brute force simulationと記載されているのはそのまま力技で期待値を計算する方法で、標準正規分布の乱数が[2, 2.1]の値をとった時のみ、期待値を計算するようにします。
当然ながらほとんどのサンプルがrejectされてしまうので、期待値が収束するまで時間がかかるはず。

一方でImportance samplingは一様分布を用いており、全てのサンプルを用いて期待値を計算することができる。
期待値の収束の過程をプロットしてみると以下。
黒線がBrute forceで、赤線がImportance samplingとなります。
黒線はなかなかサンプルが発生しないので収束まで時間がかかっていることが分かります。
150107_1

追記(21-Jan-2015)

効率の悪い書き方をしているような気がするので書き直し。結果を出力するために余分なコードが入っていますが基本的には同じ(はす)

参考

Lecture notes: Simulation
http://people.hss.caltech.edu/~mshum/gradio/simulation.pdf

パーティクルフィルタ
http://watanabe-www.math.dis.titech.ac.jp/~miki_t/pdf/s1114.pdf?bcsi_scan_4c5c01dba4894524=0&bcsi_scan_filename=s1114.pdf

為替高頻度データから任意のロウソク足チャートの作成

概要

・TrueFXにて公開されているミリ秒単位のクオートをもとに秒足データに加工する
・Pandasを用いて任意の間隔(今回は5秒間隔)のOHLCデータにリサンプリング
・matplotlibでロウソク足チャートを描画する

nysol前処理

為替データに関しては参考サイト1からダウンロード可能。
nysolは2.0を使用。

TrueFXのサイトからダウンロードできるデータは月次でcsvにまとまっているので、とりあえず日次に分解して個別ファイルに保存するようにしています。
私はこれをparseFxData.shというシェルスクリプトとして保存しています。

対象ファイルは2014年8月(古くてすみません)としています。適当に実行するとdailyというディレクトリが作成され、その中に各ファイルが生成されます。

150104_1

各csvのフォーマットは
USDJPY, yyyymmdd, yyyymmdd hh:mm:ss:sss, bid, ask
という感じになっています。
sssはミリセカンドです。どうでもいいですが、ミリセカンドってこういう場合どうやって記述するんですかね。

pandasリサンプリングとロウソク足チャートの描画

pandasにはデータフレームを任意の時間間隔でリサンプリングしてくれるメソッドがあります。(参考2)
が、ドキュメントには「これだけ書いとけばあとは分かるだろ」と非常に簡素なやっつけ仕様になっています。

11行目でread_csvを使用しcsvを読み込むまでは良いのですが、ロウソク足チャートを使用するためには時刻データをDatetimeIndex形式に変更しなければなりません。
が、この処理が非常に重いです。なので既に変数dfが定義されている場合にはその処理を割愛するようにtry exceptを使っています。
また、コメントにも書いてありますが、DataFrameのdtypeはfloatに指定しないと後々エラーになります(にっこり

22行目がリサンプリングを行なっているところになります。2014年8月1日の時刻が10時(つまり10時0分0.000秒から10時59分59.999秒まで)を対象に、ohlc形式でリサンプリングを行ないます。
脱線ですが、このhowの部分はfirst, last, min, max, sumなど基本的なものには対応しています。
closed=“right”で各バケットの右境界を含むようにしています。
label=“right”で、各バケットのindexとして、右境界の時刻が使用されます。(ex.0秒〜5秒のデータを対象とした四本値が5秒としてインデックスされる。closeが5秒なので個人的にはこっちのほうが読みやすい)
fill_method=“ffill”で、NAとなる部分は、最も最近の値で埋めます。こうすることで値動きが無い部分に関してはグラフが直線になります。

さらに23行目でインデックスに使用した時刻データを別個のカラムとして追加しています。
調べた範囲ではなぜかこうしないと後のロウソク足チャートに時刻データを渡すことが出来ませんでした。

26行目以降はmatplotlibで描画する方法。
ロウソク足のロウソクの部分の幅を力技で指定しています。どうやらこのwidthは単位がdayのようなので、5秒をday換算したのち、見やすくするために少し調整(candleBoxAdj)を行なっています。

で、でてきたものがこちら

150104_2

もう少し調整が必要な気もするが、まあとりあえず今日はこのへんで

参考

TrueFX — Tick-By-Tick Real-Time And Historical Market Rates, Clean, Aggregated, Dealer Prices
https://www.truefx.com/?page=frontpage

pandas.DataFrame.resample — pandas 0.15.2 documentation
http://pandas.pydata.org/pandas-docs/dev/generated/pandas.DataFrame.resample.html

メモ: multiprocessingを使うまで(Part 2)

概要

・クラスメソッドを並列化しようとすると謎のフリーズを起こすのでそれの回避法
・並列化する関数へ引数を複数個渡す方法
・前回(http://nekopuni.holy.jp/?p=1381)の続き

前回の問題点

前回のmultiprocessingの紹介で、何とか並列化が可能になりめでたしめでたしかと思われたのですが、クラスに組み込むと謎のフリーズが発生します。
具体的には適当なクラスメソッドを並列化する際に発生するようで、以下の様なコードを書いた時に発生します。

色々stackoverflowを徘徊すると同じような問題に陥っているケースが多いようなのですが、イマイチ自分用にまとまった回答が掲載されていなかったのでメモ。

Pythonコード

必要な手順を箇条書きにすると
・おまじないが必要(後述するmultiprocess_with_instance_method.pyをロードする)
 →5行目、これはただ単に参考サイトに載っている一部をコピペして外部ファイルとして保存しただけ

・クラスを定義するときにobjectと書く
 →7行目、どういうわけかここはobjectと書く必要があるようです。

・複数引数の渡し方もついでにメモ。これは相対的に簡単。
 →23行目。これはそのまんま。pool.map()でも複数引数は可能ですが、書き方が少しトリッキーなのでpool.apply_async()のほうが簡単。
となります。それぞれ対応する行をハイライトしておりますので、適宜ご参照ください。

だんだん’__’が顔文字に見えてくるとかいうどうでもいい情報はどうでもいいですね。

参考

Example showing how to use instance methods with the multiprocessing module
https://gist.github.com/fiatmoney/1086393

2015年のご挨拶

このサイトにどういう巡り合わせか訪れてしまった皆様あけましておめでとうございます。
特に改まってご挨拶する予定もなかったのですが、こうやって書いておけば後の自分が見返した時に少しは進捗が確認できるのではないかと考え急ごしらえで書いています。

ということで2014年初からのポストを見返してみると当時はまだ回帰分析が出来るようになったのに毛が生えた程度のものでした。
最近はPythonとRに触る機会が多いのでそれのメモや、ときどき趣味で進めてる強化学習のポストが挟まるという感じでしょうか。
自分だけのアルゴトレーディングシステムを構築し、寝てる間に収益がどんどん溜まっていく夢のメシウマライフにはまだまだ程遠いですが、精進していきたいと思います。

特にオチがあるわけでは無く、これからも毒にも薬にもならないような形で進むと思われる弊ブログですが、本年も何卒よろしくお願い致します。