💨

MCMC法のギブスサンプリングをやってみよう

2023/04/01に公開

Python

本記事の目的

言語化して浅い理解を深くしていく。今回はベイズ推定で用いるサンプリング方法について。ベイズ推定は、ベイズ統計学という言葉もあるように統計の中でもメジャーだと思われる。まず、このベイズ推定をいつも通り理論→実践で理解する。理論についてはこちらの書籍を参考にさせていただいた。（amazonに飛ぶのでご注意ください）

ベイズ理論とは

ベイズ推定の話に移る前に、その根底にあるベイズ理論について簡単に説明する。これは入力したデータを加味し、次のデータを予測する理論である。身近な一例を挙げる。ある試験を受けようとする学生がいる。試験の結果は合格か不合格の2択なので、試験を受ける前の学生の合格確率はとりあえず1/2と置くことできる。そして、試験後になると学生の手ごたえはよく感じていた。すると、学生の合格確率は3/4ぐらいに上がるだろう。このようにデータをもとに次のデータを予測するのがベイズ理論である。

ベイズ推定とは

先述のベイズ理論と絡めると、直前データをもとに次データの確率を推定する手法である。途中式は省くがベイズの定理を用いれば、次のように表せる。

\pi(\theta|D)[事後確率]=k[定数]*f(D|\theta)[尤度]*\pi(\theta)[事前確率]

ここで、Dはデータ、 $\theta$ は求める母数を表している。また、尤度は「母数がある値を取るとき、観測したデータが起きる確率」である。確率とは似て異なるが、とりあえずデータの起きる確率と思っておけばいい。この事後確率をベイズ推定では求めている。コイン投げの例題を用意したので解いてみましょう。

コインの表が出る確率をθとし、4回投げ表→裏→表→裏の順に出た。このときの最も取りうる確率が大きくなるθを推測せよ。

解説：直感的にθ=2/4=0.5と思う方もいると思うが、その通りである。ただし、もう少し深く理解していこう。θ=0.1であれ、θ=0.9であれ、起きにくいかもしれないが「4回投げ表→裏→表→裏となる」確率は0ではない。今回求めるのは0<θ<1で「4回投げ表→裏→表→裏」確率が最大のθである。
まず、1回目のデータ（表）を読み込む。1回目のデータは事前確率を持たないので、事前確率は0と1の間なら何でも良い。今回は1にして計算する。1回目の事後確率 $\pi_1(\theta|D_1)$ は次のように表せる。

\pi_1(\theta|D_1)=k*\theta*1

θを0<θ<1で積分すると総和は1となる（規格化）ので、

\int_{0}^{1}k*\theta*d\theta=1

よってk=2となるので1回目の事後確率 $\pi_1(\theta|D_1)$ は $2\theta$ 。
次に2回目のデータ（裏）を読み込む。1回目の事後確率が2回目の事前確率となるので

\pi_2(\theta|D_2)=k*(1-\theta)*2\theta

以降は上記を繰り返すので、途中式は省略。各事後確率は下記の通り。
k=3となるので2回目の事後確率 $\pi_2(\theta|D_2)$ は $6*\theta*(1-\theta)$ 。
k=2となるので3回目の事後確率 $\pi_3(\theta|D_3)$ は $12*\theta^2*(1-\theta)$ 。
k=5/2となるので4回目の事後確率 $\pi_4(\theta|D_4)$ は $30*\theta^2*(1-\theta)^2$ 。

この4回目の事後確率 $\pi_4(\theta|D_4)$ が最大値（最も起きやすい）を取る $\theta$ が最も推測値 $\theta$ として適している。途中式は省くが、 $\theta=0.5$ の時、4回目の事後確率が最大値を取るので、推測値は $\theta=0.5$ 。

このように、ある事後確率が次サンプルの事前確率となり、次サンプルの事後確率を求めていく。この事後確率が分かると、次の式を用いることで、求める母数の期待値を求めることができる。

\int \pi(\theta|D)*\theta d\theta

