マルチプロセッシングによるMCMC／Python

ＭＣＭＣのマルチプロセッシングPythonプログラミング

Markov chain Monte Carlo (MCMC)の複数のマルコフ鎖をPythonにおいて同時に走らせるためにマルチプロセッシングを使ってみた。使用した統計モデルは１変量正規分布である。すなわち、平均、分散の正規分布の母集団から個のデータ、・・・、が独立に採られたとする。

である。パラメータとして、分散の平方根である標準偏差を使うことにする。このとき、独立な個のデータ、・・・、の得られる確率は

と書くことができる。をとの関数と見たものは尤度関数と呼び、と書く。すなわち、

また、事後確率は

で与えられる。データ数が多くなると、事後確率は主に尤度関数によって決まるようになる。すなわち、データの情報は尤度関数で表される。データ数が多くないときも、事後確率が尤度関数で決まるようにするため、事前確率を定数とおく。すなわち、

ここで、

すなわち、を一様分布とする。ただし、が０でない領域（サポート）にわたる積分が１になるようにするため、サポートは十分に広い有界な領域であるとする。コンピュータでの乱数は有界な範囲での値しかとらないので、その範囲をカバーする十分な大きさのサポートであるとすれば、実質的には何ら制約を受けないことになる。本ウェブサイトのプログラムで用いている乱数は、＜別のウェブサイト＞にあるものに正規乱数生成のメソッドを加えたものである。この乱数は、お互いに独立な乱数系列が容易に得られるものである。

式（１）で与えられる事後確率を、Metropolis-within-GibbsによるMCMC（Markov chain Monte Carlo）によって求める。Metropolis-within-Gibbsについては、岡本安晴「心理学データ分析と測定」勁草書房の付録C「ベイズ的分析法」pp. 227-232で説明している。ベイズ的方法については、岡本安晴「統計学を学ぶための数学入門［下］」2009、培風館、pp．153－157で説明している。

式（１）の場合のMCMCの手順は、以下のようになる。

提案分布

を用意する。提案分布は、ウェブサイト＜ベイズデータ分析＞で説明されている形では、例えばは、とも条件に含まれる形であるが、に関しては定数であるものを選ぶので条件には表示していない。

　およびの値は、MCMCにおける更新確率（acceptance probability）が適切な範囲になるように調節する。具体的には、関数scale_sgmにより、以下のように適応的に調整している。

v, ck = scale_sgm(r_acpt_mean, ck)

params['genMuSgm'] = params['genMuSgm'] * v

v, ck = scale_sgm(r_acpt_sd, ck)

params['genSDSgm'] = params['genSDSgm'] * v

params['genMuSgm']とparams['genSDSgm']がとであり、r_acpt_meanとr_acpt_sdは直前のMCMCにおける更新比率である。関数scale_sgmは以下のように与えられているので、更新比率が0.3～0.5の範囲になるように自動的に調整される。

def scale_sgm( acpt_r, ck ):

v = 1.0

if (acpt_r > 0.5):

ck += 1

v = 1.0 + 2.0 * (acpt_r - 0.4) / 0.6

elif (acpt_r < 0.3):

ck += 1

v = 1.0 / (1.0 + 2.0 * (0.4 - acpt_r) / 0.4)

return v, ck

以上の準備の下で以下のステップを実行する。

ステップ０

とおき、初期値およびを設定する。最初のMCMCにおける初期値は、データの平均値および標準偏差を初期値として設定する。提案分布におけるおよびの値もデータの標準偏差を基に設定する。２回目以降のMCMCにおいては、直前のMCMCでの平均値および調整値を用いる。

ステップ１

の更新候補値を

ととる。

採択確率を次式

により設定する。

確率で、

確率で

とおく。

ステップ２

の更新候補値を

ととる。

採択確率を次式

により設定する。

確率で、

確率で

とおく。

ステップ３

とおく。の値が目標数に達したならば終了するが、達していなときはステップ１に戻る。

ここで、上のアルゴリズムにおける計算上の注意点をいくつか上げておく。

ステップ１における採択確率の計算において

の部分は、

であることから、

である。

確率の比

の計算において、データが多くなるとPythonの実数型（double型）で計算した場合、分子あるいは分母の確率が０になる可能性が高くなる。これは、double型の表しうる０でない数値の範囲が絶対値で以上となっていることによる。すなわち、例えば、個々の確率が0.1のデータを400個ぐらい独立に得たとき、全体の確率はとなり、double型の範囲を超えた小さい値なので０として扱われる。このことを確認するため、次のプログラムを用意した。

v = 1.0

Lv = 0.0

for i in range(0, 325):

