「ベイズ統計モデリングによるデータ分析入門」をPythonで写経 ～ Vol.1 MCMCの基本

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第１部第７章「MCMCの基本」の Python 写経活動記録です。

Python による最初の実装は MCMC です！ヽ(=´▽`=)ﾉ ﾜｰｲ

MCMC は「マルコフ連鎖モンテカルロ法」の略称です。
ベイズ統計モデリングにとって無くてはならない存在です！

書籍ではコード化されていませんが、MCMC の動きを体感したくなって、MCMC の一種「ランダムウォーク MH 法」を Python でコード化しました。

では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

■ 記事の範囲
この記事はテキスト第１部第７章の以下の節を取り扱います。

■ ライブラリのインポート
Jupyter Notebook 形式でコードを記述します。
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib 等

ランダムウォークMH法の実装

テキスト p.69 のランダムウォークMH法の手計算部分を参考にして、MCMC の手続きを想定し、Python コード化して MCMC を動かします！

MCMC はベイズ統計モデリングの中核機能だと思います。
この記事で「MCMC の雰囲気」を体感していきます！

かなりざっくり、ベイズの定理とMCMC

🔢 ベイズの定理
ベイズの定理を数式とベイズ用語で確認しましょう。
テキスト第１部第１章の記述を引用・参考にしています。
$${\theta}$$ は推定対象のパラメータ、$${X}$$ はデータです。

上式の１行目はベイズの定理です。
テキストの式（1.50）を一部改変しています。
２行目はベイズの定理の各項にベイズ用語を付記したものです。
テキストの式（1.51）を一部改変しています。
３行目は周辺尤度を取り除いてカーネルに焦点をあてた数式です。
テキストの式（1.63）の一部分を取り出して一部改変しています。

３行目の左辺に相当する「パラメータ $${\theta}$$ の事後分布」を推定するために、右辺の「尤度」と「パラメータ $${\theta}$$ の事前分布」で構成される「カーネル」を利用するのが、ベイズ統計モデリングのポイントです。

🔢 MCMC
MCMC はざっくり、カーネル比などを用いてパラメータの事後分布からの乱数生成を行なう手法です。
事後分布自体を直接計算するのではなく、取得した乱数を通じて事後分布を理解するような感じです（遠回りな感じ）。
MCMC のメリットは計算が大変な「周辺尤度（正規化定数）」を求める必要がないことです。
カーネル比を計算することは、周辺尤度を打ち消して無くすことなのです。

🔢 例題のモデルの確認
テキスト第７章の例題（6.10節の計算例）をお借りします。
標本サイズ $${5}$$ の売り上げデータ $${X}$$ が平均 $${\theta}$$、分散 $${1}$$ の正規分布に従う、という確率モデルを想定します。

データ $${X}$$ が与えられたときのパラメータ $${\theta}$$ の事後分布を推定したいです。

平均 $${\theta}$$、分散 $${\sigma^2}$$ の正規分布の確率密度関数は次のとおりです。

尤度 $${f(x \mid \theta)}$$ は 標本サイズ $${5}$$ のデータ $${X}$$ 各要素の確率密度の積であり、次のようになります。

$${\prod}$$ は各要素を「掛け算」する記号「総乗記号」です。
$${\prod_{i=1}^5}$$ は $${x_1, x_2, x_3, x_4, x_5}$$ を右側の数式の$${x_i}$$ に代入してできる「５つの項」を掛け算しています。

パラメータ $${\theta}$$ の事前分布 $${f(\theta)}$$ は平均 $${0}$$、分散 $${10000}$$ の正規分布に従うとし、確率密度関数は次のようになります。

カーネル（$${=}$$ 尤度 $${\times}$$ 事前分布）は次のようになります。

このカーネルを用いて、テキスト第７章の MCMC のステップを Python で書き、動かしてみたいと思います！

🔷🔷🔷

関数の定義

MCMC 実装（実験）に向けて、Python の関数を３つ定義します。
こちらはカーネル関数です。
上式のカーネルの数式どおりに実装しています。

# p.69～ MH法の計算例

# カーネルの関数定義 p.69のカーネルの数式
# 引数 θ: パラメータ、X: データ
def kernel(θ, X):
    return (
        np.prod([1 / np.sqrt(2*np.pi) * np.exp(-(xi - θ)**2 / 2) for xi in X])
        * (1 / np.sqrt(20000*np.pi) * np.exp(-θ**2 / 20000))
    )

続いてランダムウォークMH法の実行関数です。
p.69 の「MH法をやってみましょう。」に続く計算過程を参考にして関数化しました。
テキストの計算過程とこのコードを照らし合わせることで、ランダムウォークMH法のロジックの理解につながると思います！

# ランダムウォークMH法の関数定義 p.69のMH法の手順をコード化
# 引数 X: データ、sigma2: 正規分布の分散パラメータ、iter: 乱数生成数、seed: 乱数シード
def randomwalk_MH_method(X, sigma2, iter, seed=None):

    ## 設定
    # 乱数生成器
    rng = np.random.default_rng(seed=seed)
    # 結果を格納する配列の初期化
    θs = np.zeros(iter)
    
    ## 乱数の初期値 θ1 の生成
    θs[0] = rng.uniform(low=-2, high=2)

    ## ２つ目以降の乱数生成
    for i in range(1, iter):
        # １つ前の乱数値の取得
        θ_prev = θs[i - 1]
        # 提案値の生成
        θ_teian = θ_prev + rng.normal(loc=0, scale=np.sqrt(sigma2))
        # カーネル比の算出
        rate = kernel(θ_teian, X) / kernel(θ_prev, X)
        # 乱数値の選択
        if rate > 1:  # カーネル比が１超の場合は提案値を乱数値に設定
            θs[i] = θ_teian
        else:         # カーネル比が１以下の場合は、カーネル比の確率で乱数値を選択
            θs[i] = rng.choice(a=[θ_teian, θ_prev], p=[rate, 1 - rate])
    
    ## 戻り値：パラメータθの事後分布に従う乱数
    return θs

最後にMCMCサンプルを描画する関数です。
テキストの図 1.7.2 ～ 図 1.7.4 のトレースプロットの描画に利用します。
内部でランダムウォークMH法の実行関数を呼び出しています。

# 事後分布乱数の描画関数の定義
# 引数 X: データ、sigma2: 正規分布の分散パラメータ、iter: 乱数生成数
#     lw: 描画する線の太さ、seed: 乱数シード

def trace_plot(X, sigma2, iter, lw=1.5, seed=123):
    # 描画領域の設定
    plt.figure(figsize=(10, 3))
    # パラメータθの事後分布から生成した乱数の軌跡の描画
    plt.plot(randomwalk_MH_method(X, sigma2=sigma2, iter=iter, seed=seed), lw=lw)
    # 修飾
    plt.ylabel(r'$\theta$', fontsize=12)
    plt.grid(lw=0.5, alpha=0.5)
    plt.show()

🔷🔷🔷

ランダムウォークMH法によるMCMCサンプリング

では MCMC の実行に移りましょう。
売上データ $${X}$$ を設定します。テキストのデータを引用いたします。

# p.69 観測データの設定
X = np.array([2.4, 3.2, 2.2, 4.6, 3.3])
print(f'標本平均 = {X.mean():.3f}, 標本分散 = {X.var(ddof=1):.3f}')

【実行結果】
標本平均 3.14、標本分散 0.90 です。
正規分布の平均パラメータ、分散パラメータの目安にどうぞ。

分散１で50個の乱数を生成して「トレースプロット」を描画しましょう。
トレースプロットは「乱数を生成順にプロット」するチャートです。
テキストの図 1.7.2 に相当します。
なお、乱数がテキストと異なるため、折れ線グラフの形状などの結果も異なることに留意して下さい。

# p.70 図1.7.2 MH法による乱数生成を50回繰り返した例

## 設定
# 提案値が従う正規分布の分散
sigma2 = 1
# 乱数の生成数
iter = 50

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter)

【実行結果】
テキストと同様に「20回目以降は、3前後をとる乱数が得られて」います。

続いて 2000 個の乱数を生成してトレースプロットを描画します。
分散は１です。
テキストの図 1.7.3 に相当します。

# p.70 図1.7.3 MH法による乱数生成を2000回繰り返した例

## 設定
# 提案値が従う正規分布の分散
sigma2 = 1
# 乱数の生成数
iter = 2000

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter, lw=1)

【実行結果】

テキストと同様に「100回を超えたあたりから、安定した変動を繰り返して」います（テキスト p.70 より引用）。
同じ範囲をゲジゲジと描画する様子を見て「定常分布に収束していることがうかがえます」（テキスト p.70 より引用）。

今度は 分散 10000、乱数生成数 2000 個 のトレースプロットです。
テキストの図 1.7.4 の上側のチャートに相当します。

# p.71 図1.7.4 分散10000を指定してMH法を実行

## 設定
# 提案値が従う正規分布の分散
sigma2 = 10000
# 乱数の生成数
iter = 2000

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter)

【実行結果】
折れ線グラフは直線的になり、「ゲジゲジ」が見当たりません。
テキストは「$${\hat{\theta}_t}$$ がなかなか変化せず、ずっと同じ値を取り続ける」としています。

最後に 分散 0.0001、乱数生成数 2000 個 のトレースプロットです。
テキストの図 1.7.4 の下側のチャートに相当します。

# p.71 図1.7.4 分散0.0001を指定してMH法を実行

## 設定
# 提案値が従う正規分布の分散
sigma2 = 0.0001
# 乱数の生成数
iter = 2000

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter, lw=1)

【実行結果】
ゲジゲジに到達する前段階です。収束に至っていません。
テキストは「提案値がなかなか変化」せず、「初期値から事後確率密度が高い領域まで遷移するのに、長い時間がかかる」としています。

🔷🔷🔷

収束の確認

テキスト 7.18 節の収束判定のためのトレースプロットを描画しましょう。
乱数を４セット（chains = 4）生成して、トレースプロットを描きます。

複数の chain を描くための関数を作ります。
引数 chain で与えられた本数の乱数セットをランダムウォークMH法で生成して、トレースプロットを描画します、

# 事後分布乱数をchain個描画する関数の定義

def trace_plot(X, sigma2, iter, chain=4, lw=1.5, seed=123):
    # 描画領域の設定
    plt.figure(figsize=(10, 3))
    # パラメータθの事後分布から生成した乱数の軌跡の描画
    for i in range(chain):
        plt.plot(randomwalk_MH_method(X, sigma2=sigma2, iter=iter, seed=seed+i),
                 lw=lw)
    # 修飾
    plt.ylabel(r'$\theta$', fontsize=12)
    plt.grid(lw=0.5, alpha=0.5)
    plt.show()

４本の chain を生成して、トレースプロットを描画しましょう。
まず、テキストの図 1.7.6 左側の「収束したとき」の実行です。

# p.75 図1.7.6 収束したとき

## 設定
# 提案値が従う正規分布の分散
sigma2 = 1
# 乱数の生成数
iter = 2000

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter, chain=4, lw=0.5, seed=1)

【実行結果】
４本の chain は色を変えて描画しています。
割と早いタイミングで $${\theta=3}$$ 付近の「ゲジゲジ」が出現しました。

最後にテキストの図 1.7.6 右側の「収束しなかったとき」の実行です。

# p.75 図1.7.6 収束しなかったとき

## 設定
# 提案値が従う正規分布の分散
sigma2 = 0.0001
# 乱数の生成数
iter = 2000

## トレースプロットの描画
trace_plot(X, sigma2=sigma2, iter=iter, chain=4, lw=1, seed=1)

【実行結果】
４本の chain が混ざり合うこと無く、収束していないことがわかります。

以上でランダムウォークMH法のMCMCサンプリングは終了です。
楽しかったですね！

パラメータの事前分布と尤度を定式化して MCMC を実行し、事後分布からのサンプル＝MCMCサンプルを得る、この流れはベイズ統計モデリングの手続きの「ルーチン」です！
これからお世話になる「ルーチン」です！

そしてそして、今後繰り広げるベイズ統計モデリングでは、トレースプロットによる収束確認が頻出しますので、ぜひお楽しみに🍀

ちなみに、PyMC や Stan には MCMC サンプリング機能が備わっていますので、今回のように MCMC アルゴリズムを実装する必要はありません。
トレースプロットも既存ライブラリの専用機能を使って「サクッと」描けます。
ご安心下さい！

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で８つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．統計・データ分析とつながる

シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀

３．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

４．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

５．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

６．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

７．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

８．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。