一日一万字の感謝の写経

一日一万字の感謝の写経をして強くなります。そのうち本を置き去りにして何かを作り始める予定です。

Pythonで体験するベイズ推論 PyMCによるMCMC入門 第3章 (3.2全部) 写経

PyMC Python ベイズ理論 写経

目次

始めに

　この記事では、Pythonで体験するベイズ推論 PyMCによるMCMC入門

の第3章の3.2を写経した記録をまとめます。前回はこちらPythonで体験するベイズ推論 PyMCによるMCMC入門 第3章 (3.1全部) 写経。 実行はJupyter Labにて実行しています。

　写経は以下のようにまとめました。

  • 写経ではありますが、関数を作成し、コードの再実行をやりやすいようにした部分があります。

  • より良いと思われるコードを考えた場合は書き換えてコメントを添えるようにし、変更点をなるべく明示するようにしてあります。

  • 個人的に気になった点のメモを残すようにしてあります。同じような疑問を持った方の助けになれば幸いです。

  • 以前書いたコードと同じようなコード（例えばグラフの描写等）は効率化のために飛ばしているところもあります。

　Pythonで体験するベイズ推論 PyMCによるMCMC入門を読んでいている際のちょっとした正誤表代わりになればと思います。

この記事で使用する主なモジュール、設定

　この記事では、主に以下のモジュールや設定を使用しています。

from IPython.core.pylabtools import figsize
from matplotlib import pyplot as plt
import numpy as np
import pymc as pm
from pymc.Matplot import plot as mcplot

3 MCMCの中をのぞいてみよう

3.2 自己相関

figsize(25, 8)
plt.rcParams['font.size'] = 17

x_t = pm.rnormal(0, 1, 200)
x_t[0] = 0

y_t = np.zeros(200)
for i in range(1, 200):
    y_t[i] = pm.rnormal(y_t[i - 1], 1)

plt.plot(y_t, label="$y_t$", lw=3)
plt.plot(x_t, label="$x_t$", lw=3)

plt.xlabel("Time, $t$")
plt.ylabel("Value")
plt.title("Two different series of random values")
plt.legend()

Output:

    <matplotlib.legend.Legend at 0x27d4f39ecc0>

f:id:koheitsutsumi223:20190331232651p:plain

def autocorr(x):
    result = np.correlate(x, x, mode='full')
    result = result / np.max(result)
    # 通常の割り算だと、結果が浮動小数点になり、arrayのインデックスの引数に渡すと警告が出る。
    #なので、//により商を求めるようにしておいた。
    return result[result.size // 2:]

colors = ["#348ABD", "#A60628", "#7A68A6"]
x = np.arange(1, 200)

plt.bar(x, autocorr(y_t)[1:], width=1, label="$y_t$",
        edgecolor=colors[0], color=colors[0])
plt.bar(x, autocorr(x_t)[1:], width=1, label="$x_t$",
        edgecolor=colors[1], color=colors[1])

plt.legend(title="autocorrelation")
plt.ylabel("Measured correlation \n"
           "between $y_n$ and $y_{t-k}$.")
plt.xlabel("$k$ (lag)")
plt.title("Autocorrelation plot of $y_t$ and $x_t$ "
          "for differing $k$ lags")

Output:

    Text(0.5, 1.0, 'Autocorrelation plot of $y_t$ and $x_t$ for differing $k$ lags')

f:id:koheitsutsumi223:20190331232658p:plain

3.2.2 間引き処理

max_x = 200 // 3 + 1
x = np.arange(1, max_x)

plt.bar(x, autocorr(y_t)[1:max_x], edgecolor=colors[0],
        label="no thinning",
        color=colors[0], width=1)
plt.bar(x, autocorr(y_t[::2])[1:max_x], edgecolor=colors[1],
        label="keeping every 2nd sample",
       color=colors[1], width=1)
plt.bar(x, autocorr(y_t[::3])[1:max_x], edgecolor=colors[2],
        label="keeping every 3rd sample",
       color=colors[2], width=1)

plt.autoscale(tight=True)
plt.legend(title="Autocorrelation plot for $y_t$",
           loc="lower left")
plt.ylabel("Measured correlation \n"
           "between $y_t$ and $y_{t-k}$.")
plt.xlabel("$k$ (lag)")
plt.title("Autocorrelation of $y_t$ (no thinning versus thinning) "
          "at differing $k$ lags")

Output:

    Text(0.5, 1.0, 'Autocorrelation of $y_t$ (no thinning versus thinning) at differing $k$ lags')

f:id:koheitsutsumi223:20190331232705p:plain

p = pm.Uniform("p", 0., 1.)
assignment = pm.Categorical("assignment", [p, 1-p], 
                            size=data.shape[0])

taus = 1.0 / pm.Uniform("stds", lower=0, upper=100, size=2) ** 2 # テキストを読んだ感じ、ここはupper=100なんじゃないんかな？
centers = pm.Normal("centers", [120, 190], [0.01, 0.01], size=2)

@pm.deterministic
def center_i(assignment=assignment, centers=centers):
    return centers[assignment]

@pm.deterministic
def tau_i(assignment=assignment, taus=taus):
    return taus[assignment]

observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True) # これ、どこで使うんだ？
model = pm.Model([p, assignment, taus, centers])

mcmc = pm.MCMC(model)

mcmc.sample(25000, 0, 10)

Output:

     [-----------------100%-----------------] 25000 of 25000 complete in 9.7 sec
figsize(25, 8)
plt.rcParams['font.size'] = 10
mcplot(mcmc.trace("centers", 2), common_scale=False)

Output:

    Plotting centers_0
    Plotting centers_1

f:id:koheitsutsumi223:20190331232710p:plain

　今回の写経は以上です。 ここまで読んでいただきありがとうございました。