Python x SPY500でシストレ作るならこんな構成

資金投資を考えた時に、プログラミングと機械学習の知識を使って一定自動化できないかな?と考えたので、考えている構成などをメモメモ。

最初は長期投資枠で資金の8割くらい当てておいて、残り2割でアグレッシブな運用利かせたい想定。 ↓Udemyの講座で勉強した内容を使って、実際に手を動かしたくなった。

理論面の勉強ばかりで、ちょっと飽きて来たので、Kaggleと2本立てで実際のデータ分析のトレーニングになればと。 www.udemy.com

概念マップ

対象とする金融商品

SPY500系のレバレッジ3倍ETF --> 市場全体が相手なので、騙しが少ないらしい。(酒井さん談) --> レバレッジが効いているのでボラは大きいが、利幅は取れる --> 追証の心配なし

ロードマップ

  • デモの開発
  • バックテストの実行
  • デモトレードを1ヶ月実行
  • 実戦投入!

アルゴリズムと基本戦略

  • 基本的なテクニカル指標 / 時系列(ARIMA・ARCH・GARCHあたり) / DLベースの予測を特徴量とした、2クラス分類を想定。以下のような本で、他に使えそうな特徴量があったら、そいつも使う。どなんだろ、ほとんどテクニカル指標が使われそうな気がする。

  • "〇〇以内に、損切りせずに、一定の利幅で利確可"が正解ラベルかな。

www.amazon.co.jp

  • トレンドフォローで順張り。確実に勝てそうな時に、コツコツ利確していく感じ。損切りは必ず設定する。

ja.wikipedia.org

