実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で③離散変数のサンプリング

「PythonのMCMCライブラリPyMC」章

テキスト「PythonのMCMCライブラリPyMC」の執筆者
渡辺祥則 先生

---

この記事は、テキストの「PythonのMCMCライブラリ PyMC」章の例題３「離散変数のサンプリング」の実践を取り扱います。
変化点検出モデルと呼ばれるベイズモデリングにチャレンジします。
PyMC5 でサクッと離散変数を推論しましょう！

はじめに

---

岩波データサイエンスVol.1の紹介

この記事は書籍「岩波データサイエンス vol.1」（岩波書店、以下「テキスト」と呼びます）の特集記事「ベイズ推論とMCMCのフリーソフト」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

岩波データサイエンス Vol.1

テキストは、2015年10月に発売され、ベイズモデリングの様々なソフトウェアを用いたモデリング事例を多数掲載し、ベイズモデリングの楽しさを紹介する素晴らしい書籍です。
入門的なモデルから2次階差を取り扱う空間モデルまで、幅広い難易度のモデルを満喫できます！

このシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「岩波データサイエンス vol.1」
第９刷、編者 岩波データサイエンス刊行委員会 　岩波書店

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

かわいいフリー素材集

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。

PythonとPyMCの最新化
テキストのPyMCのバージョンはVer.3です。
また、テキスト発行の2015年から現在2024年までの約10年で、Python及び各種パッケージの利用方法が変わっています。
そこでこの記事では、テキスト掲載のPython＆PyMC3コードを自己流Python＆PyMC5コードに変換することを、実践事項の１つに取り上げようと思います。

例題３ 離散変数のサンプリング

---

1. 前説

イギリスの炭鉱事故の年間発生件数データを用いてモデリングします。
PyMC公式サイトの「Case study 2: Coal mining disasters」もこのデータを用いたモデリングを実践しています。

Introductory Overview of PyMC — PyMC dev documentation

４０年目あたりで件数減少の傾向が見られるとのこと。
この件数減少＝変化点を検出するモデルを構築します。

2. 準備

「PythonのMCMCライブラリPyMC」章で利用するパッケージをインポートします。

### インポート

# ユーティリティ
import pickle

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# 回帰分析
import scipy.stats as stats
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression

# 分類タスク
from sklearn.metrics import classification_report

# WEBアクセス
import requests

# yfinance株価データ取得
from pandas_datareader import data as pdr
import yfinance as yf

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

3. データの設定

テキストの件数データを実直に実装します。

### データの設定 テキストのPythonスクリプトから引用
data3 = np.array(
    [4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6, 3, 3, 5, 4, 5, 3,
     1, 4, 4, 1, 5, 5, 3, 4, 2, 5, 2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1,
     1, 3, 0, 0, 1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1, 0, 1,
     0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1,
     1, 2, 4, 2, 0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1]) 

print('data3.shape: ', data3.shape)

【実行結果】
111件＝111年分の時系列データです。

4. データの外観の確認

時系列にプロットします。
テキストの図6に相当します。

### データの可視化 ※図6に相当
plt.figure(figsize=(8, 4))
ax = plt.subplot()
ax.plot(range(data3.shape[0]), data3, '-o', ms=3)
ax.set(xlabel='年', ylabel='件数', title='イギリスの炭鉱事故数の推移',
       xlim=(-5, 121))
ax.grid(lw=0.5)

【実行結果】
確かに40年あたりで発生件数が減少している感じがします。

階層ベイズモデル

---

1. 前説

■ モデル数式
テキストのPyMCモデルから数式を逆引きいたします。
事前分布から尤度関数の順で書きます。

$$
\begin{align*}
switchpoint &\sim \text{DiscreteUniform}\ (\text{lower}=0,\ \text{upper}=years) \\
 \\
earlyMean &\sim \text{Exponential}\ (\text{lam}=1) \\
lateMean &\sim \text{Exponential}\ (\text{lam}=1) \\
means = &\text{pt.switch}\ (yearIdx < switchpoint,\ earlyMean,\ lateMean) \\
 \\
y &\sim \text{Poisson}\ (\text{mu}=means) \\
\end{align*}
$$

$${switchpoint}$$は変化点です。興味のあるパラメータです。
0から111年を範囲にする離散一様分布に従います。
つまり、離散変数なのです。

$${earlyMean,\ lateMean}$$は変化点前と変化点後の発生件数の平均です。
指数分布に従います。指数分布は０以上の連続値を扱う分布です。
$${means}$$（テキストのrate）は111年分の平均をまとめて扱う変数です。
変化点前の年は$${earlyMean}$$が設定され、変化点後の年は$${lateMean}$$が設定されるように、pytensor.tensor の switch 関数で計算します。

