見出し画像

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第12章「12.6 １次元の空間構造」

ネイピア DS

第12章「時間や空間を扱うモデル」

書籍の著者　松浦健太郎 先生

この記事は、テキスト第12章「時間や空間を扱うモデル」の 12.6節「１次元の空間構造」の PyMC5写経 を取り扱います。
１次元の空間構造の表現には PyMC の ICAR() を用いました。
たぶん、テキストのモデルと異なっています。

なお、12.4節「その他の拡張方法」、12.5節「時間構造と空間構造の等価性」は写経を省略いたしました。

はじめに

StanとRでベイズ統計モデリングの紹介

この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。

テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！

テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版

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

PyMC環境の準備

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


12.6 １次元の空間構造

モデリングの準備

インポート

### インポート

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

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

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

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

データの読み込みと確認
サンプルコードのデータを読み込みます。
１次元で等間隔に隣接し合う、ある植物の個体数データです。

### データの読み込み ◆data-kubo11a.txt
# Y:ある植物の個体数（インデックスは1次元区間上に等間隔に並ぶ50の観測地点を示す）

data = pd.read_csv('./data/data-kubo11a.txt')
print('data.shape: ', data.shape)
display(data.head())

【実行結果】

要約統計量を算出します。

### 要約統計量の表示
data.describe().round(2)

【実行結果】

棒グラフで可視化してみます。

### 棒グラフの描画
plt.figure(figsize=(10, 4))
plt.bar(range(1, 51), data['Y'], alpha=0.7)
plt.xlabel('観測地点')
plt.ylabel('個体数 $Y$')
plt.grid(lw=0.5);

【実行結果】
x軸が１次元の空間です。
隣接する地点の個体数は似ているような似ていないような、そんな感じです。

モデル式12-11：ICARランダム効果モデル

PyMCでモデル式12-11を実装します。
ICARを用いたランダム効果モデルです。

PyMCのモデル定義

隣接行列を作成します。
隣り合う地点に１を立てます。

### 隣接行列Wの作成
W = np.zeros((len(data), len(data)))
for i in range(1, len(data)):
    W[i, i-1] = 1
    W[i-1, i] = 1
print('W.shape: ', W.shape)
print(W)

【実行結果】
行が１次元の各地点を示し、列が「お隣さん」を示します。
列に１が立っている地点（インデックス）がお隣さんです。

モデルの定義です。 

### モデルの定義 ◆モデル式12-11 model12-11.stan ※ICAR利用