運用している時のイメージ

  • 基本API Gatewayで叩き口を用意しておいて、GASで特徴量の計算と予測を実行。(GASじゃなく、Vueの方が楽しいのでそちらにするかも)
  • バックテストで有効なしきい値を設定しておいて、特定の値を超えた場合は、Slackに通知(# 買い時)

ライブラリ群

  • numpy / pandas / altair (基本ライブラリ)
  • statsmodels (時系列分析に使用)
  • XGBoost / LightGMB (最終的な予測に使用)
  • pyti(テクニカル指標の計算)
  • keras(ディープニューラルネットワークでの特徴量生成に使用。LSTM。)
  • backtesting.py (バックテストの実行)

kernc.github.io

バックテストのハイパーパラメーター

  • モデルの更新期間
  • 採用しきい値

残り色々

  • ETFをAPI経由で売買できるのか?? 契約してる証券会社によるが、国内証券会社だと難しい気がする
  • そもそもレバ3倍のSPY500 ETFをどこで買えるかまだ調べてない
  • 試験対策 & アプリケーション実装という意味で、実際に手を動かすフェイズになってきたかなという感じ。

【学習メモ】因果推論の介入と調整の図形的な解釈

前回の記事の続きです。自分流の理解がかなり混じっているので、あくまで参考程度に。 tokuchie.com

介入を考える際に、条件付き確率を図形解釈すると理解しやすいかなと思ったので、画像を残しておきます。 確率の図形解釈については、以下の書籍を参考にしています。

www.amazon.co.jp

問題設定

書籍のP3にある、新薬についての調査結果の集計表を用います。

薬投与 薬投与無し
男性 87人中81人が回復(93%) 270人中234人が回復(87%)
女性 263人中192人が回復(73%) 80人中55人が回復(69%)
合計 350人中273人が回復(78%) 350人中289人が回復(83%)

以下のグラフィカルモデルは以下になります。Z = 性別 / X = 薬投与有無 / Y = 回復有無となります。 この例ですと、性別(Z)が、X及びY両方に影響を及ぼしています。このようなデータからX->Yの影響を素直に分析しようとしても、Z->X->YのChainでの因果が入ってくるので、単純にはいきません。

こうした場合、Xをある特定の値に固定し、Z->Xの相関を無くしてしまう事で、純粋にX->Yの因果のみを明らかにする事ができます。

このXの値を固定する行為を介入と呼び、親に当たるZの影響を除外する事になるので、Zによる調整、Zのコントロールなどと呼びます。 要は、分析対象のXにインプットとして入ってくる有向線を削除したいという意図で、今回のように親要素(PA)が観測できている時には、Xに介入->PAによる調整というお決まりの形になります。

合わせて、図形も記してありますが、各変数間の従属関係を、"塗り方の歪み"で表現されるのがポイントです。

介入後は、Z->Xの相関が無くなるので、Z内部でのXの塗分けが均一になります。(今回の例ですと、全てのZでXが均一なので分かりにくいですが...)

f:id:osapii:20200623202135p:plain
グラフィカルモデルと図形的解釈

因果推論については、以下が同値である事を意識しながら、その時々に自分が理解しやすい方法で考えるのがポイントかなと思っています。

  • グラフィカルモデルにおいてX->Yで有向線が引かれている
  • YをXの関数で表す事ができる
  • YとXの相関係数が0ではない (ちなみに検定可)
  • XとYには因果がある
  • P(X|Y) = P(X)が成り立たない値の組がある
  • XとYは従属である
  • 図形的にZ毎にXで塗分けた際に、それぞれの面積が異なる

メモ

  • グラフで親要素が存在するときは、その要素に調整を加える(つまり、子要素に介入する)事で、子要素 -> 応答変数の因果を分析できる
  • PAが観測できていない場合にはバックドア基準により介入対象の変数を特定しますが、理解が進んだら記事にします。
  • Latex書くのしんどかったので、数式は飛ばしました。

バックドア基準はこのスライドがとっても分かりやすそう。感謝。

www.slideshare.net

スライドの講義verも滋賀大さんがアップロードしてくれていた。大感謝。 www.youtube.com

【学習メモ】d分離性について、実際の問題を解いて理解する

前回の記事の続き。

tokuchie.com

実際のグラフィカルモデルは、基本3類型には当てはめれないものの、逐次的に1つ1つの伝搬を見ていくと3類型で得た知識を使って、どんな複雑なモデルでも、従属性が存在するかどうか?を確かめる事ができる。というお話。

冷静に見ていくと、そこまで複雑では無いので、実際の例題を解いてみながら、d分離性についての理解を深めます!

www.amazon.co.jp

概念マップ

f:id:osapii:20200622163905p:plain

drive.google.com

d分離性とは?

とりあえず、書籍に記載してある定義を少し、わかりやすくして書きます。

(d分離) XがノードXとYの間の全ての道をブロックする時、Zが与えられた下でXとYはd分離されている。すなわちZが与えられた下でXとYは条件付き独立である。

そして、道pがノードの集合Zによりブロックされている事は以下と同値である。

  1. pは連鎖A->B->Cまたは、分岐A <- B -> Cを含み、Bが条件付けされている
  2. pは合流 A->B<-Cを含み、Bが条件付けされていない。さらにいかなるBの子孫も条件付けされていない。

何やら一見複雑そうですが、重要な点は以下です。 1. 条件付けされているかいないかで、従属性が流れるかどうかが変わる 2. 全ての道を確認し、1つの道であっても従属性が流れるならばその場合はd連結

上記の定義に従って、2つのノード間がd分離であるかどうかを確かめるには、3類型で確認した規則を、逐次的に経路に対して適用することになります。

具体的な例題1

実際の問題を解いてみます。

f:id:osapii:20200622164923p:plain 今回はこのようなグラフィカルモデルを題材に、ZとYがd分離であるかどうかを考えていきます。 Wには、特に何も条件を課していないことに注意します。 f:id:osapii:20200622164937p:plain まずは、ZとYの間に存在する経路を全て洗い出します。今回の例題ですと、Z->W->X->Yが唯一の経路となるので、この途中でブロックされるかどうかを逐次的に確認していきます。 f:id:osapii:20200622165222p:plain まずは、紫の円で囲まれた部分の伝搬を見ていきます。ここで、前回の記事で取り扱っている、合流点の形をしていますが、この時ZとXは独立となるので、Z->X間で従属性は流れません。

d分離の定義に戻ると、2に当てはまるのでブロックされている事になります。この時点で、ZとYがd分離である事が分かります。(色々な言い換えをすると、ZとYは独立 / ZとYの相関係数の絶対値は0に近い / ZとYは全く関係ないよと言う事です。)

具体的な例題

さて、同じグラフィカルモデルですが、今度はWに条件を課してみます。つまり、"Wが〇〇になる"と分かっている状況下において、ZとYに従属性が認められるか?を考えていきます。 f:id:osapii:20200622165713p:plain 今回は、条件付けされているので、Wを赤い点で示しています。 f:id:osapii:20200622165723p:plain 伝搬のステップは、先ほどの例題と同じですが、条件付けによって、従属性の流れが変わってきます。

f:id:osapii:20200622165733p:plain この時、Z->Xは従属の関係があるので、従属性が流れる事になります。(画像中では太い黒矢印で示しています。) Wが判明しているという状況下では、Zを知る事で、Xの値がある程度分かる。つまり従属性が認められます。具体的な例は以下です。

*1

f:id:osapii:20200622165745p:plain 次のステップに進みます。(オレンジの円) ここは単純に因果の関係が認められるので、従属性が流れます。 これで、Z->W->X->Yがブロックされる事なく繋がり、ZとYがd連結している事が分かりました。

メモ

  • 従属性が流れる部分は、ニューラルネットの逆伝搬と近いものがあって共通点を感じられた。同じように連鎖率を使う事で、ZがYに与える影響なんかも定量的に把握できるようになっているのだと思う。今後出てきそう。
  • d分離性の概念を用いる事で、グラフィカルモデルが与えられた時に、変数間の独立性を判断する事ができる。
  • ただし、重要なのはそのグラフィカルモデルが正しいかどうか?なので、その辺りの検証方法はこの先に出てくる、多分。

*1:音楽の成績がずば抜けて成績が優秀である(X) OR 運動の成績がずば抜けて成績が優秀である(Y)場合に、奨学金が与えられる(Z)とする。 この時、Zを観測する前であれば、ある学生が音楽の成績が良くても、運動の成績を推測する要素にはならない。(周辺独立が成立している)

一方で、奨学金を受け取ったという事象を観測した際には、音楽の成績が低い(X) -> 運動の成績がずば抜けて優秀(Y)などというように、Xを知る事で、Yの値を当てることが出来る。つまり、Zの下で、XとYは条件付き従属である。

【学習メモ】グラフィカルモデルの基本3類型のまとめ

こちらの書籍で、因果推論を学習したので、内容のメモ。 www.amazon.co.jp

パス図(グラフィカルモデルとほぼ同義なのかな??)やSCM(構造的因果モデル)については、こちらの書籍を読んだ事があったが、新しい気付きが色々あって面白い。特に、グラフィカルモデル単体で、関数の具体的な数式記述が無くても、現象のメカニズムを掴むのに有益だという部分は興味津々で、基本3類型はその基礎となる部分。 www.amazon.co.jp

概念マップ

f:id:osapii:20200622151903p:plain

drive.google.com

  • Ux / Uy / Uzは、それぞれの外生変数を表す。外生変数とはモデル内部でメカニズムが説明されない、外部から与えられる所与の変数の事。モデル内では説明されない、"諸々いろんなノイズを全てまとめたもの"としての誤差項は、外生変数として与えられる。

  • 通常のモデルは、各ノード間での経路が唯一であるという事は無いので、上記のグラフィカルモデルを単純に当てはめる事で、因果を分析する事は無い。(が、工夫をする事で当てはめられる。keyword: d分離性)

連鎖経路 (chain) のグラフィカルモデル

f:id:osapii:20200622152426p:plain

当てはまる特徴

  1. ZとYは(おそらく)従属である
  2. XとYは(おそらく)従属である
  3. ZとXは(おそらく)従属である
  4. ZとXは、Yの下で条件付き独立である ※Ux/Uy/Uzがそれぞれ独立である事を仮定

分岐経路 (fork) のグラフィカルモデル

f:id:osapii:20200622152416p:plain

当てはまる特徴

  1. XとYは(おそらく)従属である
  2. XとZは(おそらく)従属である
  3. YとZは(おそらく)従属である
  4. YとZは、Xの下で条件付き独立である ※Ux/Uy/Uzがそれぞれ独立である事を仮定

合流点(collider) のグラフィカルモデル

f:id:osapii:20200622152441p:plain

当てはまる特徴

  1. XとZは(おそらく)従属である
  2. YとZは(おそらく)従属である
  3. XとYは独立である
  4. XとYは、Zの下で条件付き従属である ※Ux/Uy/Uzがそれぞれ独立である事を仮定

メモ

  • 従属は独立で無い事を指し、相関がある状態を意味する。日本語の字面的になんとなく因果関係を想像してしまいそうになるので、その点注意。

VScodeでAnaconda仮想環境の実行とデバッグを行う時のチェックポイント

PyCharmや、JupyterNotebookでコードを書くことが多かったが、やはりエディタを1つに統一したい欲があり、VScodeに移行する事に。

  初め、Import Errorが頻出して非常に手間取ったので、備忘録的に対応内容を記載しておく。

手こずったのはそもそもエディタの設定に対する理解がザルだったのが大きいので、当たり前だろ!という事も書いてあると思います。笑

今回色んな記事をググって時間を浪費してしまったが、恐らく以下の対応は全て行えていれば問題ないはずです。

自分の場合は、これでちゃんと動くなるようになりましたというリストなので、もしかしたら過剰な設定内容も含まれているかもしれません!

 

バージョン

VScode April 2020 (version 1.45)

Anaconda 2020.02-Windows-x86_64

記事内の仮想環境名 zero_to_deep_2

公式ドキュメント

Using Python Environments in Visual Studio Code

①ターミナル上で仮想環境に切り替えているか?

まず前提として、VScode内では①実行環境と②実行用のPythonの参照先③デバッグ用のPythonの参照先の3つを設定する必要がある。 実行環境は今回仮想環境と想定しているので、VScode内のターミナルで以下のコマンドを入力して、あらかじめ仮想環境に入っておく。

activate your_env_name

f:id:osapii:20200610091749p:plain

この時、VScodeの統合ターミナル内で、正しく仮想環境に切り替えできているか確認する。(画像赤線の部分)

②Pythonの実行環境のPathが通っているか?

続いて、実行環境のPathを通していく。設定から行っていく。(CTRL + , ) ここに今回使用したい仮想環境のPythonPathを記述する。これで実行時に該当のPathのPythonをVScodeが参照してくれる。 f:id:osapii:20200610092229p:plain

後は、condaPathにAnacondaのパスを設定する箇所もある。こちらも、必要かどうか分からなかったが一応設定しておいた。 f:id:osapii:20200610092502p:plain

ここで重要なポイントとして、設定はワークスペース設定の方に記述するという事。VScodeの設定内容の優先度は、 ○○.json > ワークスペース設定 > ユーザー設定 となっている。

実際、UIからワークスペース設定を変更すると、setting.jsonが更新される仕組みになっている。

別の仮想環境で実行したいよ~という場合は、別のワークスペースを新規に設定してpythonPathを設定する必要があると思われる。 f:id:osapii:20200610092922p:plain

ちなみに、Python拡張機能を導入するとエディタ下部のバーから実行環境を変更できるが、これは正に、ワークスペース設定のpythonPathを変更している。なので、実際はパスを調べて入力しなくても、バーから任意の仮想環境を選択する事で事足りるはず。 f:id:osapii:20200610093358p:plain

↑①と②が正しく済んでいれば、Pythonファイルの実行が通るのではないか?と思う。ちなみに、この状態だとデバッグは動かない。

③Pythonのデバッグ環境のPathが通っているか?

最初に書いたとおり、VScodeでは実行環境とデバッグ環境は扱いが異なる。デバッグ環境の設定はlaunch.jsonに記載していく。

ここでは、実行環境と同様にpythonPathを変更すればOK。②が済んでいれば、setting.jsonに正しいpythonPathが入っていると思うので、そいつをコピペしておく。

↓今、手元で動いているデバッグ設定を貼っておくとこんな感じです。

{
    // Use IntelliSense to learn about possible attributes.
    // Hover to view descriptions of existing attributes.
    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
    "version": "0.2.0",
    "configurations": [
        {
            "name": "Debug",
            "type": "python",
            "request": "launch",
            "program": "${file}",
            "console": "integratedTerminal",
             # ここにpythonのパスを書こう
            "pythonPath": "C:\\Users\\m_osa\\anaconda3\\envs\\zero_to_deep_2\\python.exe"
        }
    ]
}

これでデバッグ環境の設定が終わったので、F5でのデバッグも機能するようになる。はずだが、自分の場合はここでつまづいた。

f:id:osapii:20200610093957p:plain
やったー

④VScodeを起動してからAnacondaを読み込むまで十分に時間が経過しているか?

①~③まで完了したのに、なぜかImport Errorが発生。30分ほど設定変更と再起動を繰り返していましたmm

色々試してみると、最初の読み込みに時間がかかるようで、エディタの立ち上げ直後は実行もデバッグも正しく環境を参照してくれません。

大体1分程待つと、設定した内容で動作してくれます。結構忘れがちなポイントだと思うので気を付けて下さい。

おわり

環境構築に時間取られている時、悲しい気持ちになる( ;∀;)

検出力(検定力)の勉強をする時に、初めに知っておきたい事!

検定力の勉強をしているときに、こんなイメージで理解しておくと分かりやすいかもと思ったので、メモ的にまとめておきます。 以下、単純化のために分散既知の母平均の片側検定を行うとして書きます。実際は、検定の種類によって検出力及び効果量の算定方式が異なりますが、根本的なアイデアは同じです。

モチベーション

  • 検定力(検出力)周辺の概念を抑える事で、普段の仮説検定の質をより高める事ができると思っておけば良い
  • 帰無仮説を採択する場合(棄却できなかった場合)は、ただ棄却できませんでした!ではなくて、"少なくとも○○以上の差は無いようです!"と、範囲を示すことができる
  • 対立仮説を採択する場合(棄却できた場合)は、ただ棄却できました!ではなくて、"〇〇以上の差はあります!"と、影響度の幅を示すことができる

小難しそうな様相をしているが、ツールボックスとして使えるようにしておくと、役に立つ概念だと思う。ビジネスの現場では、ROIを検討する場面が多く、検出力の概念を考慮する事により合理的な意思決定に貢献できる場面がありそう。

検定力について知っておきたい事柄

学習するときにこのイメージを持っておくと、諸々スムーズに行きそうです!

1.仮説検定はサンプルサイズを増やすと、どんな仮説でもカンタンに棄却できてしまう

母平均を μ、母分散(既知とする)を σ^{2}、標本平均を \overline{X}、サンプル数を Nとおくと、母平均の検定で用いる統計量  zは以下の式で表すことができる。

 z = \frac{\overline{x} - μ_0}{\sqrt\frac{σ^{2}}{N}}

この式の分母に注目し、Nを ∞にした時、つまりサンプル数を数えきれないほど大量に集めた場合を考えてみる。この時分母は0に収束するので、 z全体は非常に大きい値に落ち着く。 そして、統計量zの値が大きくなると、棄却域に達し棄却されることになる。ここで、問題なのは観測した標本平均の\overline{X}値に関係なく、仮説が棄却されてしまう事。

現象を理解するために、わざわざ検定を行っているのに、現象と全く関係のないサンプルサイズの大きさという要因で、検定結果が決まってしまう。これではまずい。

2. 検定力は、【ある特定の】対立仮説に対する議論

検定力 というワードが出て来た時点で、ある特定の対立仮説を想定している事をイメージしておく。母平均の検定の例で言えば、対立仮説は μ = 34 μ = 42など、具体的な値が入った上で、それぞれの対立仮説にたいする検出力の議論が成り立つ。 関数の形で書くと、 検出力 = f( μ = 42の対立仮説)といった雰囲気。

という訳で、現在行っている仮説検定一般に対する検定力というものは存在せず、ある特定の対立仮説に対して定義される量である事を意識する。

3. 有意水準は"帰無仮説の元で仮説を棄却する確率"、検定力は"ある特定の対立仮説の元で仮説を棄却する確率"

例えば、とある帰無仮説を棄却できたとする。この時もし本当は対立仮説(例えば、 μ = 30)が正しかったとしたら、その対立仮説が成立している世界の元で帰無仮説を棄却できる確率はどれくらいだろう?←この確率が検定力の正体です。 検定力という言葉を、確率に無理やり頭の中で変換するのがポイントだと思います。

そして、1-検定力は、 第二種過誤と定義されます。先に、第二種過誤から入って、1-第二種過誤 = 検定力 とするテキストもありますが、理解しやすい方で覚えておけば良いかと思います。 上記の例、仮に μ = 30の検定力が0.80だったとしましょう。逆に言うと、0.20(20%)の確率で、帰無仮説を棄却できません。←本当は対立仮説を採択すべき場面で正しく採択できない事を意味するので第二種過誤と言われます。

↓具体的なケースで考えてみると以下のイメージです。

①帰無仮説を棄却できず、 μ = 30の検定力が十分に高いとき

 μ = 30の世界だと高い確率で棄却できるのにも関らわず棄却できなかった ---> 少なくとも、 μ = 30に対する差分よりも大きな差は無いだろう。

②帰無仮説を棄却できず、 μ = 30の検定力が低いとき

 μ = 30の世界でも棄却できない可能性が高い --->  μ = 30に対して、確かな情報は得られない。

③帰無仮説を棄却、 μ = 30の検定力が十分に高いとき

 μ = 30の世界だと高い確率で棄却できる&帰無仮説が棄却できた --->  μ = 30の元でも棄却される確率が高いので、この程度の差分は存在しているだろう。

④帰無仮説を棄却、 μ = 30の検定力が低いとき

 μ = 30の世界だと棄却できない可能性が高いのに帰無仮説が棄却できた ---> 検出力が足りていないので、 μ = 30の対立仮説が成り立っているかどうかは不明慮

-->もし、 μ = 30の対立仮説を評価したいならば、十分な検定力を獲得できるまでに、サンプルサイズ数を増やす必要がある。

途中で、差分というワードが出てきましたが、実際は効果量という単位で表します。"どれくらいの差があるか?"は検定毎に単位が異なりますし、母集団の分散によっても評価が異なるので、離れ具合を定義する指標として効果量があるというイメージです。

4. "ある特定の対立仮説の元で仮説を棄却する確率"を計算する橋渡しは、観測値

以下、抽象的ですが、記載しておきます。帰無仮説の元で、統計量 zが棄却される可能性というのは、標準正規分布の棄却域を例えば以下のような不等式で表されます。  1.96 > z = \frac{\overline{x} - μ_0}{\sqrt\frac{σ^{2}}{N}}

観測値 \overline{x}を未知数として展開していくと、 \overline{x} < 25などの式に書き直すことができます。つまり、"観測値が○○以下だったら帰無仮説を棄却しますよ~!”という範囲です。 ↑は式変形をしているだけなので、観測値が○○以下になる = 帰無仮説を棄却する という関係になっていることに注意してください。

後は、対立仮説の元で計算した  zの元で、観測値が○○以下になる確率を計算する事で、求めたい "その対立仮説が成立している世界の元で帰無仮説を棄却できる確率"を計算する事ができます。 このように、観測値をハブにして、計算を行うイメージを持っておくと何をやっているのか理解しやすいかと思います。

おすすめの学習ソース

私もまだ学習中の身ですが、記事を書くにあたって参考にした学習教材を置いておきます。今後、ExcelやPythonを使った具体的な検出力計算も記事にしたいと思います。

www.amazon.co.jp

2010 年度前期 情報統計学 第12回 検定と検出力 http://racco.mikeneko.jp/Kougi/10a/IS/IS12pr.pdf

メトロポリス・ヘイスティング法をPythonで簡単に実装してみる

メトロポリス・ヘイスティングス法とは?

ja.wikipedia.org

  • ある特定の確率分布から、その確率密度の大小にしたがってサンプル値を抽出するための基本的なアルゴリズムの一種
  • サンプリング自体は、ベイズ統計で、事前分布と尤度の掛け合わせからなる事後分布からサンプル抽出を行いたい時によく用いられる。
  • しかし、一般にメトロポリス・ヘイスティングス法によるサンプリングは収束に時間がかかるという課題があるため、改良版として、ハミルトニアンモンテカルロ(HMC)という改良版がある。しかし、基礎は共通。

サンプリング手法は沢山あるが、メトロポリス・ヘイスティング法は諸々の手法のベースになっているので、抑えておく必要がある。

ステップ数が無限ならば、正確に対象の確率分布をトレースしたサンプリングを実施できる事は理論的に証明されているらしいです。

アルゴリズム のアイデア

直感的に書くと以下のような仕組みでサンプリングを行う事で、対象とする確率分布に従うサンプル群を得ることができます。

対象とする確率変数を  x とし、確率分布を  f(x)とおきます。また、サンプリング結果を  sampleと表します。(例えば、2番目のサンプリング結果は  sample_2と表す事にします。)

STEP1: 初期値 x_1を決め、確率密度の値  f(x_1)を求めておきます

STEP2:  x_1に、ランダムな値を足し、遷移先の値候補 x_2及び、確率密度の値  f(x_2)を求めておきます。

STEP3:  f(x_2) > f(x_1)ならば、 sample_2 = x_2とし x_2をサンプルとして採用(受容と言います。) 一方で、  f(x_1) > f(x_2)ならば、確率密度の比率 \frac{f(x_2)}{f(x_1)} の確率で受容しますが、それ以外の場合は sample_2 = x_1として、既存の値をサンプル値として採用します。

STEP4: STEP2に戻り、決められたステップ数の数だけSTEP2->3を繰り返す。

現在位置よりも、移動先の確率密度が低い場合であっても、一定確率で受容する事で、確率分布に沿ったサンプリングを実行できるのが肝です。より詳しい事を知りたい場合は、以下の書籍や動画が理解の助けになるかと思います。

www.youtube.com

Pythonで実装してみる

サンプリングを実行する関数metropolis()の定義

実装したコードは以下になります。なお、本コードは先に紹介している書籍の実装を参考にしています。というよりかは、ほぼそのままなので、是非本書をご覧になってください。

def metropolis(func, steps = 10000):
    # サンプリング結果保存用のリスト
    samples = np.zeros(steps)

    # STEP1: 初期値の設定 及び 確率密度の計算
    # 与えられた確率分布の期待値で初期値に設定
    old_x = func.mean()
    old_prob = func.pdf(old_x)

    # stepsの数だけループ
    for i in range(steps):
       # STEP2: 正規分布に従う乱数を足して、新しい遷移先候補を決定
        new_x = old_x + np.random.normal(0, 0.5)
       # 遷移先候補の確率密度を計算
        new_prob = func.pdf(new_x)
       # STEP3: 確率密度の比を計算
        acceptance = new_prob / old_prob
       # acceptanceを基に受容 or 棄却を決定
        if acceptance >= np.random.random():
            samples[i] = new_x
            old_x = new_x
            old_prob = new_prob
        else:
            samples[i] = old_x
    return samples

文章で書くとややこしそうでしたが、とても短いコードでサンプリングを実行する関数を記述できました。acceptanceの判定に、np.random.random()で[0~1]の乱数を用いている点がポイントかと思います。

STEP2では、移動先の提案に正規分布に従う乱数を用いていますが、この場合は、正規分布を提案分布として採用している事になります。提案分布の選択は、幾分かサンプリング結果にも影響を与える事になります。

サンプリングの実行

サンプリング対象となる確率分布を設定 -> 関数に渡してあげます。

func = stats.beta(1.2,3)
samples = metropolis(func, steps = 10000)

今回の例では、 α = 1.2  β = 3のベータ分布からのサンプリングを実行しています。stats.~で定義する確率分布を正規分布や、ポアソン分布やその他の分布に変える事で、同様にサンプリングを行う事ができるでしょう。

次は?

より効率良くサンプリングを行うための、HMC法やNUTSなどのアルゴリズムについても書きたいと思います。