読者です 読者をやめる 読者になる 読者になる

ともにゃん的データ分析ブログ

勉強したことの備忘録とかね

ハミルトニアン・モンテカルロ(HMC)法のざっくりとした解説とPythonによる実装

python ベイズ統計学

ベイズモデリングが流行っている中で多くのRユーザーはStanを使って解析をしているんではないかと思います。そして、Stanはハミルトニアンモンテカルロ(HMC)法と呼ばれる方法で事後分布からのサンプルを得ています。色々と解説記事はありますが、超ざっくりとHMCの原理をメモとして残しておくことにします。

ここで、基本的に私はHMCを伊庭先生のハミルトン4.pdf - Google ドライブこの資料で勉強したので、伊庭先生の資料を読めば私のこの記事は必要ないことをあらかじめ断っておきます((´^ω^)
なんでこれでいいの?という疑問は伊庭先生の資料で解決することでしょう((´^ω^)


記号の定義

まず、以下のように記号を定義します:

X: データ
\theta=(\theta_{1},\ldots,\theta_{d})^{'}: d次元パラメータベクトル。'は転置を表します。
r(\theta|X): 事後分布(後にpが別に出てくるので、ここではrとしました)

上記の事後分布からHMCでサンプリングをすることを考えます。


ハミルトニアンとハミルトン方程式

まっっっったく詳しくないのですが、ハミルトニアンとは運動エネルギーV(p)とポテンシャルエネルギー(位置エネルギー)U(\theta)の和で定義される以下の物理量のことです。ここでp=(p_{1},\ldots,p_{d})^{'}\thetaと同じd次元のベクトルとします:


H(\theta, p) = V(p) + U(\theta)


また、ハミルトン方程式と呼ばれるものが以下のように表現されます:


\begin{eqnarray}
\frac{d\theta_{j}(t)}{dt} &=& \frac{\partial V(p)}{\partial p_{j}}\\
\frac{d p_{j}(t)}{dt} &=& -\frac{\partial U(\theta)}{\partial \theta_{j}}
\end{eqnarray}

物理の言葉を借りると、pは運動量で\thetaはある物体の位置を表しています。上記のハミルトン方程式はある物質に運動する力pを与えてやることで、その物質の位置\thetaが変化するという現象を表したモデルになっています。このアナロジーは以下の説明でHMCを理解する上でとっても役立ちます。


さて、このとき以下のようにp\thetaの同時確率密度関数q(p,\theta)ハミルトニアンH(\theta, p)だけの関数だけの関数、つまり

 
q(p,\theta) = f(H(p,\theta))

という形で定義したとき、上記のハミルトン方程式に従ってpの値と\thetaの値を変化させれば、それは確率密度関数 q(p,\theta) = f(H(p,\theta))からのサンプリングをしたこととみなせるという素晴らしい性質を持っています。ここがHMCのポイントだと思います。

ハミルトニアン方程式の解き方

一般に、リープ・フロッグ法と呼ばれる方法で解くようです。ここで次のように記号を導入します:

p(t)t番目のステップにおけるpの値
\theta(t)t番目のステップにおける\thetaの値

また、次のようにパラメータを導入します:

T:状態変化を行う回数
\epsilon:1ステップで状態変化をどの程度許すか



詳しい解説はグーグル先生にまかせて、以下に計算の流れを記します:

===========================================================

Step 0. t \leftarrow 0、またp(0)\theta(0)に初期値を設定。

Step 1. 以下のように\thetaを更新:
\theta_{j}(t+0.5) \leftarrow \theta_{j}(t) - \frac{\epsilon}{2}\frac{\partial U(p)}{\partial p_{j}(t)}

Step 2. 以下のようにpを更新:
p_{j}(t+1) \leftarrow p_{j}(t) + \frac{\epsilon}{2}\frac{\partial V(\theta)}{\partial \theta_{j}(t)}

Step 3. 以下のように\thetaを更新:
\theta_{j}(t+1) \leftarrow \theta_{j}(t+0.5) - \frac{\epsilon}{2}\frac{\partial U(p)}{\partial p_{j}(t+1)}

Step 4. t \leftarrow t + 1 としてStep. 1 - Step. 4をt=Tになるまで繰り返す。

===========================================================


ここで、t+0.5ステップという普段見ない表記を導入していますが、中間のステップということを表現したいがために便宜的に導入したものです。


具体的なq(p,\theta)の形

上記のハミルトニアンの右辺に出てくるV(p)U(\theta)を、具体的に以下のように指定してやります:

\begin{eqnarray}
V(p) &=& \sum_{i=1}^{N}\frac{p_{i}^{2}}{2\tau^{2}}\\
U(\theta) &=& -\log{r(\theta|X)}
\end{eqnarray}


また q(p,\theta) = f(H(p,\theta))fについて、指数関数の逆数

f(y) = \exp(-y)
を指定してやります。


このときq(p,\theta)は以下のように:


\begin{eqnarray}
q(p,\theta) &=& \exp\left(-\sum_{i=1}^{N}\frac{p_{i}^{2}}{2\tau^{2}} + \log{r(\theta|X)}\right)\\
&=& \exp\left(-\sum_{i=1}^{N}\frac{p_{i}^{2}}{2\tau^{2}}\right)\exp\left(\log{r(\theta|X)}\right)\\
&=& \left(\prod_{i=1}^{N}\exp\left(-\frac{p_{i}^{2}}{2\tau^{2}}\right)\right)r(\theta|X)\\
&=& \left(\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi}\tau}\exp\left(-\frac{p_{i}^{2}}{2\tau^{2}}\right)\right)r(\theta|X)\\
\end{eqnarray}


上記の分布をよくみてください。平均0、分散\tau^{2}正規分布と目的の事後分布の掛け算になっています。掛け算で表現されるということは、つまりp\thetaは独立しているということです。

さて、q(p,\theta)ハミルトニアン方程式の状態変化の定常分布でした。加えて、上記の式をみるとpは平均0、分散\tau^{2}正規分布に従います。つまり、\thetaとは独立にpを平均0、分散\tau^{2}正規分布からサンプリングし、その値を用いてハミルトン方程式を解くことで新たな\thetaの値を得ます。

ここで、ハミルトン方程式と解く際に、現在の\thetaの値から状態を変化させるために以下のようにリープ・フロッグ法のStep 0. を変更します:

===========================================================

Step 0.
t \leftarrow 0と設定
p(0)に平均0、分散\tau^{2}正規分布から発生させた値を設定
\theta(0)を現在の\thetaの値を設定

Step 1. 以降は上で述べた方法と同じで、ハミルトン方程式を解いて新たな\thetaを得ます。

===========================================================

以上の流れで得た\thetaを目的の事後分布からのサンプルとみなしたいところですが、上記のハミルトン方程式をリープ・フロッグ法を使って得ということは、コンピュータで”近似的”に解くということで、厳密解ではありません。そこで、ハミルトン方程式を解いて得た\thetaを、目的の分布からのサンプルの”候補”とみて、メトロポリス法のように以下のようなサンプルの受理・棄却のステップを設けてやります。ここで\thetapを現在の\thetaの値、\theta^{\ast}p^{\ast}をハミルトン方程式を解いて\thetapの状態変化を行った後の値とします:

===========================================================

確率 min\left(\frac{\exp(H(p^{\ast},\theta^{\ast}))}{\exp{H(p,\theta)}},1\right) で状態変化後の値\theta^{\ast}に移動

===========================================================


このように、ハミルトン方程式を決定論的に値を決定するパートと上記の受理・棄却パートのように確率的に値を決定するパートがあるので、HMCはハイブリッド・モンテカルロと呼ばれたりもしています。

HMCのパラメータ

今まで見てきた中で、分析者が指定しなければならないパラメータが以下のように3つあります:

\tau:運動量pが従う正規分布の分散
T:リープ・フロッグ法でハミルトン方程式を解くときの状態変化の回数
\epsilon;リープ・フロッグ方でハミルトン方程式を解くときの1ステップあたりの状態変化の大きさ


\tauを大きくすれば0から離れた値のpが出やすくなり(物理的には\thetaの運動が大きくなり)、それにともなって\thetaの値は大きく変わりやすく(物理的には\thetaの位置が大きく移動しやすく)なります。Tの値を大きくすれば、状態変化をする回数が多くなるということなので、前の値の\thetaから異なる値へより変化していきます。\epsilonの値を変更すると、状態変化の具合が変わります。


Stanに実装されているNUTSという手法は、上記のパラメータTを自動化するもので、HMCの敷居を大きく下げました。


Pythonによる実装

PythonでHMCを実装します。
今回は、簡単に平均100、分散2^2正規分布からサンプリングをしたいと思います。

pythonのコードは以下のようになりました:

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Log probability
# ここでは平均mu, 標準偏差sigmaの正規分布
# -1をかけると物理でいうポテンシャルエネルギー
def log_normal(x, mu, sigma):
    return -0.5*np.log(2*np.pi*sigma**2) - (x-mu)**2/(2*sigma**2)

# Log probabilityの導関数
def d_log_normal(x, mu, sigma):
    return -(x-mu)/sigma**2

# 運動エネルギー
def momentum(p, tau):
    return p**2/(2*tau**2)

# 運動エネルギーの導関数
def d_momentum(p, tau):
    return  p/tau**2

# ハミルトニアン
# 運動エネルギーとポテンシャルエネルギーの和
def Hamiltonian(x, p, tau):
    global mu, sigma
    
    return momentum(p, tau) + (-1.*log_normal(x, mu, sigma))

# リープ・フロッグ法を事項する関数
def proceed_leapflog(epsilon, x, p, tau):
    global mu, sigma
    
    x += -0.5*epsilon*(-1.*d_momentum(p, tau))
    p += epsilon*d_log_normal(x, mu, sigma)
    x += -epsilon*(-1.*d_momentum(p, tau))
    
    return x, p

# HMCを1ステップ実行するサブルーチン
def proceed_HMC_iteration(x, tau, epsilon, T):
    p = np.random.normal(0, tau, size=1)[0]
    p_new = p
    x_new= x
    for t in range(T):
        x_new, p_new = proceed_leapflog(epsilon, x_new, p_new, tau)
        
    alpha = np.exp(Hamiltonian(x, p, tau) - Hamiltonian(x_new, p_new, tau))
    u = np.random.uniform()
    if u < alpha:
        x_accepted = x_new
    else:
        x_accepted = x
        
    return x_accepted

# HMCを実行する関数
def proceed_HMC(tau, epsilon, T, ite, init):
    # Initialize
    x = [init]
    for i in range(ite):
        x.append(proceed_HMC_iteration(x[i], tau, epsilon, T))
    
    return x

以下のように目的の正規分布の平均を初期値と設定して実行してみます。

init = 100
theta = proceed_HMC(tau=2, epsilon=0.1, T=10, ite=2000, init=init)

plt.plot(theta)

f:id:kefits:20170314002838p:plain

ちょっと自己相関が高い感じになってるのでリープ・フロッグ法でハミルトニアン方程式を解く際のステップ数Tを30にして、もっと\thetaの状態を変化を待ってみましょう。

init = 100
theta = proceed_HMC(tau=2, epsilon=0.1, T=30, ite=2000, init=init)

plt.plot(theta)

f:id:kefits:20170314003947p:plain

自己相関が小さくなった気がしますね。

ここでちょっと意地悪をして、初期値の値を目的の分布のテールがほとんどない0からスタートしてみます。

init = 0
theta = proceed_HMC(tau=2, epsilon=0.1, T=30, ite=5000, init=init)

plt.plot(theta)

f:id:kefits:20170314004610p:plain

すごいぜHMC。メトロポリス・ヘイスティングではこうはいかないと思います。


最後に

ベイズ推論をするにあたって、多くの場合で事後分布が解析的に解けないために事後分布に従う乱数(事後サンプル)を発生させ、それ使用して色々な分析をします。よく聞くマルコフ連鎖モンテカルロ(MCMC)法は、事後サンプルを得る方法群の名前です。ギブスサンプラーメトロポリス・ヘイスティング法といった色々な種類のMCMC法が存在します。ベイズ推論において、Stanに実装されているHMCも事後サンプルを得るために使われる一つの方法です。要はどんな方法であれ、最終的に手元に"理論的に正しい"事後サンプルが手に入っていれば解析ができるので、それでよいのです。ただしそれぞれの方法には特徴があって、問題に合わせて使い分ける必要があるように思います。

世の中にはStanやBUGS、JAGSといったベイズ推論を実行することのできるとっても便利な言語があって、我々はそれらを使ってすぐにベイズ推論を行うことができるようになっています。モデリングの試行錯誤も楽ちんです。一から実装していては、モデルが変わるたびにコードを書き直さなければならないので試行錯誤が大変です。ただし、例えば最終的に構築されたモデルをに対してギブスサンプラーが構築できる時にはギブスサンプラーを使ったほうがサンプリング時間が圧倒的に短く、収束も早いことがあります。昔、以下の本の内容をStanで実装していて思い知らされました(この本の著者はマーケティング分野ではかなり有名な研究者です)。

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)

特定のモデルについてですが、Rossiは線形代数や分布の特性を駆使してギブスサンプラーによって効率的にサンプリングできる関数をパッケージ化して{bayesm}というパッケージを作っています。Stanとは比べ物にならないスピードで収束します。
Stan便利だぜ、すげえぜ、と思いますが、他の方法の方が良いこともあるということは頭に入れて置かなければなりません。

【python】beautifulsoupでYahoo! ファイナンスから日経平均のデータをスクレイピング

python

ずっっっっっと前にbeautifulsoupでスクレイピングしたことがあったけど、使い方を完全に忘れてたので再び入門的なことをやってみた。

とりあえずYahoo! Financeから日経225に関するデータを引っ張ってこようかと思います。

from bs4 import BeautifulSoup
import urllib.request
import time

# 今回は30ページ分のデータを取得してみる。
page_num = 30
stock_temp = []
for i in range(page_num):
    # Yahoo Financeのページ。url末尾の数字を変更すると日経225の過去のデータが取得できる。
    url = "http://info.finance.yahoo.co.jp/history/?code=998407.O&sy=2010&sm=12&sd=4&ey=2017&em=3&ed=4&tm=d&p=" + str(i+2)
    html = urllib.request.urlopen(url)
    soup = BeautifulSoup(html, "lxml")

    # 上記urlのソースをみると<td>~~~</td>にほしい数値が入っているっぽいから、soup.find_all("td")でその部分を抽出する。
    # soup.find_all("td")では<td>~~~</td>といったタグと一緒にリスト型で結果を抽出してくるので、
    # リストのそれぞれの要素に対してget_textメソッドを使って数値だけにする。
    stock_extract = [value.get_text() for value in soup.find_all("td")[3:103]]
    stock_temp.extend(stock_extract)
    
    time.sleep(0.5)
    
stock_temp = np.array(stock_temp)


stock = stock_temp.reshape(int(len(stock_temp)/5), 5)
stock = pd.DataFrame(stock[:,1:5], columns=["start", "high", "low", "end"], index=stock[:,0])

# 株価のカラムが文字列になっていて、かつカンマが入っているのでカンマを除去してfloat型にする。
for i in range(4):
    stock.ix[:,i] = stock.ix[:,i].str.replace(",", "").astype(float)


結果として

stock

f:id:kefits:20170305123205p:plain

こんな感じ。いいんじゃないですかね?

【R】ダミー変数を一度に生成する関数

R package 【R】自作関数

探せばあるんだろうけど、データフレームを引数に複数列を一度にダミー変数化する関数を作りました。
よければ使ってください。


使い方:

(0) {dummies} パッケージをインストールする

(1) 引数 data にダミー変数化したいデータフレームを入れる

(2) 回帰等で多重共線性を避けるために1列除きたい場合は is.drop = TRUE とする

convDummies <- function(data, is.drop = FALSE){
  library(dummies)
  
  N <- ncol(data)
  row_names <- names(data)
  
  names_list <- c()
  new_data <- rep(NA, nrow(data))
  for(n in 1:N){
    unique_value <- sort(unique(data[,n]))
    dummied_data <- dummy(data[,n])
    
    if(is.drop == TRUE){
      new_data <- cbind(new_data, dummied_data[,-ncol(dummied_data)])
      names_list <- c(names_list, 
                      paste(row_names[n], unique_value, sep = ".")[-ncol(dummied_data)])
    } else {
      new_data <- cbind(new_data, dummied_data)
      names_list <- c(names_list, paste(row_names[n], unique_value, sep = "."))
    }
  }
  
  new_data <- as.data.frame(new_data)
  names(new_data) <- c("temp", names_list)
  
  return(new_data[,-1])
}

ベイジアン仮説検定

ベイズ統計学

前に学内の勉強会でベイジアン仮説検定について発表したので、その時のスライドをアップロードします。

ベイズを使って頻度主義の区間推定的なものだけでなく、点推定値の検定を行う方法も紹介しています。


vec 演算子とトレース

線形代数

行列を以下の様に定義する。


A=(\textbf{a}_{1},\ldots,\textbf{a}_{m})
ここで \textbf{a}_{j},\ j=1,\ldots,m は縦ベクトルとする。

このとき


vec(A)=
\begin{bmatrix}
\textbf{a}_{1}\\ 
\vdots\\
\textbf{a}_{m}
\end{bmatrix}


また基本的な操作として



\begin{eqnarray*}
&&tr(A^{T}B)=(vec(A))^{T}vec(B)\\
&&tr(ABC)=(C^{T}\otimes A)vec(B)
\end{eqnarray*}


がある。

{rpart}{partykit}のプロットで日本語を使用する方法

R package

自分用メモ。実行例は後日載せるかも。

fit.dt <- rpart(y ~., data)
plot(as.party(fit.dt), gp = gpar(fontfamily = "Osaka"))

pythonによる粒子フィルタの実装

python ベイズ統計学 時系列解析

2階差分トレンド+季節(週)トレンドを考慮した以下の状態空間モデル(線形ガウス状態空間モデル)

{ \displaystyle
\begin{eqnarray*}
y_{t}&=&\mu_{t}+s_{t}+w_{t},\qquad w_{t}\sim N(0,\sigma^{2})\\
\mu_{t}&=&2\mu_{t-1}-\mu_{t-2}+v_{t},\qquad v_{t}\sim N(0,\alpha_{\mu}^{2}\sigma^{2})\\
s_{t}&=&-\sum_{i=1}^{6}s_{t-i}+z_{t},\qquad z_{t}\sim N(0,\alpha_{s}^{2}\sigma^{2})
\end{eqnarray*}\\
}

について、粒子フィルターをpythonで実装しました。
ここでy_{t},\ \mu_{t}, s_{t}はそれぞれ時刻tにおける観測値、トレンド(平均)、季節トレンドです。
状態空間モデルの詳しい説明については

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)