目的変数である発生件数$${y}$$は、平均パラメータに$${means}$$を用いるポアソン分布に従います。ポアソン分布は０以上の整数値を扱う分布です。

2. モデルの定義

数式モデルどおり実直にPyMCのモデルを定義しましょう。

### モデルの定義

## データの定義
num_years = len(data3)
year_idx = list(range(num_years))

## モデルの定義
with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    # 観測データのインデックス
    model3.add_coord('data', values=year_idx, mutable=True)
    
    ## dataの定義
    # 観測値Y :事故数
    Y = pm.ConstantData('Y', value=data3, dims='data')
    # 年のインデックス
    yearIdx = pm.ConstantData('yearIdx', value=year_idx, dims='data')

    ### 事前分布
    # 変化点
    swithpoint = pm.DiscreteUniform('switchpoint', lower=0, upper=num_years)
    # 変化点前の平均
    earlyMean = pm.Exponential('earlyMean', lam=1)
    # 変化点後の平均
    lateMean = pm.Exponential('lateMean', lam=1)
    # 全年の平均の要素
    means = pm.Deterministic(
        'means',
        pt.switch(yearIdx <= swithpoint, earlyMean, lateMean),
        dims='data')
    
    ### 尤度:ポアソン分布
    y = pm.Poisson('y', mu=means, observed=Y, dims='data')

モデルの内容を表示・可視化してみましょう。
モデル名 model3 とタイプするとモデルの数式の概要が表示されます。

### モデルの表示
model3

【実行結果】

続いてモデルを可視化します。
１行のコードでグラフ形式の可視化ツール「Graphviz」に自動連携して、変数の関係を図示してくれます！

### モデルの可視化
pm.model_to_graphviz(model2)

【実行結果】
データやパラメータの繋がりが示されます。
丸い形状の変数には「確率分布」が併記されます。
データ値を示す変数はグレー色が付いています。

3. MCMCの実行と収束の確認

■ MCMC
いよいよベイズプログラミングの真骨頂 MCMC の実行です。
マルコフ連鎖=chainsを4本、バーンイン期間=tuneを1000、利用するサンプル=drawを１chainあたり1000に設定して、合計4000個の事後分布からのサンプリングデータを生成します。テキストの指定数とは異なっています。

サンプリングデータは idata3 に格納されます。

このモデルは離散変数の事前分布を含んでいるため、サンプル方法に numpyro を指定できません。

### 事後分布からのサンプリング 1分10秒
with model3:
    idata3 = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1234,
                       target_accept=0.9)

【実行結果】
処理時間は約1分10秒です。
サンプル方法は、離散変数の switchpoint にPyMC標準の Metropolis が適用、その他の連続変数に PyMC標準の NUTS が適用されました。

■ $${\hat{R}}$$で収束確認
収束の確認をします。
指標$${\hat{R}}$$（あーるはっと）の値が$${1.1}$$以下のとき、収束したこととします。
次のコードではしきい値を$${1.01}$$にキツく設定して、$${\hat{R}>1.01}$$のパラメータの個数をカウントします。
個数が０ならば、$${\hat{R} \leq1.01}$$なので収束したとみなします。

### r_hat>1.1の確認
# 設定
idata_in = idata3        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R}>1.01}$$のパラメータは０個でした。
$${\hat{R} \leq1.1}$$なので収束したとみなします。

■ $${\hat{R}}$$の値把握と事後統計量の表示

### 事後統計量の表示
var_names = ['beta', 's', 'r', 'q']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names)

【実行結果】
MCMCでサンプリングデータを生成したすべてのパラメータについて、統計量や診断情報が一覧表示されました。

右端の$${\hat{R}}$$（r_hat）はすべて$${1.0}$$です。
最上行が注目するパラメータ$${switchpoint}$$です。
事後分布の統計量は、平均$${38.902}$$、95%HDI区間$${[34, 44]}$$です。
少々幅広ですが、変化点は39年前後、39年までが変化点前、40年以降が変化点後という感じです。

■ トレースプロットで収束確認
トレースプロットを描画します。
テキストの図7に該当します。
こちらの図でも収束の確認を行えます。