with pm.Model() as model1:

    ### データ関連定義
    # coordの定義
    model1.add_coord('data', values=data.index, mutable=True)
    # dataの定義
    y = pm.ConstantData('y', value=data.Y.values, dims='data')

    ### 事前分布
    # 標準偏差
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    # 空間構造を取り入れたランダム効果Sのモデル
    # 非中心配置のパラメータ化：ICARの標準偏差を1にして、結果のSにsigmaを乗ずる
    S = pm.ICAR('S', W=W, sigma=1, dims='data')
    r = pm.Deterministic('r', S*sigma, dims='data')

    ### 尤度
    Y = pm.Poisson('Y', mu=pt.exp(r), observed=y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model1

【実行結果】

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

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 20秒
with model1:
    idata1 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['sigma', 'S', 'r']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata1, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】
Sとrに外れ値？的な推論データが含まれているようです。

推論結果の解釈

推論データを用いてテキスト図12.7左に相当するグラフを描画します。
描画関数を作成します。

### Yの観測値rの事後分布を描画する関数の定義

def plot_pred(idata):
    ## 描画用データの作成
    # 推論データからrのMCMCサンプルデータを取り出して指数をとる
    exp_r_samples = np.exp(idata.posterior.r.stack(sample=('chain', 'draw')).data)
    # exp(r)の中央値・80%CI・50%CIの算出
    exp_r_median = np.median(exp_r_samples, axis=1)
    exp_r_80ci = np.quantile(exp_r_samples, q=[0.10, 0.90], axis=1)
    exp_r_50ci = np.quantile(exp_r_samples, q=[0.25, 0.75], axis=1)
    # x軸の値の算出
    xvals = list(range(1, len(data)+1))

    ## 描画処理
    # 描画領域の設定
    fig, ax = plt.subplots(figsize=(7, 4))
    # Yの観測値の描画
    ax.plot(xvals, data['Y'], 'o', ms=7, color='blue', alpha=0.5, 
            label='$Y$ の観測値')
    ax.plot(xvals, exp_r_median, color='tab:red',
            label='$exp[r]$ の事後分布・中央値')
    ax.fill_between(xvals, exp_r_50ci[0], exp_r_50ci[1], color='tomato',
                    alpha=0.4, label='$exp[r]$ の事後分布・50%CI')
    ax.fill_between(xvals, exp_r_80ci[0], exp_r_80ci[1], color='tomato',
                    alpha=0.2, label='$exp[r]$ の事後分布・80%CI')
    ax.set(xlabel='観測地点 $i$', ylabel='個体数 $Y[i]$')
    ax.grid(lw=0.5)
    ax.legend(bbox_to_anchor=(1, 1))
    plt.show()

描画します。

### rの推定結果の描画 ◆図12.7左
plot_pred(idata1)

【実行結果】
テキストと比べるとベイズ信頼区間が広い感じがします。
別のモデルを検討しましょう。

モデル式12-11：切片β＋ICARランダム効果モデル

こちらもPyMCでモデル式12-11を実装します。
切片とICARを用いたランダム効果のモデルです。
前のモデルに切片$${\beta}$$を加えています。

PyMCのモデル定義

モデルの定義です。 
正規分布に従う切片 beta を追加しています。

### モデルの定義 ◆モデル式12-11 model12-11.stan ※岩波データサイエンスvol1参考

with pm.Model() as model1_2:

    ### データ関連定義
    # coordの定義
    model1_2.add_coord('data', values=data.index, mutable=True)
    # dataの定義
    y = pm.ConstantData('y', value=data.Y.values, dims='data')

    ### 事前分布
    # 切片
    beta = pm.Normal('beta', mu=0, sigma=100)
    # 標準偏差
    sigma = pm.Uniform('sigma', lower=0, upper=100) 
    # 切片β＋空間構造を取り入れたランダム効果Sのモデル
    # 非中心配置のパラメータ化：ICARの標準偏差を1にして、結果のSにsigmaを乗ずる
    S = pm.ICAR('S', W=W, sigma=1, dims='data')
    r = pm.Deterministic('r', beta + S*sigma, dims='data')

    ### 尤度
    Y = pm.Poisson('Y', mu=pt.exp(r), observed=y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model1_2

【実行結果】

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

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model1_2:
    idata1_2 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['beta', 'sigma', 'S', 'r']
pm.summary(idata1_2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata1_2, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】
前のモデルと比べるとキレイな分布になっていると思います。

推論結果の解釈

推論データを用いてテキスト図12.7左に相当するグラフを描画します。

### rの推定結果の描画 ◆図12.7左
plot_pred(idata1_2)

【実行結果】
テキストに近いプロットだと思います。
前のモデルと比べて滑らかさがあります。

（参考：ICARランダム効果モデル）

モデル式12-11：ARモデル

時間構造と空間構造には等価性があるとのことです。
モデル式12-11を時間構造で用いられる AR() で表現してみましょう。

PyMCのモデル定義

モデルの定義です。 
自己回帰係数 rho=1 の１次のARモデルです。

### モデルの定義 ◆モデル式12-11 model12-11.stan ※AR利用

with pm.Model() as model1_3:

    ### データ関連定義
    # coordの定義
    model1_3.add_coord('data', values=data.index, mutable=True)
    # dataの定義
    y = pm.ConstantData('y', value=data.Y.values, dims='data')

    ### 事前分布
    # 標準偏差
    sigma = pm.Uniform('sigma', lower=0, upper=100) 
    # AR(1)
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    r = pm.AR('r', rho=1, sigma=sigma, constant=False,
              init_dist=init_dist, ar_order=1, dims='data')

    ### 尤度
    Y = pm.Poisson('Y', mu=pt.exp(r), observed=y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model1_3

【実行結果】

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

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model1_3:
    idata1_3 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['sigma','r']
pm.summary(idata1_3, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata1_3, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

推論結果の解釈

推論データを用いてテキスト図12.7左に相当するグラフを描画します。

### rの推定結果の描画 ◆図12.7左
plot_pred(idata1_3)

【実行結果】
切片＋ICARランダム効果モデルとよく似たプロットになりました。
時間構造と空間構造の等価性があることも良く分かりました。

（参考：切片＋ICARランダム効果モデル）

（参考：ICARランダム効果モデル）

モデル式12-12：ARモデル

隣$${r[i-1]}$$と隣の隣$${r[i-2]}$$を用いて滑らかな空間構造を推論します。
こちらも AR() で実践します。
model12-11.stan の 12行目$${2r[i-1]-r[i-2]}$$に準拠して、自己回帰係数$${2,\ -1}$$を AR の rho 引数に与えます。

PyMCのモデル定義

モデルの定義です。 

### モデルの定義 ◆モデル式12-12 model12-12.stan ※AR利用

with pm.Model() as model2:

    ### データ関連定義
    # coordの定義
    model2.add_coord('data', values=data.index, mutable=True)
    # dataの定義
    y = pm.ConstantData('y', value=data.Y.values, dims='data')

    ### 事前分布
    # 標準偏差
    sigma = pm.Uniform('sigma', lower=0, upper=100) 
    # AR(2)
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    r = pm.AR('r', rho=[2, -1], sigma=sigma, constant=False,
              init_dist=init_dist, ar_order=2, dims='data')

    ### 尤度
    Y = pm.Poisson('Y', mu=pt.exp(r), observed=y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model2

【実行結果】

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

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model2:
    idata2 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['sigma', 'r']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata2, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】
発散を示す黒いバーコードが見られます。
いったん見なかったことにします。

推論結果の解釈

推論データを用いてテキスト図12.7右に相当するグラフを描画します。

### rの推定結果の描画 ◆図12.7右
plot_pred(idata2)

【実行結果】
テキストとよく似たプロットになりました。
１つ前のモデル式12-11のARモデルと比べて、かなり滑らかな線になりました。

（参考：モデル式12-11 ARモデル）

12.6 節は以上です。

シリーズの記事

次の記事

前の記事

目次


ブログの紹介

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

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

この記事が参加している募集

#新生活をたのしく

47,145件

#うちの積読を紹介する

1,724件

この記事が気に入ったらサポートをしてみませんか？