cythonを使ってMCMCを実装 | サンプリング編

この記事では,cythonを用いたMCMCの実装を行う.
なぜcythonを使うかと言われれば,ただでさえ時間のかかるMCMCをpythonで実装したところで,収束まで時間が掛かりすぎてしまう.そこで,cythonを使ってMCMCを実装しようという流れになった.

cythonのコードを外部からimportする方法 を先に読むとcython独特の記法が把握出来る. 
コードは,git clone https://github.com/akitoshiblog/blogsiteでmcmc_cythonの下にある.

*2020年1月19日に大幅に書き換えました.

使い方

メインのコードは.myMCMC.pyx/myMCMC.pxdとmyStats.pxd/myStast.pyxである.図示に関しては cythonを使ってMCMCを実装 | 図示編 を参照のこと. 

動作はjupyter notebook上を想定するので, cythonのセルマジックを使えるように,%load_ext Cython を回す.

次に,myMCMCからクラスを引き継いで次の三点を行う.
1. super().__init__ でmcmcのアウトプットの保存用のパスと変数名,ランダムウォーク用の標準偏差の値を渡す.
2. パラメータを受け渡しするための,setParametersとsetInitialを定義する.
3. 評価関数用のlogLikelihoodを定義する. 

Rosenbrock’s banana functionを尤度関数として用いてパラメータを a=1,b=100の下で,(x,y)=(1,1)という大域最小解が得られるか確かめる.

モデルのためのcythonのコードは以下のようにする.MCMCは関数を最大にしようとするので,返り値は負の値にしておく.

jupyter notebook

%%cython -a

cimport cython
# custom function
cimport myStats as ms
cimport myMCMC 

cdef class model(myMCMC.MCMC):
    ''' Log-Likelihood function is defined.
    '''
    cdef public double x
    cdef public double y
    
    def __init__(self,path,parName,sdRandomWalk):
        super().__init__(path,parName,sdRandomWalk)
        
    cpdef int setParameters(self):
        self.x = self.par[0]
        self.y = self.par[1]
        
    cpdef int setInitial(self): 
        self.par[0] = self.x        
        self.par[1] = self.y        
        
    ########## Data setting part ##################
    cdef public double a
    cdef public double b
    def setInitParams(self):
        self.a = 1
        self.b = 100
        self.x = 2
        self.y = 2 
        
    ############# calculate log likelihood ###############
    @cython.cdivision(True)
    cpdef double logLikelihood(self):
        cdef double v 
        cdef double lv
        v = ms.pow(self.a - self.x,2) + self.b*ms.pow(self.y  - self.x*self.x,2) + 1e-5
        lv = ms.log(v)
        return(- lv)

大事なのは,logLikelihoodの関数がきっちりc言語としてコンパイルされること.それ以外の部分はMCMCの最中には呼ばれないので,大量のパラメーター推定の際には,名前や初期値はforやlist,dictを用いて設定することもできる.
modelを定義したら,必要なパラメータを与えてrunさせる.

parName = ["x","y"]
sdRandomWalk = [0.5,0.5]
mcmc = model(path ="output/result",
             parName = parName,
             sdRandomWalk = sdRandomWalk)
# set initial values for each parameter
mcmc.setInitParams()

mcmc.nSim = 1000000
mcmc.nThin = 1000
print(mcmc.path2Save)
mcmc.run()

これでMCMCのアウトプットが得られる.

ソースコードの準備の仕方

myMCMC.pyx/myMCMC.pxdにmcmcのアーキテクチャーが,myStats.pxd/myStast.pyxに基本的なc言語の関数と正規分布に関する定義を記してある.

一度,これらのファイルをいじったら再度cythonを用いてコンパイルする必要がある.その場合は,python setup.py build_ext --inplace を用いて実行すること.
setup.py のコードは以下の通りである.

from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
import numpy as np

# if you run this code, type tye following code
# >>> python setup.py build_ext --inplace

ext_modules = [ 
            Extension("myStats", ["myStats.pyx"]),
            Extension("myMCMC", ["myMCMC.pyx"]),
            ]

setup( name = "Myproject", 
        ext_modules =  cythonize(ext_modules,annotate=True),
        include_dirs=[np.get_include()]
)

python
スポンサーリンク
あきとしのスクラップノート

コメント

  1. […] この記事は cythonを使ってMCMCを実装 | サンプリング編 の続きです.前回の記事でoutputされたデータを図示する. […]

タイトルとURLをコピーしました