v *= 0.1

print("%d v = %0.17e"%((i + 1), v))

上のプログラムの実行結果を図１に示す。

・

図１

を過ぎた行から精度が低下し、を過ぎたあたりで０になっていることがわかる。double型の有効桁数は１５である。

確率の積が０になることを避けるため、対数をとって計算する。すなわち、

である。このとき、

となる。

対数をとって計算した値を戻すときは指数関数を使うが、このときもdouble型の扱える数値の範囲による制約がある。これを調べるプログラムとして次のものを用意した。

import math

v = 700.0

for i in range(0, 10):

try:

print("exp(%f) = %0.15e exp(%f) = %0.15e"%(\

v, math.exp(v), -v, math.exp(-v)))

v += 3.0

except Exception as e:

print(e)

print("v = %f (exp overflow) exp(%f) = %0.15f"%(\

v, -v, math.exp(-v)))

上のプログラムの実行結果が図２である。

図２

のの値が負の方向に絶対値が大きくなると、の値は０になっているが、逆にの値が正の方向に絶対値が大きくなるとoverflowエラーで例外が生成されている。

対数をとって計算した値を採択確率に戻すときは、確率の計算なので、において、の値が１以上になるときは、の値がいくら大きくなっても関数の値は１になる関数を用意して用いることにする。関数Exp0はそのための関数である。

def exp0( x ):

if (x > 0.0):

return 1.0

else:

return math.exp(x)

ステップ２における採択確率の計算もステップ１の場合と同様である。対数をとっての計算は以下のようになる。

ただし、

のときは、

となるので、上の対数による計算は行わず、採択確率を０とおいている。採択の判定は、一様乱数との比較で行っているので、比較の基準値 a の値を０以下の－１とおいて、実質的に採択確率を０としている。

上のMCMCのアルゴリズムを表すコードは次のsubprg.pyで宣言さている関数mcmc_funcとして実装されている。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ subprg.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

import math

from RNGen import RNp2to191 as RN

def sqr( v ):

return v * v

def exp0( x ):

if (x > 0.0):

return 1.0

else:

return math.exp(x)

def scale_sgm( acpt_r, ck ):

v = 1.0

if (acpt_r > 0.5):

ck += 1

v = 1.0 + 2.0 * (acpt_r - 0.4) / 0.6

elif (acpt_r < 0.3):

ck += 1

v = 1.0 / (1.0 + 2.0 * (0.4 - acpt_r) / 0.4)

return v, ck

def mcmc_func( params ):

"""

params = {'jump': 'n': 'x': 'NSimu':

'init_mean': 'init_sd': 'genMuSgm': 'geSDSgm': }

"""

jump = params['jump']

n = params['n']

x = params['x']

NSimu = params['NSimu']

init_mean = params['init_mean']

init_sd = params['init_sd']

genMuSgm = params['genMuSgm']

genSDSgm = params['genSDSgm']

sumX = 0.0

ssumX = 0.0

for i in range(0, n):

sumX += x[i]

ssumX += x[i] ** 2.0

mean = [init_mean]

sd = [init_sd]

acpt_mean = 0

acpt_sd = 0

rn = RN()

# Independent sequences can be obtained for different values of jump

for i in range(0, jump):

rn.jump(100)

# NSimu samples from MCMC

for t in range(0, NSimu):

if (t % 1000) == 0:

print("t/NSimu = %d/%d for jump = %d"%(t, NSimu, jump))

y = rn.normalMS(mean[t], genMuSgm)

log_a = (-2.0 * sumX * (y - mean[t]) + n * (sqr(y) - sqr(mean[t]))) / (-2.0 * sqr(sd[t]))

a = exp0(log_a)

if (rn.uni() < a):

mean.append(y)

acpt_mean += 1

else:

mean.append(mean[t])

y = rn.normalMS( sd[t], genSDSgm )

a = 0.0

if (y <= 0.0):

a = -1.0

else:

log_a = n * math.log(sd[t] / y) \

+ 0.5 * ((1.0 / sqr(sd[t])) - (1.0 / sqr(y))) \

* (ssumX - 2.0 * mean[t + 1] * sumX + n * sqr(mean[t + 1]))

a = exp0(log_a)

if (rn.uni() < a):

sd.append(y)

acpt_sd += 1

else:

sd.append(sd[t])

return mean, acpt_mean, sd, acpt_sd

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

提案分布、の調整は、 (‘genMuSgm’) および（’genSDSgm’）の調整として、以下のコードで自動的に行っている。

ck = -1

n_adap = 0

while True:

print("ck = %d n_adap = %d"%(ck, n_adap))