この $\theta$ の期待値のことをEAP推定量という。

区間推定と点推定について

統計を勉強している方は既知のことだと思うが、区間推定では信頼区間を用いて母数を推定している。例えば、95%信頼区間で $\mu$ を求めるとき、 $\mu=\bar{x} \pm 1.96\hat{\sigma} / \sqrt{n}$ のような範囲で求めていると思う。ベイズ推定では確率を使うので最も確率が高いところが母数となる。つまり点推定の1種である。上記の例題でも推定値は範囲でなく、一個であることに是非注目してほしい。

最尤推定とベイズ推定の違い

似た点推定に最尤推定という手法がある。最尤推定とは、その名の通り「最も尤もらしい値を推定」する手法。ここでもコイン投げを例に取る。

コインの表が出る確率をpとし、4回投げ表となった回数は2回だった。この事象が最も起きやすい確率pを求めよ。

解説：「4回投げ2回、表となる」確率をL(p)とすると

L(p)=6*p^2(1-p)^2

これを解いていく。指数の場合は両辺を対数でとってやると計算しやすくなる。

logL(p)=2log(p)+2log(1-p)+log6

両辺をpで微分すると

(logL(p))'=2/p-2/(1-p)=(2-4p)/p(1-p)

よって、logL(p)の極値はp=0.5である。計算は省くが0<p<1で(logL(p))''<0となるため、p=0.5でlogL(p)は最大値を取る。また、logL(p)は単調増加なので、L(p)もp=0.5で最大値を取る。よって最も「4回投げ2回、表となる」確率(L(p))が大きくなるのはp=0.5。

以上から分かるようにベイズ推定は事前情報を使うのに対し、最尤推定は事前情報を使わない。

簡単にベイズ推定を行うには

上記のベイズ推定の例題は、たまたま簡単にベイズ推定できたが、本来ならばこうはいかない。EAP推定量 $\int \pi(\theta|D)*\theta d\theta$ の式を思い出してほしいが、この事後確率 $\pi(\theta|D)$ が全てのサンプルで決まった型の分布でないと計算が難しそうなのは想像つくだろう。また、求める母数が増えてくるとさらに計算が難しくなり、コンピューターでも計算するのが難しくなる。なので、ベイズ推定を簡単に行うには、事前分布と事後分布が同じ形の分布を取らなければならない。そのためには、事前分布と尤度の組み合わせは下記の通りである必要がある。

事前分布	尤度	事後分布
ベータ分布	二項分布	ベータ分布
正規分布	正規分布	正規分布
逆ガンマ分布	正規分布	逆ガンマ分布
ガンマ分布	ポアソン分布	ガンマ分布

※このような事前分布を尤度の自然な共役分布という。

上記の例題では事前分布がベータ分布、尤度が二項分布であったため、事後分布がベータ分布となり、繰り返し処理を行うだけで良かった。

複雑でもベイズ推定を行うには

現実世界では、上記のようなケースは稀である。コンピューターでできる範囲で計算が複雑になってもよいから、ベイズ推定を行うには後述するMCMC法（マルコフ連鎖モンテカルロ法）という乱数発生アルゴリズムを用いる。この手法は、マルコフチェーンとモンテカルロ法を組み合わせた手法であり、サンプリング手法の1種である。この手法を用いることで、 $\pi(\theta|D)$ で重みづけした疑似母数 $\theta_p$ を何個でも生み出すことができる。これの良いところは、サンプルを多く生み出すのが難しいケース（化学実験など）のとき、疑似的に実験回数を嵩増しできるということ。また、EAP推定量を求めるときには、 $\int \pi(\theta|D)*\theta d\theta$ を用いなくても、 $\theta_p$ の平均を取るだけで良くなるのも良いところだ。

マルコフチェーンとは

事後データは、前データの中でも直前のデータのみに依存するモデルを指す。数学的に書くなら、時刻t+1のデータ $X_{t+1}$ は $X_t$ のみに依存し、 $X_{t-1}$ は関係ない。ランダムウォークともいう。マルコフチェーンについては下記リンクで簡単にまとめている。

