見出し画像

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第10章「10.2.4 分散共分散行列」

ネイピア DS

第10章「収束しない場合の対処法」

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

この記事は、テキスト第10章「収束しない場合の対処法」・10.2節「弱情報事前分布」の 10.2.4項「分散共分散行列」の PyMC5写経 を取り扱います。

テキストは第10章で 収束に向けた工夫を取り扱っています。
座学や数式モデルのみの掲載項があれば、Stanコードの掲載項もあります。
PyMC化は主にStanコードが明示されているモデル式を対象にして実施します。
したがいまして、10.2.3「確率など[0, 1]の範囲のパラメータ」の写経は省略しました。

はじめに

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


10.2.4 分散共分散行列

インポート

### インポート

# 数値・確率計算
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')

データの読み込み・確認

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

### データの読み込み ◆データファイル8.1 data-salary-2.txt
# X:年齢-23, Y:年収, KID:勤務会社ID

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

【実行結果】

モデル式10-5 PyMCのモデル定義

PyMCでモデル式10-5を実装します。
データの前処理を行います。

### 会社インデックスの作成（0始まりにする）
k_idx = data.KID - 1

【実行結果】なし

モデルの定義です。
分散 var と共分散 cov に一様分布を用いて、分散共分散行列 covMtx を表現します。

### モデルの定義 ◆モデル式10-5 model10-5.stan

with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data.index, mutable=True)

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス
    KIdx = pm.ConstantData('KIdx', value=k_idx.values, dims='data')

    ### 事前分布
    # 会社別係数の２変量正規分布の平均 ab0
    ab0 = pm.Uniform('ab0', lower=-10000, upper=10000, shape=2)
    # 会社別係数の２変量正規分布の分散共分散行列 covMtx
    var = pm.Uniform('var', lower=0, upper=10000, shape=2)
    cov = pm.Uniform('cov', lower=-10000, upper=10000)
    covMtx = pm.Deterministic('covMtx',
                              pt.stacklists([[var[0], cov],[cov, var[1]]]))
    # 会社別係数：２変量正規分布
    ab = pm.MvNormal('ab', mu=ab0, cov=covMtx, shape=(4, 2))
    # 尤度の正規分布の標準偏差
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=10000)

    ### 線形予測子
    mu = pm.Deterministic('mu', ab[KIdx, 0] + ab[KIdx, 1] * X, dims='data')

    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigmaY, observed=Y, dims='data')

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

### モデルの表示
model1

【実行結果】

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

【実行結果】

モデル式10-5 MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 2分20秒
with model1:
    idata1  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
                        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 = ['ab0', 'covMtx', 'ab', 'sigmaY', 'mu']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

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

### トレースプロットの表示
var_names = ['ab0', 'covMtx', 'ab', 'sigmaY', 'mu']
pm.plot_trace(idata1, compact=True, divergences=None, var_names=var_names)
plt.tight_layout();

【実行結果】

モデル式10-5 推定結果の解釈

平均パラメータ ab0、分散共分散行列 covMtx、尤度関数の標準偏差 sigmaY について、事後分布の要約統計量を算出します。
算出関数を定義します。

### 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

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

### 要約統計量の算出・表示
vars = ['sigmaY']
param_samples = idata1.posterior[vars].to_dataframe().reset_index(drop=True)
ab0_samples = pd.DataFrame(
            idata1.posterior.ab0.stack(sample=('chain', 'draw')).data.T,
            columns=[f'ab0[{i+1}]' for i in range(2)])
covMtx1_samples = pd.DataFrame(
            idata1.posterior.covMtx.stack(sample=('chain', 'draw')).data[0].T,
            columns=[f'covMtx[1,{i+1}]' for i in range(2)])
covMtx2_samples = pd.DataFrame(
            idata1.posterior.covMtx.stack(sample=('chain', 'draw')).data[1].T,
            columns=[f'covMtx[2,{i+1}]' for i in range(2)])
param_samples = pd.concat([ab0_samples, covMtx1_samples,
                           covMtx2_samples, param_samples], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】
テキストに事後分布の推定値が掲載されていないので、PyMCモデルによる推論の適否は不明です。

モデル式10-6 PyMCのモデル定義

PyMCでモデル式10-6を実装します。
モデルの定義です。
標準偏差 sigma に半ｔ分布を、相関係数 rho に一様分布をそれぞれ用いて、分散共分散行列 covMtx を表現します。

### モデルの定義 ◆モデル式10-6 model10-6.stan

with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data.index, mutable=True)

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス
    KIdx = pm.ConstantData('KIdx', value=k_idx.values, dims='data')

    ### 事前分布
    # 会社別係数の２変量正規分布の平均 ab0
    ab0A = pm.Normal('ab0A', mu=400, sigma=200)
    ab0B = pm.Normal('ab0B', mu=15, sigma=15)
    ab0 = pm.Deterministic('ab0', pt.stack([ab0A, ab0B]))
    # 会社別係数の２変量正規分布の分散共分散行列 covMtx
    sigmaA = pm.HalfStudentT('sigmaA', nu=4, sigma=200)
    sigmaB = pm.HalfStudentT('sigmaB', nu=4, sigma=20)
    rho = pm.Uniform('rho', lower=-1, upper=1)
    covMtx = pm.Deterministic('covMtx',
                              pt.stacklists(
                                  [[pt.pow(sigmaA, 2), sigmaA * sigmaB * rho],
                                   [sigmaA * sigmaB * rho, pt.pow(sigmaB, 2)]]))
    # 会社別係数：２変量正規分布
    ab = pm.MvNormal('ab', mu=ab0, cov=covMtx, shape=(4, 2))
    # 尤度の正規分布の標準偏差
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=10000)

    ### 線形予測子
    mu = pm.Deterministic('mu', ab[KIdx, 0] + ab[KIdx, 1] * X, dims='data')

    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigmaY, observed=Y, dims='data')

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