が分かりやすいです。
行列とベクトルの表記を使って綺麗に状態空間モデルのモデリングや種々のアルゴリズムの導出を行っていますし、P.89には粒子フィルタの疑似コードが載っています。


ただしこの本に従って実際に実装するとなると、モデルによっては膨大な量の行列とベクトルを実装しなければなりません。
特に拡大状態ベクトルを用いて固定ラグ平滑化を導出していますが、いざ固定ラグ平滑化を実装しようとなると拡大状態ベクトルの実装が面倒だと感じます。また固定ラグ平滑化については実装例が載っていません。


そこで、固定ラグ平滑化まで含めて粒子フィルタをpythonで実装しましたので参考にしていただけたらと思います。

github.com


使用したデータは上記書籍のサポートページ
http://daweb.ism.ac.jp/yosoku/
からダウンロードできます。



固定ラグ平滑化の概要は以下であると解釈しています。予測とフィルタリングは普通の粒子フィルタです。縦に5つ並んでいる丸が粒子で、この図では5つの粒子が描かれています。x_{i,t|t-1}t-1時点でのデータで条件付けられたt時点でのi番目の粒子ということを表しています。粒子フィルタはシステムモデルによってt-1時点のデータからt時点の値を"予測"し、実際にt時点でのデータが手に入った時点で観測モデルの尤度でもって状態xをフィルタリング(更新)します。

