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()]
)