L_mean, acpt_mean, L_sd, acpt_sd = mcmc_func( params )

sum_mean = 0.0

sum_sd = 0.0

for t in range(1, NSimu + 1):

sum_mean += L_mean[t]

sum_sd += L_sd[t]

params['init_mean'] = sum_mean / NSimu

params['init_sd'] = sum_sd / NSimu

r_acpt_mean = acpt_mean / float(NSimu)

r_acpt_sd = acpt_sd / float(NSimu)

ck = 0

v, ck = scale_sgm(r_acpt_mean, ck)

params['genMuSgm'] = params['genMuSgm'] * v

v, ck = scale_sgm(r_acpt_sd, ck)

params['genSDSgm'] = params['genSDSgm'] * v

if (ck == 0):

print("ck = %d"%ck)

break

n_adap +=1

提案分布と初期値は、適応的に自動化されて設定されるので、気にせずにプログラムを実行して分析を行うことができる。

MCMCによるサンプリングの結果は、ヒストグラム（図３、図４）として表示される。

図３

図４

事後分布の四分位数および９５％CIは、出力ファイルに出力されているので、プログラムの実行終了後、エディタで開いて見ることができる（図５）。

図５

メインプログラムのコードは以下のとおりである。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ mainprg.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

from multiprocessing.pool import Pool

import matplotlib.pyplot as plt

import math

from subprg import scale_sgm, mcmc_func

if __name__ == "__main__":

# Prepare the input data file

s_in = input("Input data file = ")

f_in = open(s_in, "r")

# Set the contents of the input data file in the object data

data_f = f_in.readlines()

f_in.close()

# Prepare the output file

s_out = input("Output data file = ")

f_out = open(s_out, "w")

f_out.write("Input data file = " + s_in + "\n\n")

# Set the data in list x

pos = 0

while True:

if data_f[pos].find("/", 0) == 0: break

pos += 1

n = 0

x = []

while True:

pos += 1

if data_f[pos].find("/", 0) == 0: break

x.append(float(data_f[pos]))

n += 1

print("X...")

i = 1

for v in x:

f_out.write("%5d:"%i + " %f\n"%v)

i+=1

f_out.write("\n")

sum = 0.0

for v in x:

sum += v

mean = sum / n

print("Mean = %f"%mean)

f_out.write("Mean = %f\n"%mean)

ssum = 0.0

for v in x:

ssum += (v - mean) ** 2.0

var = ssum / n

sd = math.sqrt(var)

print("Variance = %f"%var + " sd = %f"%sd)

f_out.write("Variance = %f"%var + " sd = %f\n"%sd)

uvar = ssum / (n - 1.0)

sd_uvar = math.sqrt(uvar)

print("Unbiased Var. = %f"%uvar + " sd(uvar) = %f"%sd_uvar)

f_out.write("Unbiased Var. = %f"%uvar + " sd(uvar) = %f\n"%sd_uvar)

NSimu = 3000

# Initial values of the parameters of MCMC

params = {'jump': 0,

'n': n,

'x': x,

'NSimu': NSimu,

'init_mean': mean,

'init_sd': sd,

'genMuSgm': 0.1 * sd,

'genSDSgm': 0.1 * sd }

# Adjust the parameters of MCMC

ck = -1

n_adap = 0

while True:

print("ck = %d n_adap = %d"%(ck, n_adap))

L_mean, acpt_mean, L_sd, acpt_sd = mcmc_func( params )

sum_mean = 0.0

sum_sd = 0.0

for t in range(1, NSimu + 1):

sum_mean += L_mean[t]

sum_sd += L_sd[t]

params['init_mean'] = sum_mean / NSimu

params['init_sd'] = sum_sd / NSimu

r_acpt_mean = acpt_mean / float(NSimu)

r_acpt_sd = acpt_sd / float(NSimu)

ck = 0

v, ck = scale_sgm(r_acpt_mean, ck)

params['genMuSgm'] = params['genMuSgm'] * v

v, ck = scale_sgm(r_acpt_sd, ck)

params['genSDSgm'] = params['genSDSgm'] * v

if (ck == 0):

print("ck = %d"%ck)

break

n_adap +=1

# The main stage of MCMC

pool = Pool()

NSimu = 10000

params['NSimu'] = NSimu

n_mcmc = 4 # The number of concurrent MCMCs

ary_params = []

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

# Run the processes of MCMCs

print("The main MCMCs started.")

results = pool.map(mcmc_func, ary_params)

print("The main MCMCs ended.")

# Merge the samples from MCMCs

mrg_mean = []

mrg_sd = []

for t_mean, t_acpt_mean, t_sd, t_acpt_sd in results:

print("t_acpt_mean = %d t_acpt_sd = %d"%(t_acpt_mean, t_acpt_sd))

mrg_mean += t_mean[1:NSimu+1]

mrg_sd += t_sd[1:NSimu+1]

# Calculate the quantiles

mrg_mean.sort()

L_mean = mrg_mean[int(NSimu * n_mcmc * 0.025)]

Q1_mean = mrg_mean[int(NSimu * n_mcmc * 0.25)]

Med_mean = mrg_mean[int(NSimu * n_mcmc * 0.5)]

Q3_mean = mrg_mean[int(NSimu * n_mcmc * 0.75)]

U_mean = mrg_mean[int(NSimu * n_mcmc * 0.975)]

print("\nQ1_mean = %f Med_mean = %f Q3_mean = %f"%(Q1_mean, \

Med_mean, Q3_mean))

print("95%%CI = [%f, %f]"%(L_mean, U_mean))

f_out.write("\nQ1_mean = %f Med_mean = %f Q3_mean = %f\n"%(Q1_mean, \

Med_mean, Q3_mean))

f_out.write("95%%CI = [%f, %f]\n"%(L_mean, U_mean))

plt.figure(figsize = (12, 9))

plt.title("Posterior Distribution of the Mean")

plt.hist(mrg_mean, bins = 50)

plt.show()

mrg_sd.sort()

L_sd = mrg_sd[int(NSimu * n_mcmc * 0.025)]

Q1_sd = mrg_sd[int(NSimu * n_mcmc * 0.25)]

Med_sd = mrg_sd[int(NSimu * n_mcmc * 0.5)]

Q3_sd = mrg_sd[int(NSimu * n_mcmc * 0.75)]

U_sd = mrg_sd[int(NSimu * n_mcmc * 0.975)]

print("\nQ1_sd = %f Med_sdn = %f Q3_sd = %f"%(Q1_sd, \

Med_sd, Q3_sd))

print("95%%CI = [%f, %f]"%(L_sd, U_sd))

f_out.write("\nQ1_sd = %f Med_sd = %f Q3_sd = %f\n"%(Q1_sd, \

Med_sd, Q3_sd))

f_out.write("95%%CI = [%f, %f]\n"%(L_sd, U_sd))

plt.figure(figsize = (12, 9))

plt.title("Posterior Distribution of the SD")

plt.hist(mrg_sd, bins = 50)

plt.show()

f_out.close()

print(s_out + " was saved.")

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

主MCMCにおいて、以下のようにPoolを用いてマルチプロセスによる計算を行っている。Poolによるマルチプロセスの説明は、ウェブサイト＜Pythonマルチプロセシング＞で行っている。

pool = Pool()

NSimu = 10000

params['NSimu'] = NSimu

n_mcmc = 4 # The number of concurrent MCMCs

ary_params = []

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

# Run the processes of MCMCs

results = pool.map(mcmc_func, ary_params)

各プロセスにおけるＭＣＭＣのマルコフ鎖は’jump’の値が以下のように

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

とお互いに異なるように設定されている。これにより、subprg.pyで宣言されている関数mcmc_funcにおいて

for i in range(0, jump):

rn.jump(100)

と乱数のメソッドjump(100)の実行回数が異なることになり、各プロセスのマルコフ鎖は独立なものとなる。メソッドjump(100)については、＜別のウェブサイト＞で説明している。

プログラムを実行すると、入力データファイル名の設定が求められる（図６）。

図６

入力データファイルは、図７のように、スラッシュ / で始まる行でデータを挟む形で用意する。データは、１行に１つ書く。

図７

入力データファイル名を設定して、Enterキーを押すと、次に出力用ファイル名を聞いてくる（図８）。

図８

出力ファイル名は、適当なテキストファイル名を設定する。出力ファイル名を設定してEnterキーを押すと計算が始まる。MCMCの計算が終了すると、まず、平均の事後分布のヒストグラムが表示される（図３）。平均の事後分布のヒストグラムのウィンドウの右上角のX印をクリックして閉じると、次にの事後分布のヒストグラムが表示される（図４）。このヒストグラムのウィンドウを閉じるとプログラムの実行終了となる。四分位数などは、出力用ファイルに書き出されているので、エディタで開いて見ることができる。

子プロセスでの関数printによる出力は、Shellのウィンドウには表示されないが、端末（コンソール）で実行すると、端末に表示される。図９に実行例を示す。

図９

プログラムなどは圧縮ファイルmcmcMultiProcPrg.zipとしてまとめた。ファイル名mcmcMultiProcPrg.zipをクリックしてダウンロードすることができる。