モンテカルロ法とは

確率で重みづけした乱数を用いて、データサンプリングを行う。モンテカルロ法については下記リンクで簡単にまとめている。

MCMC法の代表的な手法：ギブスサンプリング

MCMC法にはギブスサンプリングやモンテカルロ法など、さらに細分化されている。本記事ではギブスサンプリングに焦点を当てて説明する。他の手法についても、いずれはまとめたい。

書くにあたって、こちらのサイトを参考にさせていただいた。

多次元のサンプルA $(x_1,x_2,x_3,･･･)$ が1つあるとする。このサンプルの次をBとするとき、そのサンプルの座標は $(x_1',x_2,x_3,･･･)$ とする手法をギブスサンプリングという。つまり、1回のサンプリングで変数は1つしか変わらない。このようなサンプルBの確率を条件付き確率（事後確率）といい、 $p(x_1'|x_2,x_3,･･･)$ と表す。この手法の良いところは、次元数を削減できるところである。今回の実装では、条件付き確率で重みづけした乱数を簡単に取得できるよう、条件付き確率が必ずある形の分布になるようにする。その一例が多次元正規分布である。

D次元正規分布の式

N(x|\mu,\sigma)=\frac{1}{\sqrt{(2\pi)^D|\Sigma|}}exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))

2次元正規分布でμとΣを自分で定義し、条件付き確率が正規分布になるか確かめよう。

解説：各変数を次のように定義する。

\mu=0 \\\\ \Sigma^{-1}= \begin{pmatrix} 1 & -a \\ -a & 1 \end{pmatrix} (-1<a<1)

記載はしないが、逆行列の公式

A= \begin{pmatrix} a & b \\ c & d \end{pmatrix}の時、 A^{-1}=\frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix}

が成り立つことを用いて、 $\Sigma$ も求めておこう。
時刻tの時、サンプルが $(x_1,x_2)$ の状態にあるとし、同時分布確率 $p(x_1,x_2)$ を求める。

\begin{align*} p(x_1,x_2)&=\frac{\sqrt{1-a^2}}{2\pi} exp( -\frac{1}{2} (x_1,x_2) \begin{pmatrix} 1 & -a \\ -a & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} ) \\ &=\frac{\sqrt{1-a^2}}{2\pi} exp( -\frac{1}{2} (x_1-ax_2,-ax_1+x_2) \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} ) \\ &=\frac{\sqrt{1-a^2}}{2\pi} exp( -\frac{1}{2} (x_1^2-2ax_1x_2+x_2^2)) \\ &=\frac{\sqrt{1-a^2}}{2\pi} exp(-\frac{1}{2} (1-a^2) x_2^2) exp(-\frac{1}{2} (x_1-ax_2)^2) \end{align*}

ここで、条件付き確率 $p(x_1|x_2)$ は次のように表せる。

p(x_1|x_2)=\frac{p(x_1,x_2)}{p(x_2)}=\frac{p(x_1,x_2)}{\int p(x_1,x_2) dx_1}-(☆)

分母の $\int p(x_1,x_2) dx_1$ からも分かる通り、 $x_1$ 以外の変数は動かないことに注目。次に、分母の $\int p(x_1,x_2) dx_1$ に、求めた $p(x_1,x_2)$ を代入する。

\int p(x_1,x_2) dx_1=\int \frac{\sqrt{1-a^2}}{2\pi} exp(-\frac{1}{2} (1-a^2) x_2^2) exp(-\frac{1}{2} (x_1-ax_2)^2)dx_1

ここで $t=x_1-ax_2$ とすると

\int p(x_1,x_2) dx_1=\frac{\sqrt{1-a^2}}{2\pi} exp(-\frac{1}{2} (1-a^2) x_2^2) \int exp(-\frac{1}{2} t^2)dt

ガウス積分より $\int exp(-ax^2)dx=\sqrt{\frac{\pi}{a}}$ なので、

