「Stan入門」章

テキスト「Stan入門」の執筆者
松浦健太郎 先生

この記事は、テキストの「Stan入門」章の例題３「空間構造のあるベイズモデル」の実践を取り扱います。
2次元・2階階差のCARモデルです。記事第９回・10回に通じます！
テキストの最難関のモデルです（個人の感想です）。
たのしくPyMCモデリングを進めましょう！

この章を執筆された松浦先生は「StanとRでベイズ統計モデリング (Wonderful R 2)」を執筆されています。
こちらもお楽しみくださいませ。

はじめに

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

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

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

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

引用表記

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

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

PyMC環境の準備

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

PythonとPyMC
テキストで利用するツールは、R、Stan です。
この記事では、PythonとPyMCを用いたコードに変換してベイズモデリングを実践いたします。

イントロ

1. モデルのイメージ

■ データ概要
生物の実験を想定します。
縦16行・横24列・合計384個の穴があるプラスチックプレート「384 plate」を用いて、どの処理が有効かどうかのデータ$${Y}$$を得ます。
$${96}$$種類の処理を 384 plate で４回ずつ処理して繰り返し数４のデータを得ています。

■ モデル
２つ隣までのプレートの位置の影響（近隣の影響）を考慮する「２階階差のCARモデル」を利用します。
数式はテキストでご確認ください！

2. インポート

「Stan入門」章で利用するパッケージをインポートします。

### インポート

# ユーティリティ
import pickle

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

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

# 最適化
from scipy.optimize import minimize_scalar

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

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

3. データの読み込み

CSVファイルをpandasのデータフレームに読み込みます。

### 観測値Yのデータの読み込み
data3 = pd.read_csv('./data/data-plate.txt', header=None)
print('data3.shape: ', data3.shape)
display(data3.head())

【実行結果】
16 X 24 の格子状に観測値$${Y}$$が格納されています。

4. データの外観の確認

観測値$${Y}$$のヒートマップを描画します。

### 観測値のヒートマップの描画
fig, ax = plt.subplots(figsize=(12, 5))
sns.heatmap(data3, ax=ax, cmap='Oranges', annot=True, fmt='.1f',
            xticklabels=range(1, data3.shape[1]+1),
            yticklabels=range(1, data3.shape[0]+1),
            cbar_kws={'label': 'Y'})
ax.set_xlabel('Plate Column', fontsize=14)
ax.set_ylabel('Plate Row', fontsize=14)
plt.tight_layout();

【実行結果】
11行6列あたりなど幾つかのホットスポット（？）が見られます。

5. もう１つのデータの読み込み

プレートのどの位置にどの処理を行ったのかを表すデータを読み込みます。

### 処理Tのデータの読み込み　96種類の処理×4を割付けたplateデータ
data3T = pd.read_csv('./data/data-plate-design.txt', header=None)
print('data3T.shape: ', data3T.shape)
display(data3T.head())

【実行結果】
プレートの１行目・１列目の処理は処理番号88、のように読みます。

では、PyMCモデリングに進みます。

ベイズモデル

1. モデルの定義

２階階差のCARモデルをどのように組めばよいのか・・・
以前のICARモデルでは重みが１階分しか使えそうにないし・・・
という事で、変数$${r}$$について、PyMC定番の「全組み合わせを１つずつ事前分布を書くモデル」で描いてみました。for文の箇所です。
なお、変数$${r}$$に設定した事前分布（平均パラメータの値）は正確ではない可能性があります。予めご了承くださいませ。

### モデルの定義

## パラメータの設定
n_row = data3.shape[0]        # plateの行数
n_col = data3.shape[1]        # plateの列数
row_idx = list(range(n_row))  # plateの行インデックス
col_idx = list(range(n_col))  # plateの列インデックス