### モデルの表示
model2

【実行結果】

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

【実行結果】

モデル式10-6 MCMCの実行と収束確認

MCMCを実行します。

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

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

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

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

### 推論データの要約統計情報の表示
var_names = ['sigmaA', 'sigmaB', 'rho', 'ab0', 'ab', 'sigmaY', 'mu']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

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

### トレースプロットの表示
var_names = ['sigmaA', 'sigmaB', 'rho', 'ab0', 'ab', 'sigmaY', 'mu']
pm.plot_trace(idata2, compact=True, divergences=None, var_names=var_names)
plt.tight_layout();

【実行結果】

モデル式10-6 推定結果の解釈

平均パラメータ ab0、分散共分散行列 covMtx、尤度関数の標準偏差 sigmaY について、事後分布の要約統計量を算出します。

### 要約統計量の算出・表示
vars = ['sigmaY']
param_samples = idata2.posterior[vars].to_dataframe().reset_index(drop=True)
ab0_samples = pd.DataFrame(
            idata2.posterior.ab0.stack(sample=('chain', 'draw')).data.T,
            columns=[f'ab0[{i+1}]' for i in range(2)])
covMtx1_samples = pd.DataFrame(
            idata2.posterior.covMtx.stack(sample=('chain', 'draw')).data[0].T,
            columns=[f'covMtx[1,{i+1}]' for i in range(2)])
covMtx2_samples = pd.DataFrame(
            idata2.posterior.covMtx.stack(sample=('chain', 'draw')).data[1].T,
            columns=[f'covMtx[2,{i+1}]' for i in range(2)])
param_samples = pd.concat([ab0_samples, covMtx1_samples,
                           covMtx2_samples, param_samples], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】
テキストに事後分布の推定値が掲載されていないので、PyMCモデルによる推論の適否は不明です。

（参考：モデル式10-5の事後分布の推定値）

モデル式10-7 PyMCのモデル定義

PyMCでモデル式10-7を実装します。
モデルの定義です。
LKJ 分散相関を使用した分散共分散行列 LKJCholeskyCov() を利用します。

### モデルの定義 ◆モデル式10-7 model10-7.stan

## LKJCholeskyCovのパラメータηの値の設定
Nu = 2

## モデルの定義
with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    model3.add_coord('data', values=data.index, mutable=True)

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス
    KIdx = pm.ConstantData('KIdx', value=k_idx.values, dims='data')

    ### 事前分布
    # 会社別係数の２変量正規分布の平均 ab0
    ab0A = pm.Normal('ab0A', mu=400, sigma=200)
    ab0B = pm.Normal('ab0B', mu=15, sigma=15)
    ab0 = pm.Deterministic('ab0', pt.stack([ab0A, ab0B]))
    # 会社別係数の２変量正規分布の分散共分散行列の代わりに
    # LKJ分布のコレスキー共分散分解 chol を利用
    # 異なる標準偏差の確率分布の情報: 
    # https://discourse.pymc.io/t/lkjcholeskycov-with-different-standard-deviation-distributions/9617
    sd_dist = pm.HalfStudentT.dist(nu=[4, 4], sigma=[200, 20]) # ２つの半ｔ分布
    chol, _, _ = pm.LKJCholeskyCov(
                'chol_cov', eta=Nu, n=2, compute_corr=True, sd_dist=sd_dist)
    # 会社別係数：２変量正規分布
    ab = pm.MvNormal('ab', mu=ab0, chol=chol, shape=(4, 2))
    # 尤度の正規分布の標準偏差
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=10000)

    ### 線形予測子
    mu = pm.Deterministic('mu', ab[KIdx, 0] + ab[KIdx, 1] * X, dims='data')

    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigmaY, observed=Y, dims='data')

    ### 計算値
    covMtx = pm.Deterministic('covMtx', chol @ chol.T)

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

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

【実行結果】

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

【実行結果】

モデル式10-7 MCMCの実行と収束確認

MCMCを実行します。

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

【実行結果】省略

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

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

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

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

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

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

【実行結果】

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

### トレースプロットの表示
var_names = ['ab0', 'chol_cov', 'ab', 'sigmaY', 'mu']
pm.plot_trace(idata3, compact=True, divergences=None, var_names=var_names)
plt.tight_layout();

【実行結果】

モデル式10-7 推定結果の解釈

平均パラメータ ab0、分散共分散行列 covMtx、尤度関数の標準偏差 sigmaY について、事後分布の要約統計量を算出します。

### 要約統計量の算出・表示
vars = ['sigmaY']
param_samples = idata3.posterior[vars].to_dataframe().reset_index(drop=True)
ab0_samples = pd.DataFrame(
            idata3.posterior.ab0.stack(sample=('chain', 'draw')).data.T,
            columns=[f'ab0[{i+1}]' for i in range(2)])
covMtx1_samples = pd.DataFrame(
            idata3.posterior.covMtx.stack(sample=('chain', 'draw')).data[0].T,
            columns=[f'covMtx[1,{i+1}]' for i in range(2)])
covMtx2_samples = pd.DataFrame(
            idata3.posterior.covMtx.stack(sample=('chain', 'draw')).data[1].T,
            columns=[f'covMtx[2,{i+1}]' for i in range(2)])
param_samples = pd.concat([ab0_samples, covMtx1_samples,
                           covMtx2_samples, param_samples], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】
テキストに事後分布の推定値が掲載されていないので、PyMCモデルによる推論の適否は不明です。

（参考：モデル式10-6の事後分布の推定値）

（参考：モデル式10-5の事後分布の推定値）

10.2.4 項は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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,013件