f:id:kefits:20160805172611p:plain


さて、t時点での最新のデータが手に入った時点で例えばt-1t-2の過去のデータを更新するというのが平滑化という操作で、平滑化の方法はMCMCといった方法で実行することもできますが、粒子フィルターの枠組みでは固定ラグ平滑化という方法が一般に用いられます。その固定ラグ平滑化の概要は上図のように、最新のデータでもって現在の状態xをフィルタリングした後に、そのフィルタリングされた粒子xは前の時点t-1ではどの粒子から派生したものかというのを逆向きに辿ることで平滑化ができる(と私は解釈しています)。

青い四角で囲まれた数字は、前の時点のどの(上から何番目の)粒子から派生したかのインデクスを表しています。上図ではラグ幅が1の固定ラグ平滑化を実行していることになっていますが、プログラムを実装する上では上図のように各時点の粒子がどの粒子から派生したものかのインデックスを記憶しておいて、それをデータが手に入ってフィルタリングをするたびに記憶しておいたインデックスを頼りに過去に辿っていけば平滑化ができるということです。




今回は冒頭で紹介した書籍を参考に、pythonで実際に粒子フィルタを実装をしてみたわけですが、これをStanで実装した例として@berobero11さんのスライドやブログ

www.slideshare.net


statmodeling.hatenablog.com


は大変参考になります。





間違いや質問等ございましたら@kefismまで遠慮なくご連絡ください。