StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第９章「9.3.1 多変量正規分布」

2024年3月11日 11:26

第９章「一歩進んだ文法」

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

この記事は、テキスト第９章「一歩進んだ文法」・9.3節「ベクトルや行列の数学的性質の利用」の 9.3.1項「多変量正規分布」の PyMC5写経 を取り扱います。

テキストは第９章で Stan の文法上の工夫を取り扱っています。
Stanの文法とPyMCの文法は異なっており、Stan の工夫点をPyMCに取り入れることができない、または、取り入れる方法が分からない場合があります。
したがいまして第９章では、私個人のスキルと相談して、PyMC化に意味を見いだせるテーマを選択して取り上げています。
選択の結果、9.1「型とインデックス」および 9.2「ベクトルによる高速化」の写経は省略しました。

はじめに

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を動かすまでの準備」章をご覧ください。

9.3.1 多変量正規分布

テキストは多変量正規分布に関する２点の Stan の工夫を取り扱っています。
①分散共分散行列：Stan の「cov_matrix型」の利用
②多変量正規分布：Stan の「multi_normal関数」のベクトル化

２つの工夫点に関するPyMCモデリングの対応内容は以下のとおりです。
①の分散共分散行列：LKJCholeskyCov() を利用
　なお10.2.4項の写経にて、分散共分散行列の複数案に取り組みました。
②のベクトル化：PyMC化の方策が分からないので未実施

インポート

### インポート

# 数値・確率計算
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
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

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

データの読み込み・確認

サンプルコードのデータを読み込みます。

### データの読み込み ◆データファイル9.1 data-mvn.txt
# Y1:50m走のタイム[秒], Y2:走り幅跳びで飛んだ距離[m]

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

【実行結果】

データの外観を確認しましょう。
まずは要約統計量から。

### 要約統計量の表示
data1.describe().round(1)

【実行結果】

続いて50m走のタイムと走り幅跳びの距離の相関係数を確認します。

### 相関係数の算出
data1.corr().round(3)

【実行結果】
強い負の相関＝タイムが短いほど距離が長くなる関係のようです。

相関係数はベイズ推論結果で確かめます。
$${-0.859}$$の値を覚えておきます。

散布図で可視化して相関関係を見てみましょう。
テキスト図9.1に相当します。

### 散布図の描画 ◆図9.1
plt.figure(figsize=(5, 4))
ax = plt.subplot()
sns.scatterplot(data=data1, x='Y1', y='Y2', s=100, alpha=0.5, ax=ax)
ax.set(xlabel='$Y1$: 50m走のタイム[秒]', ylabel='$Y2$: 走り幅跳びの跳躍距離[m]',
       yticks=(2.0, 2.5, 3.0, 3.5))
plt.grid(lw=0.5);

【実行結果】
右肩下がりの傾向がはっきり見られます。

PyMCのモデル定義

PyMCでモデル式9-2を実装します。
モデルの定義です。
分散共分散の表現に LKJCholeskyCov() を用いています。

### モデルの定義 ◆モデル式9-2 model9-2.stan

with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data1.index, mutable=True)
    model1.add_coord('param', values=['Y1', 'Y2'], mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data1[['Y1', 'Y2']].values,
                        dims=('data', 'param'))

    ### 事前分布
    # 平均ベクトル（要素数2）
    mu = pm.Uniform('mu', lower=-100, upper=100, dims='param')    
    # コレスキー共分散分解（2x2行列）
    sd_dist = pm.Uniform.dist(lower=0, upper=100, size=2)
    chol, _, _ = pm.LKJCholeskyCov('chol_cov', eta=1, n=2, sd_dist=sd_dist)

    ### 尤度関数
    obs = pm.MvNormal('obs', mu=mu, chol=chol, observed=Y,
                      dims=('data', 'param'))
    
    ### 計算値 分散共分散行列
    cov = pm.Deterministic('cov', pt.dot(chol, chol.T))

モデルの定義内容を見ます。
LKJCholeskyCov() を用いた関係で変数名にアンダースコアが入り、KaTeX エラーが発生するので、別の方法で確率変数の確率分布を確認します。

### モデルの表示
model1.basic_RVs

【実行結果】

### モデルの可視化 2分10秒
pm.model_to_graphviz(model1)

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 15秒
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.01         # しきい値

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

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

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

### 推論データの要約統計情報の表示
var_names = ['mu', 'cov', 'chol_cov', 'chol_cov_corr', 'chol_cov_stds']
pm.summary(idata1, hdi_prob=0.95, round_to=3, var_names=var_names)

【実行結果】

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

### トレースプロットの表示
var_names = ['mu', 'cov', 'chol_cov', 'chol_cov_stds']
pm.plot_trace(idata1, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

推定結果の解釈

事後分布の要約統計量を算出します。
算出関数を定義します。

### mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

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

### 要約統計量の算出・表示 ◆テキスト１５９ページ パラメータの事後分布
mu_samples = pd.DataFrame(
            idata1.posterior.mu.stack(sample=('chain', 'draw')).data.T,
            columns=[f'mu[{i+1}]' for i in range(2)])
cov1_samples = pd.DataFrame(
            idata1.posterior.cov.stack(sample=('chain', 'draw')).data[0].T,
            columns=[f'cov[1, {i+1}]' for i in range(2)]
            )
cov2_samples = pd.DataFrame(
            idata1.posterior.cov.stack(sample=('chain', 'draw')).data[1].T,
            columns=[f'cov[2, {i+1}]' for i in range(2)]
            )
param_samples = pd.concat([mu_samples, cov1_samples, cov2_samples], axis=1)
stats_df = make_stats_df(param_samples)
display(stats_df.round(2))

【実行結果】
平均 mu の平均値はテキストと同じ結果になりました。
共分散 cov の平均値はテキストと少し異なる結果になりました。

共分散 cov を用いて、Y1（50m走タイム）とY2（走り幅跳び距離）の相関を計算します。

### Y1とY2の平均的な相関 ◆テキスト１５９ページ
corr = (stats_df.iloc[3, 0]
        / (np.sqrt(stats_df.iloc[2, 0]) * np.sqrt(stats_df.iloc[5, 0])))
print(f'{corr:.2f}')

【実行結果】
テキストの$${-0.86}$$に近い結果$${-0.84}$$を得られました。
データから算出した相関係数$${-0.859}$$と比べて若干ずれています（残念）。

9.3.1 項は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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

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

いいなと思ったら応援しよう！

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

#仕事について話そう

115,662件