\int p(x_1,x_2) dx_1=\sqrt{\frac{1-a^2}{2\pi}} exp(-\frac{1}{2} (1-a^2) x_2^2)

以上より☆式は

p(x_1|x_2)=\frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} (x_1-ax_2)^2)

これは、 $\mu=ax_2$ 、 $\sigma=1$ の一次元正規分布であることが分かる。（終）

バーンインについて

MCMC法では、求めるパラメータの初期値は、そのパラメータの取りうる値ならばどんな値でも良いという特徴がある。とてもありがたい特徴なのだが、あまりにも母数から違いすぎると、求めるデータに影響を与える場合がある。そこで、MCMC法の最初の試行回数の何個かは省きましょうというのがバーンインという。下にサンプルの平均値(mu)10、標準偏差(sigma)8の正規分布データからギブスサンプリングした結果をバーンイン無し・有りでの推定結果を記す。

最初のデータだけ値が大きく違う

左：バーンイン無し、右：最初の100データをバーンイン

バーンインの数はデータによりけりなので、色々な値で試してみる必要がある。これについても後ほど実装してみる。

理論：ギブスサンプリングをやってみよう

では、ギブスサンプリングをやってみよう。といっても上記の例題と同様にパラメータの条件付き確率がどのような分布を取るか確認し、それをプログラミングするだけ。今回は既に持っているサンプル $y_i(i=1,･･･,N)$ が一次元ガウス分布（正規分布）である $N(\mu,\sigma)$ に従っているとする。このパラメータ $\mu,\sigma$ をギブスサンプリングしてみよう。まず、各データ、尤度、事前分布は次の通りである。
<既に持っているサンプルの分布>

y_i \sim N(\mu,\sigma^2)(i=1,･･･,N)

<既に持っているサンプルの確率分布>

p(y_1,･･･,y_N|\mu,\sigma)=\Pi_{i=1}^{N}N(y_i|\mu,\sigma^2)

<初期値>

\mu \sim uniform(-\infty,\infty) \\ \sigma \sim uniform(0,\infty)

最初に $\mu$ の条件付き確率を求めてみよう。

p(\mu|\sigma,y)=\frac{p(\mu,\sigma,y)}{p(\sigma,y)}

ここで分子 $p(\mu,\sigma,y)$ は $\mu,\sigma$ がそれぞれ独立していることより、 $p(\sigma|\mu)=p(\sigma)$ なので

p(\mu,\sigma,y)=p(y|\mu,\sigma)*p(\sigma|\mu)*p(\mu)=p(y|\mu,\sigma)*p(\sigma)*p(\mu)

よって、 $p(\mu|\sigma,y)\propto p(y|\mu,\sigma)$ となる。以降は $\mu$ の条件付き確率の分布の形を見たいだけなので、 $\mu$ が入っていない項はどんどん省いていく。（先にも述べたが、条件付き確率では対象パラメータ以外の項は定数とみなしている。）

\begin{align*} p(\mu|\sigma,y) &\propto p(y|\mu,\sigma) \\ &=\Pi_{i=1}^{N}N(y_i|\mu,\sigma^2) \\ &=\Pi_{i=1}^{N}(\frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y_i-\mu)^2}{2\sigma^2})) \\ &\propto \Pi_{i=1}^{N} exp(-\frac{\mu^2}{2\sigma^2}+\frac{y_i\mu}{\sigma^2}) \\ &=exp(-\frac{N}{2\sigma^2}\mu^2+\frac{\Sigma_{i=1}^{N}y_i}{\sigma^2}\mu) \\ &=exp(-\frac{(\mu-\frac{\Sigma_{i=1}^{N}y_i} {N})^2-\frac{(\Sigma_{i=1}^{N}y_i)^2} {N^2}} {\frac{2\sigma^2} {N}}) \\ &\propto \frac{1} {\sqrt{\frac{2\pi\sigma^2} {N}}} exp(-\frac{(\mu-\frac{\Sigma_{i=1}^{N}y_i} {N})^2} {\frac{2\sigma^2} {N}}) \end{align*}