## モデルの定義
with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    # 観測データ・処理割付データのインデックス
    model3.add_coord('row', values=data3.index, mutable=True)
    model3.add_coord('col', values=data3.columns, mutable=True)
    # 処理パターン
    model3.add_coord('treat', values=range(1, data3T.max().max()+1),
                     mutable=True)
    
    ## dataの定義
    # 観測値
    y = pm.MutableData('y', value=data3.values, dims=('row', 'col'))
    # 処理割付
    T = pm.MutableData('T', value=data3T.values, dims=('row', 'col'))

    ### 事前分布
    ## 標準偏差たち
    sigmaBeta = pm.Uniform('sigmaBeta', lower=0, upper=100)
    sigmaGamma = pm.Uniform('sigmaGamma', lower=0, upper=100)
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=100)
 
    ## 処理の影響 β
    beta = pm.Normal('beta', mu=0, sigma=sigmaBeta, dims='treat')
    
    ## プレートの位置の影響 γ
    # 16x24のリストを準備
    r =  [[0 for j in range(n_col)] for i in range(n_row)]

    # 2次階差を取れない0行,1行 x 0列,1列の事前分布を定義
    for i in range(2):
        for j in range(2):
            r[i][j] = pm.Normal(name='r' + str(i) + '.' + str(j),
                                mu=0, sigma=sigmaGamma)
    # 2次階差を取れるセルの事前分布を定義
    for i in range(n_row):
        for j in range(n_col):
            if (i<2)&(j<2):     # 2次階差を取れないケース
                pass
            elif (i<2)&(j>=2):  # 列jのみ2次階差を取れるケース
                r[i][j] = pm.Normal(
                    name='r' + str(i) + '.' + str(j),
                    mu=2*r[i][j-1] - r[i][j-2],
                    sigma=sigmaGamma)
            elif (i>=2)&(j<2):  # 行iのみ2次階差を取れるケース
                r[i][j] = pm.Normal(
                    name='r' + str(i) + '.' + str(j),
                    mu=2*r[i-1][j] - r[i-2][j],
                    sigma=sigmaGamma)
            else:               # 行iと列jの両方が2次階差を取れるケース
                r[i][j] = pm.Normal(
                    name='r' + str(i) + '.' + str(j),
                    mu=(2*r[i][j-1] - r[i][j-2] + 2*r[i-1][j] - r[i-2][j]) / 2,
                    sigma=sigmaGamma)
    # rのリストからgammaを作成
    gamma = pt.stacklists(r)

    ### 尤度
    Y = pm.Normal('Y', mu=gamma + beta[T[row_idx]-1], sigma=sigmaY, observed=y,
                  dims=('row', 'col'))

モデルの内容を可視化してみましょう。

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

【実行結果】
モデルの一部分です。やばい感満載です。

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

■ MCMC
ではMCMC実行！

### 事後分布からのサンプリング 6分ほど処理した後にエラー終了します
# テキスト：iter=5200, warmup=200, thin=5, chains=3
with model3:
    idata3 = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1234,
                       nuts_sampler='numpyro', target_accept=0.8)

【実行結果】
処理時間６分経過したあたりで次の表示が出ます（抜粋）。

モデルのコンパイルに用いている「g++.exe」の実行時に発生したエラーのようです。
stack overflow 等のWebサイトの過去問い合わせでは、「Theano」「Keras」に問題がある事例（両パッケージの再インストールでエラー解消）が掲載されていました。
TheanoはPyMC3でバックエンド処理に利用していたパッケージです。
TheanoもKerasもインストールしていないし、どうしたらいいのでしょう？
無念ですが、ここで記事を終了します（強制終了）

【緩募】
もしもエラー回避コード、または、Python＆PyMCの実行可能コードをご存じの方がいらっしゃいましたら、そのコードを教えてください。

ベイズ記事は以上です。

ちなみに、データ量を減らした実装を次の記事に掲載しています。
よかったらご一読くださいませ。

次回予告

今回が最終回です。

結び

最後の最後で強制終了。
なんだか私らしくて（苦笑い）

おわり

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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