### トレースプロットの描画 ※図7に相当
pm.plot_trace(idata3, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】
左の曲線は、４本のマルコフ連鎖の分布です。
４本の曲線がほぼ同じ分布を描いているので、サンプリングデータが同一の分布から生成されたことが推察できます。
右のノイズみたいな図は、均等に乱雑に描画されています。
この状態は収束していると推察できます。
均等・乱雑ではなくて、何らかの傾向を示すような描画の場合は収束を疑います。

$${\hat{R}}$$とトレースプロットから、分布の収束が確認できましたので、MCMCのサンプリングデータを用いて、推論を継続できます！
やったね！

■ 事後予測プロットでモデルの良し悪しを確認
もうひとつのプロットで推論の良し悪しを確認してみましょう。
目的変数$${y}$$の事後予測サンプリングを行ってグラフで可視化します。
事後分布からのサンプリング（MCMCのこと）では事後分布を得て「パラメータ」を推論します。
事後予測では、得られた事後分布を元にして「目的変数」を予測しちゃおうってことなのです！

### 事後予測サンプリングの実行
with model3:
    idata3.extend(pm.sample_posterior_predictive(idata2, random_seed=1234))

【実行結果】

4000個のサンプリングデータのうち100個をランダムに選んで描画します。

### 事後予測プロットの描画
pm.plot_ppc(idata3, num_pp_samples=100, random_seed=1234);

【実行結果】
発生件数$${y}$$は階段状で描画されています。

黒線が目的変数$${Y}$$の観測値、オレンジの点線が目的変数の事後予測の平均値です。
平均値は小数もありなので、整数値の観測値とズレているのかもしれません。

まとめますと、今回のベイズモデルは目的変数をまあまあ表現できていると思います！

4. 分析

■ 平均値$${means}$$の推論値のプロット
MCMCのサンプリングデータに含まれる「ポアソン分布の平均パラメータ」$${means}$$は発生件数の変化の傾向を捉えています。
プロットしましょう。テキストの図8に相当します。
サンプリングデータ4000個から400系列を取り出して描画します。

### meansの事後分布プロット ※図8に相当

## データの準備
# 変化点前後の平均meansの取り出し
means_samples = idata3.posterior.means.stack(sample=('chain', 'draw')).data
# 変化点の事後平均値
switchpoint_mean = idata3.posterior.switchpoint.data.flatten().mean()

## 描画
# 描画領域の指定
fig, ax = plt.subplots(figsize=(8, 4))
# meansの事後分布サンプリングデータを1サイクルづつ繰り返し描画（10ごとの間引き）
for i in range(0, means_samples.shape[1], 10):
    # 折れ線グラフの描画
    ax.plot(year_idx, means_samples[:, i], alpha=0.05, color='tab:blue')
# 観測値のドットプロット
ax.plot(range(data3.shape[0]), data3, 'o', ms=5, color='red', alpha=0.5,
        label='観測値')
# 変化点の事後平均値の垂直線の描画
ax.axvline(switchpoint_mean, color='black', ls='--',
           label='変化点の平均値')
# 修飾
ax.set(xlabel='年', ylabel='件数',
       title='ポアソン分布の平均パラメータの事後予測プロット')
ax.grid(lw=0.5)
ax.legend()
plt.show()

【実行結果】
400系列の推論値の青い線の束はクランクのような形状です。
黒い点線は変化点$${switchpoint}$$の事後分布の平均値です。
平均$${means}$$の推論値は変化点の平均値前後で減少しています。
赤点の観測値の中心あたりを青い線が捕捉しているように感じます。

5. 推論データの保存

pickle で MCMCサンプリングデータ idata3 を保存できます。

### 推論データの保存
file = r'idata3_ch1.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata3, f)

次のコードで保存したファイルを読み込みできます。

### 推論データの読み込み
file = r'idata3_ch1.pkl'
with open(file, 'rb') as f:
    idata3_ch1_load = pickle.load(f)

ベイズ記事は以上です。

次回予告

「PythonのMCMCライブラリ PyMC」章
　例題４「シークレット」

結び

---

変化点検出モデルは多くの書籍で紹介されているベイズモデルの代表例です。
Stanは離散変数のパラメータをそのまま扱えず、「周辺化」というテクニックを利用するそうです。
離散変数のその他の代表例には「隠れマルコフモデル」というモデルがあるようです。
さまざまなモデリングのノウハウが先達の活動によって蓄積されているのでしょう。
ベイズの世界は奥が深いですね。
ベイズがわかった！と言える日がくることを信じて、ベイズを楽しもうと思います。

おわり

---

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

---

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

１．のんびり統計

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

２．実験！たのしいベイズモデリング１＆２を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の教科書です。
よかったらぜひ、お試しくださいませ。

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