最後の式ではガウス分布の式に合わせるために無理やり定数 $\frac{1} {\sqrt{\frac{2\pi\sigma^2} {N}}}$ をかけた。これにより、 $p(\mu|\sigma,y)\sim N(\frac{\Sigma_{i=1}^{N}y_i} {N},\frac{\sigma^2} {N})$ であり、正規分布であることが分かる。

次に $\sigma$ を求めるが、これが少し厄介。確率変数xに対し、尤度がexp(xの二次式)の形を取るとき、確率変数xは正規分布の形になるのだが、下記の式の通り、 $\sigma$ だとexpの中身が $\sigma$ の二次式にならない。

\Pi_{i=1}^{N}(\frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y_i-\mu)^2}{2\sigma^2}))

この件は下記のリンク参照した。

では、どう解くか。expの式の場合は、とりあえず対数をとってみるのが定石。また、 $\sigma=1/\sqrt{\tau}$ としておく。理由は後述する。

\begin{align*} \log p(\sigma|\mu,y) &\propto \log p(y|\mu,\sigma) \\ &=N\log (\frac{1} {\sqrt{2\pi\sigma^2}})-\frac{\Sigma_{i=1}^{N}(y_i-\mu)^2} {2\sigma^2} \\ &=-N\log \sigma - \frac{\Sigma_{i=1}^{N}(y_i-\mu)^2} {2\sigma^2} -\frac{N}{2}\log 2\pi \\ \log p(\tau|\mu,y)&=\frac{N}{2}\log \tau - \frac{\tau}{2}\Sigma_{i=1}^{N}(y_i-\mu)^2 - \frac{N}{2} \log 2\pi \end{align*}

ここでガンマ分布 $Gamma(\alpha,\beta)$ の確率密度関数f(x)を見てみよう。

\begin{align*} f(x|\alpha,\beta)&=\frac{\beta^\alpha x^{\alpha-1} exp(-\beta x)} {\Gamma(\alpha)} \\ &\propto \beta^\alpha x^{\alpha-1} exp(-\beta x) \\ \log f(x|\alpha,\beta)&\propto \alpha\log \beta + (\alpha-1)logx - \beta x \end{align*}

ここで $\Gamma(\alpha)$ はガンマ関数であり、 $\Gamma(\alpha+1)=\alpha\Gamma(\alpha)$ という関係が成り立つ。 $\log p(\tau|\mu,y)$ と $\log f(x|\alpha,\beta)$ の形を見比べてみよう。xと $\tau$ が同じだとすると

\frac{N}{2} = (\alpha-1) \\\\ \frac{1}{2}\Sigma_{i=1}^{N}(y_i-\mu)^2 = \beta

よって、 $p(\tau|\mu,y) \sim Gamma(\frac{N}{2}+1,\frac{1}{2}\Sigma_{i=1}^{N}(y_i-\mu)^2)$ となり、ガンマ分布であることが分かる。これで、正規分布とガンマ分布から確率で重みづけしたサンプリングをすれば、ギブスサンプリングができそうだ。この考えをもとに、プログラミングしてみよう。

実践：ギブスサンプリングをやってみよう

繰り返しになるが、こちらのサイトを参考にしている。

今回はモデルの評価や途中データの詳細な確認等は行わないためプログラムから省いているが、興味がある方は上記のサイトを是非参考にしてほしい。

まずは環境構成とディレクトリ構成図。

環境構成

Python：3.9.13
VSCode：1.76.2

ディレクトリ構成図

GIBBS-SAMPLING
├ backend.py
└ frontend.py

まずは、backend.pyのコード。

backend.py

import numpy as np
import pandas as pd

def sample_mu(y, N, sigma):
    #平均
    mean = np.sum(y) / N
    #分散
    variance = sigma * sigma / N
    return np.random.normal(mean, np.sqrt(variance))

