Up

PyMC3‚ΜŠΦ”sample‚Μtarget_acceptˆψ”‚ΜŽw’θ–@•ΟX‚Ι‚Β‚’‚āi2021.09j

 

PyMC3‚ΜƒTƒ“ƒvƒŠƒ“ƒOŠΦ”sample‚Μˆψ”nuts_kwargs‚ͺ–³‚­‚Θ‚Α‚Α‚Δ‚’‚½i2021.09jB‚±‚Μ‚½‚߁AΩ’˜ƒ‰ͺ–{ˆΐ°u‚’‚ά‚³‚η•·‚―‚Θ‚’Python‚Εƒf[ƒ^•ͺΝvŠΫ‘Po”Ł„‚Μp. 232‚ΜŠΦ”sample‚π—p‚’‚Δ‚’‚ιƒXƒNƒŠƒvƒg‚ΝˆΘ‰Ί‚Μ‚ζ‚€‚ɕύX‚·‚ι•K—v‚ͺ‚ ‚ιB

 

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

                         target_accept = 0.95)

                        #nuts_kwargs = dict(target_accept = 0.95))

 

‚·‚Θ‚ν‚ΏAtarget_accept‚Μ’l‚πŽ«‘Œ^‚Ι‚¨‚―‚ι’l‚Ζ‚΅‚Δ‚Ε‚Ν‚Θ‚­A’ΌΪŽw’θ‚·‚κ‚Ξ‚ζ‚’B

ƒvƒƒOƒ‰ƒ€‘S‘Μ‚ΜƒXƒNƒŠƒvƒg‚πAƒŠƒXƒg‚P‚Ζ‚΅‚ΔŽ¦‚·Bƒtƒ@ƒCƒ‹‚Νˆ³kƒtƒ@ƒCƒ‹PoissonRegressionPyMC3.zip‚Ζ‚΅‚Δ‚ά‚Ζ‚ί‚½B

‚Θ‚¨APyMC3‚ΜƒvƒƒOƒ‰ƒ€‚πWindows‚Ε—˜—p‚·‚ι‚Ζ‚«‚́AJupyter Notebookγ‚ΕŽΐs‚Ε‚«‚ιB‚±‚ΜƒEƒFƒuƒTƒCƒg‚πŽQl‚³‚κ‚½‚’B

 

ƒŠƒXƒg‚P

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