Up

PyMC3の関数sampleのtarget_accept引数の指定法変更について（2021.09）

PyMC3のサンプリング関数sampleの引数nuts_kwargsが無くなっっていた（2021.09）。このため、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞のp. 232の関数sampleを用いているスクリプトは以下のように変更する必要がある。

trace = pm.sample(draws = 1000, tune = 2000,

target_accept = 0.95)

#nuts_kwargs = dict(target_accept = 0.95))

すなわち、target_acceptの値を辞書型における値としてではなく、直接指定すればよい。

プログラム全体のスクリプトを、リスト１として示す。ファイルは圧縮ファイルPoissonRegressionPyMC3.zipとしてまとめた。

なお、PyMC3のプログラムをWindowsで利用するときは、Jupyter Notebook上で実行できる。このウェブサイトを参考されたい。

リスト１

import matplotlib.pyplot as plt

import numpy as np

import pymc3 as pm

RawData = [[1, 1, 0, 0, 1, 0, 0],

[1, 1, 0, 0, 0, 1, 0],

[10, 1, 0, 0, 0, 0, 1],

[3, 0, 1, 0, 1, 0, 0],

[2, 0, 1, 0, 0, 1, 0],

[10, 0, 1, 0, 0, 0, 1],

[6, 0, 0, 1, 1, 0, 0],

[2, 0, 0, 1, 0, 1, 0],

[4, 0, 0, 1, 0, 0, 1]]

F = []

Xr1 = []

Xr2 = []

Xr3 = []

Xc1 = []

Xc2 = []

Xc3 = []

for v in RawData:

F.append(v[0])

Xr1.append(v[1])

Xr2.append(v[2])

Xr3.append(v[3])

Xc1.append(v[4])

Xc2.append(v[5])

Xc3.append(v[6])

for i in range(9):

print(' ', F[i], ' ', Xr1[i], ' ', Xr2[i], ' ', Xr3[i],

' ', Xc1[i], ' ', Xc2[i], ' ', Xc3[i])

with pm.Model() as PoiModel:

mu = pm.Normal('mu', mu = 0.0, sd = 100.0)

a2 = pm.Normal('a2', mu = 0.0, sd = 100.0)

a3 = pm.Normal('a3', mu = 0.0, sd = 100.0)

b2 = pm.Normal('b2', mu = 0.0, sd = 100.0)

b3 = pm.Normal('b3', mu = 0.0, sd = 100.0)

g22 = pm.Normal('g22', mu = 0.0, sd = 100.0)

g23 = pm.Normal('g23', mu = 0.0, sd = 100.0)

g32 = pm.Normal('g32', mu = 0.0, sd = 100.0)

g33 = pm.Normal('g33', mu = 0.0, sd = 100.0)

LogLmbd = mu + a2 * Xr2 + a3 * Xr3 + b2 * Xc2 + b3 * Xc3 \

+ g22 * Xr2 * Xc2 + g23 * Xr2 * Xc3 \

+ g32 * Xr3 * Xc2 + g33 * Xr3 * Xc3

Lmbd = np.exp(LogLmbd)

F_obs = pm.Poisson('F_obs', mu = Lmbd, observed = F)

trace = pm.sample(draws = 1000, tune = 2000,

target_accept = 0.95)

#nuts_kwargs = dict(target_accept = 0.95))

pm.traceplot(trace)

plt.show()

summary = pm.summary(trace)

print(summary)

f = open('summary.txt', 'w')

f.write('Summary...\n{}'.format(summary))

smpl_a3 = trace['a3']

a3_p2p5 = np.percentile(smpl_a3, 2.5)

a3_p97p5 = np.percentile(smpl_a3, 97.5)

print('\n95%CI for a3 = [{0:.5f}, {1:.5f}]'.format(a3_p2p5, a3_p97p5))

f.write('\n\n95%CI for a3 = [{0:.5f}, {1:.5f}]\n'.format(a3_p2p5, a3_p97p5))

cnt_a3_pos = (smpl_a3 > 0.0).sum()

print('\nP(a3 > 0.0) = {0:.3f}'.format(cnt_a3_pos / len(smpl_a3)))

f.write('\nP(a3 > 0.0) = {0:.3f}\n'.format(cnt_a3_pos / len(smpl_a3)))

smpl_b3 = trace['b3']

cnt_b3_pos = (smpl_b3 > 0.0).sum()

print('P(b3 > 0.0) = {0:.3f}'.format(cnt_b3_pos / len(smpl_b3)))

f.write('P(b3 > 0.0) = {0:.3f}\n'.format(cnt_b3_pos / len(smpl_b3)))

smpl_g33 = trace['g33']

cnt_g33_neg = (smpl_g33 < 0.0).sum()

print('P(g33 < 0.0) = {0:.3f}'.format(cnt_g33_neg / len(smpl_g33)))

f.write('P(g33 < 0.0) = {0:.3f}\n'.format(cnt_g33_neg / len(smpl_g33)))

f.close()

Up