def sample_sigma(y, N, mu):
    #アルファ
    alpha = N / 2 + 1
    #残差(yi-mu)
    residuals = y - mu
    #ベータ
    beta = np.sum(residuals * residuals) / 2
    #タウ
    tau = np.random.gamma(alpha, 1/beta)
    return (1 / (np.sqrt(tau)))

def model1(y, iters, init):
    mu = init["mu"]
    sigma = init["sigma"]
    N = len(y)
    
    trace = np.zeros((iters, 2))
    
    for i in range(iters):
        mu = sample_mu(y, N, sigma)
        sigma = sample_sigma(y, N, mu)
        trace[i, :] = np.array((mu, sigma))
    
    trace = pd.DataFrame(trace)
    trace.columns = ['mu', 'sigma']
    
    return trace

np.random.normalで何をやっているか
事後分布である正規分布の確率密度関数f(x)から乱数を発生させている。事後分布が正規分布になるものを用いた理由は、乱数発生をnp.random.normal()だけでできるためである。この乱数はf(x)の値に比例してxを抽出している。（=f(x1)大ならばx1は乱数抽出されやすい）今回は省略しているが、第三引数は乱数抽出数を表す。省略しているので、乱数は1つしか抽出されない。

tau = np.random.gamma(alpha, 1/beta)で何をやっているか
numpyではガンマ関数を下記の通り設定してある。

p(x) = x^{k-1}\frac{e^{-x/\theta}}{\theta^k\Gamma(k)}

理論編で述べたガウス関数と比較すると、 $\theta=1/\beta$ であるので、第二引数は1/betaとしていることに注意。

次にfrontend.pyのコード。

frontend.py

import matplotlib.pyplot as plt
import numpy as np
from backend import model1
import streamlit as st

with st.form(key="gibbs-sampling"):
    iters:int=st.slider("試行回数", 0, 10000, 0, 100)
    size:int=st.slider("現時点のサンプル数", 0, 10000, 0, 10)
    MU:int=st.number_input("サンプルの平均値",0)
    SIGMA:int=st.number_input("サンプルの標準偏差",0)
    burn_in:int=st.number_input("バーンイン数",0)
    submit_button=st.form_submit_button("ギブスサンプリング")
    my_bar = st.progress(0)
    if submit_button:
        plt.rcParams['figure.figsize'] = (16, 8)
        np.random.seed(11)

        init = {"mu": np.random.uniform(-10000, 10000), "sigma": np.random.uniform(0, 10000)}
        my_bar.progress(20)

        y = np.random.normal(MU, SIGMA, size=size)
        trace = model1(y, iters, init)

        fig1 = plt.figure()
        ax1 = fig1.add_subplot(2, 1, 1)
        ax2 = fig1.add_subplot(2, 1, 2)
        trace.plot(ax=ax1)
        ax1.legend(loc='upper right', prop={'size': 15})

        kwargs = dict(histtype='stepfilled', alpha=0.3, density=True, bins=60, ec="k")
        my_bar.progress(40)
        ax2.hist(trace['mu'], label='mu', **kwargs)
        ax2.hist(trace['sigma'], label='sigma', **kwargs)
        ax2.legend(loc='upper right', prop={'size': 15})

        row_data=trace.describe()
        st.markdown("#### データ推移と事後確率密度分布")
        st.pyplot(fig1)
        st.markdown("#### データ表")
        st.write(row_data)

        burn_in_trace=trace[burn_in:]
        fig2 = plt.figure()
        ax3 = fig2.add_subplot(2, 1, 1)
        ax4 = fig2.add_subplot(2, 1, 2)
        burn_in_trace.plot(ax=ax3)
        ax3.legend(loc='upper right', prop={'size': 15})
        my_bar.progress(60)

        kwargs = dict(histtype='stepfilled', alpha=0.3, density=True, bins=60, ec="k")
        my_bar.progress(80)
        ax4.hist(burn_in_trace['mu'], label='mu', **kwargs)
        ax4.hist(burn_in_trace['sigma'], label='sigma', **kwargs)
        ax4.legend(loc='upper right', prop={'size': 15})

        burn_in_data=burn_in_trace.describe()
        st.markdown("#### データ推移と事後確率密度分布（バーンイン込み）")
        st.pyplot(fig2)
        st.markdown("#### データ表（バーンイン込み）")
        st.write(burn_in_data)
        my_bar.progress(100)

itersとsizeは何か
itersはギブスサンプリングの試行回数、つまり求める事後確率密度分布を何個の点で作るかを示している。sizeは現時点で獲得しているサンプル数を示している。

my_bar = st.progress(0)で何をやっているか
プログレスバーの表示を今回は導入。計算に時間がかかるときは、これを入れておくとよいだろう。使い方としては、my_bar = st.progress(0)を入れた後のコードの各部分にmy_bar.progress(20)、my_bar.progress(40)、･･･、と入れていき、最後の処理完了時にmy_bar.progress(100)を入れることでプログレスバーが完成する。

plt.rcParams['figure.figsize'] = (16, 8)で何をやっているか
rcParamsで以降のmatplotlibの設定がすべて変わる。図を複数同じ形で出力するときは初めにplt.rcParams[]を行うとよい。

init = {"mu": np.random.uniform(-10000, 10000), "sigma": np.random.uniform(0, 10000)}で何をやっているか
np.random.uniform(下限値,上限値)を用いることで、初期値を一様分布から乱数発生させている。

本来ならば、 $\mu \sim uniform(-\infty,\infty)$ 、 $\sigma \sim uniform(0,\infty)$ として初期値を設定しなければならないが、np.random.uniform()の引数にnp.float('inf')を使ってしまうと、下記のエラーが発生する。

OverflowError: Range exceeds valid bounds

とりあえず大きい数を引数にすることで暫定処置を行っているが、良い方法がまだ思いつかない。もう少し調査する。

ax1 = fig1.add_subplot(2, 1, 1)で何をやっているか
add_subplot()でfigの表示方法を指定している。第一引数が行の分割数、第二引数が列の分割数を表しており、今回の場合は2×1表示している。第三引数は通し番号を表しているが、ここが1から採番していないと下記のエラーが発生する。

#第三引数に3を入力したとき
num must be an integer with 1 <= num <= 2, not 3

kwargs = dict(histtype='stepfilled', alpha=0.3, density=True, bins=60, ec="k")で何をやっているか
kwargs（ケーダブルアーグス）は複数のキーワード引数を持つデータを辞書型で受け取る際に用いられる。下記のリンクから分かるようにグラフの設定引数は多く存在するため、kwargsを用いることで、見やすくなる。

burn_in_trace=trace[burn_in:]で何をやっているか
序盤に述べたバーンイン部分をカットしている

結果：ギブスサンプリングをやってみよう

初期画面は下のようになるはず。

streamlit初期画面

試しにやってみよう。

実行画面

出力結果①

出力結果②

無限母集団のデータであっても、100個のデータ数だけで母数の確率密度関数が割り出せることが分かった。また、色々なデータで試してみると、当たり前のことかもしれないが下記のことが分かる。

・サンプル数が少ないと母数の確率密度関数の分散が大きくなるので、サンプル数が増えるたびにギブスサンプリングを行うとよい
・ギブスサンプリングの試行回数を増やしてみると母数の確率密度関数がより滑らかになる
・バーンインの範囲が小さくても、良い精度のサンプリングができる。

皆さんも色々なデータで試してみてほしい。

まとめ

ギブスサンプリングの素晴らしさを理解してもらえただろうか。長く書いてしまったが、やっていることは「事前データから求める一つの母数の事後確率を求める→事後確率から疑似母数を発生させる→異なる一つの母数の事後確率を求める→事後確率から疑似母数を発生させる→･･･」である。実サンプルのサンプリングが難しくても、ギブスサンプリングを用いれば疑似サンプルを生み出せるのは良いことである。一つ欠点を挙げるとしたら、事後確率が分かりやすい分布でないと疑似母数を発生させることが難しいことである。これを補っているのがメトロポリス法であるがまた今度。お疲れさまでした。