StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第４章「4.4 単回帰」

ネイピア DS

2024年2月4日 13:30

第４章「StanとRStanをはじめよう」

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

この記事は、テキスト第４章「StanとRStanを始めよう」の4.4節「単回帰」の PyMC5写経 を取り扱います。

はじめに

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

4.4 単回帰

インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np
import statsmodels.api as sm

# PyMC
import pymc as pm
import arviz as az
import bambi as bmb

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

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

4.4.2 データの分布の確認

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

### データの読み込み ◆データファイル4.1 data-salary.txt
# X:社員の年齢(歳), Y:社員の年収(万円)

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

【実行結果】

散布図を描画します。図4.2相当です。

### 散布図の描画 by seaborn ◆図4.2
sns.lmplot(data=data, x='X', y='Y', line_kws={'color': 'tomato'})
plt.title(f'相関係数: {data.corr().iloc[0, 1]:.3f}')
plt.grid(lw=0.5);

【実行結果】

4.4.4 Rのlm関数で推定

statsmodels.api の OLS で推定します。
4.4節で最も難しかったです。

### 回帰直線にあてはめ by statsmodels ◆lm.R 1~2行目

# Y, Xの取り出し、Xに定数項設定
Y = data.Y
X = sm.add_constant(data.X)
# 線形回帰モデルの設定
model_lm = sm.OLS(Y, X)
# 回帰分析の実行
res_lm = model_lm.fit()
# 回帰分析の結果表示（係数部分のみ）
display(res_lm.summary2().tables[1].round(2))

【実行結果】
const の Coef. が切片、X の Coef. が傾きです。

予測を実行します。

### 予測の実行(get_prediction) ◆lm.R 3~5行目

# 予測するXの値の設定（23歳～60歳）とXに定数項設定
X_new = list(range(23, 61))
X_new_const = sm.add_constant(X_new)
# 予測の実行
pred = res_lm.get_prediction(X_new_const)
# 予測結果をデータフレーム化（区間は95%）
pred_df = pred.summary_frame(alpha=0.05)
pred_df.index = X_new
display(pred_df.round(1))    # mean_ciは信頼区間、obs_ciは予測区間

【実行結果】

予測結果の描画を行います。
まず描画用データを整えます。

### 予測結果predより描画用データの取得 ◆lm.R 3~5行目

# 予測の平均値
pred_mean = pred.summary_frame()['mean']
# 予測の平均値の信頼区間 95%, 50%
pred_mean_95l = pred.summary_frame(alpha=0.05)['mean_ci_lower']
pred_mean_95u = pred.summary_frame(alpha=0.05)['mean_ci_upper']
pred_mean_50l = pred.summary_frame(alpha=0.50)['mean_ci_lower']
pred_mean_50u = pred.summary_frame(alpha=0.50)['mean_ci_upper']
# 予測（ノイズ込み）の予測区間 95%, 50%
pred_obs_95l = pred.summary_frame(alpha=0.05)['obs_ci_lower']
pred_obs_95u = pred.summary_frame(alpha=0.05)['obs_ci_upper']
pred_obs_50l = pred.summary_frame(alpha=0.50)['obs_ci_lower']
pred_obs_50u = pred.summary_frame(alpha=0.50)['obs_ci_upper']

基本年収の信頼区間を描画します。図4.3の左です。

# 基本年収の信頼区間の描画 ◆図4.3左

# axesの設定(ラベル、タイトルをset()１行で書きたい)
ax = plt.subplot()
# 観測値の散布図の描画
ax.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax.plot(X_new, pred_mean, color='tomato')
# 信頼区間の描画（95%、50%）
ax.fill_between(X_new, pred_mean_95l, pred_mean_95u, color='tomato', alpha=0.2)
ax.fill_between(X_new, pred_mean_50l, pred_mean_50u, color='tomato', alpha=0.4)
# 修飾
ax.set(xlabel='$X$', ylabel='$Y$', title='基本年収の信頼区間')
plt.grid(lw=0.5);

【実行結果】

年収の予測区間を描画します。図4.3の右です。

# 年収の予測区間の描画 ◆図4.3右

# axesの設定(ラベル、タイトルをset()１行で書きたい)
ax = plt.subplot()
# 観測値の散布図の描画
ax.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax.plot(X_new, pred_mean, color='tomato')
# 予測区間の描画（95%、50%）
ax.fill_between(X_new, pred_obs_95l, pred_obs_95u, color='tomato', alpha=0.2)
ax.fill_between(X_new, pred_obs_50l, pred_obs_50u, color='tomato', alpha=0.4)
# 修飾
ax.set(xlabel='$X$', ylabel='$Y$', title='年収の予測区間')
plt.grid(lw=0.5);

【実行結果】

4.4.5 Stanで実装

PyMC Ver.5 で実装します。
モデルの定義です。

### モデルの定義 ◆model4-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')

    ### 事前分布
    a = pm.Uniform('a', lower=-10000, upper=10000)
    b = pm.Uniform('b', lower=-10000, upper=10000)
    sigma = pm.Uniform('sigma', lower=0, upper=10000)
    mu = pm.Deterministic('mu', a + b * X, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

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

### モデルの表示
model1

【実行結果】

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

【実行結果】

4.4.6 Rから実行

PythonでMCMCを実行します。

### 事後分布からのサンプリング 5秒 ◆run-model4-5.R
# テキスト: iter=2000, warmup=1000, chain=4, thin=1
with model1:
    idata1 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

4.4.7 RStanの見方

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 = ['a', 'b', 'sigma']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=2)

【実行結果】

4.4.8 収束診断をファイルへ出力する

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

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

【実行結果】

【番外編：Bambiでベイズモデリング】

R の Formula 式ライクな記述でベイズモデリングできる PyMC のお友達 Bambi でモデルの定義～トレースプロットの描画を実践します。

### モデルの定義
model1_bmb = bmb.Model(
        formula='Y ~ 1 + X',   # フォーミュラ式
        data=data,             # データ
)

# モデルの表示
model1_bmb

【実行結果】

### モデルの可視化
model1_bmb.build()
model1_bmb.graph()

【実行結果】

### 事後分布からのサンプリング 5秒
# テキスト: iter=2000, warmup=1000, chain=4, thin=1
idata1_bmb = model1_bmb.fit(draws=1000, tune=1000, chains=4, target_accept=0.8,
                            nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

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

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

【実行結果】

### 推論データの要約統計情報の表示
pm.summary(idata1_bmb, hdi_prob=0.95, round_to=2)

【実行結果】

（参考：PyMCのmodel1の事後統計量）

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

【実行結果】

■ Bambiまとめ
事前分布やそのパラメータ値をBambiが自動的に設定します。
通常のPyMCとよく似た結果になったと思います。

以上で番外編を終了いたします。

4.4.9 MCMCの設定の変更

パラメータの初期値の設定、乱数シードの設定、chain・warmup・iterの設定に取り組みます。

### 事後分布からのサンプリング 5秒
# ◆rstan-modify-MCMCsettings.R
# テキスト: iter=1000, warmup=200, chain=3, thin=2

# 初期値に用いる乱数生成器の設定
rng = np.random.default_rng(seed=123)
# 初期値の設定 by 辞書
init_vals = {'a': rng.uniform(low=-10, high=10, size=1)[0],
             'b': rng.uniform(low=0, high=10, size=1)[0],
             'sigma': 10}

# MCMCの実行
with model1:
    idata1_mod = pm.sample(draws=800, tune=200, chains=3,  # draws=iter-warmup
                           initvals=init_vals,             # 初期値
                           random_seed=123,                # 乱数シード
                           target_accept=0.8,
                           nuts_sampler='numpyro')

【実行結果】省略
MCMC実行時にthinを指定できません。PyMCのポリシーらしいです。

上記の結果で収束確認をしましょう（収束していないです）。

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

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

【実行結果】

### 推論データの要約統計情報の表示
pm.summary(idata1_mod, hdi_prob=0.95, round_to=2)

【実行結果】

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

【実行結果】

■ thin の対応
model1の事後分布からのサンプリングデータを格納した idata1 に対して thin （間引き）を実行してみます。
xarray.Dataset の thin を利用して２つづつ間引きし、トレースプロットを描画します。

### thinの実行 xarrayのDataset型のthinメソッドを利用 ◆thinの代用
pm.plot_trace(az.extract(idata1).thin(2))
plt.tight_layout();

【実行結果】

4.4.10 並列計算の実行方法

pm.sample() の引数 cores で並列計算の設定（最大４）ができます。
引数を省略する場合、システムのCPU数が設定されます。

### 並列計算:coresで指定 40秒
with model1:
    idata1_core = pm.sample(cores=4)

【実行結果】

4.4.11 ベイズ信頼区間とベイズ予測区間の算出

事後分布からのサンプリングデータを収めた idata1 （arviz の InferenceData形式）を閲覧してみます。

### 事後分布からのサンプリングデータの表示
idata1

【実行結果】

パラメータ b （傾き）にアクセスします。
stack() では chain と draw の２次元を１次元に変換します。

### idataの変数bへアクセス ※InferenceDataのstackで次元を圧縮
idata1.posterior.b.stack(sample=('chain', 'draw')).data

【実行結果】

パラメータ b （傾き）の95%ベイズ信頼区間を算出します。

### idataの変数bの95%ベイズ信頼区間 ※valuesでnumpy配列に変換して平坦化
np.quantile(idata1.posterior.b.values.flatten(),
            [0.025, 0.975])

【実行結果】

別の方法を試してみます。
arviz の extract で chain と draw の２次元を１次元に変換します。

### arvizのextractでchainとdrawの次元を解除 ◆rstan-extract-MCMCsamples.R 1~3行目
idata1_ext = az.extract(idata1)
idata1_ext

【実行結果】

パラメータ b （傾き）にアクセスして95%ベイズ信頼区間を算出します。

### extract後のデータで変数bにアクセス
print('bのサンプル:\n', idata1_ext.b.data)
print('bの95%ベイズ信頼区間: ', np.quantile(idata1_ext.b.values, [0.025, 0.975]))

【実行結果】

テキストの d_mcmc （データフレーム）を作成します。

### a,b,σのデータフレーム化
d_mcmc = (idata1.posterior[['a', 'b', 'sigma']]
          .to_dataframe().reset_index(drop=True))
d_mcmc.head()

【実行結果】

図4.7の a, b の散布図を描画します。

### a, bの散布図の描画 ※図4.7
sns.jointplot(data=d_mcmc, x='a', y='b', kind='reg', marginal_ticks=True,
              scatter_kws={'alpha': 0.3, 'ec': 'white'},
              marginal_kws={'ec': 'white'})
plt.grid(lw=0.5)

【実行結果】

50歳の基本年収の事後分布と予測分布のデータフレームを作成します。

### 50歳の基本年収の事後分布と予測分布 ※rstan-extract-MCMCsamples.R 4~7行目

## 初期値設定
# 乱数生成器の設定
rng = np.random.default_rng(seed=1234)
# MCMCサンプリングのサイズ
N_mcmc = d_mcmc.shape[0]

## 計算
# 50歳の基本年収の事後分布の算出
y50_base = d_mcmc.a + d_mcmc.b * 50
# 50歳の年収の事後予測の算出
y50 = rng.normal(loc=y50_base, scale=d_mcmc.sigma, size=N_mcmc)
# データフレームに統合
d_mcmc2 = pd.concat([d_mcmc, pd.DataFrame({'y50_base': y50_base, 'y50': y50})],
                     axis=1)
display(d_mcmc2.head())

【実行結果】

テキストの図4.8の23～60歳における基本年収のベイズ信頼区間と年収のベイズ予測区間を描画します。

### 23～60歳における基本年収のベイズ信頼区間と年収のベイズ予測区間 ◆図4.8

### 初期値設定
# 乱数生成器の設定
rng = np.random.default_rng(seed=1234)

## 23~60歳の基本年収の事後分布と年収の予測分布の算出
posts, preds = [], []
# 年齢ごとに事後分布・予測分布の算出を繰り返し処理
for i in X_new:
    # 基本年収の事後分布の算出
    post_dist = d_mcmc.a + d_mcmc.b * i
    # 年収の予測分布の算出
    pred_dist = rng.normal(loc=post_dist, scale=d_mcmc.sigma, size=N_mcmc)
    # 算出結果をリストに追加
    posts.append(post_dist)
    preds.append(pred_dist)
# リストからnumpy配列に変換
posts = np.array(posts)
preds = np.array(preds)

### 中央値値・95%区間・50%区間の計算
# 基本年収の事後分布
posts_means = np.median(posts, axis=1)
posts_95ci = np.quantile(posts, [0.025, 0.975], axis=1)
posts_50ci = np.quantile(posts, [0.250, 0.750], axis=1)
# 年収の予測分布
preds_means = np.median(preds, axis=1)
preds_95ci = np.quantile(preds, [0.025, 0.975], axis=1)
preds_50ci = np.quantile(preds, [0.250, 0.750], axis=1)

### 描画
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
## 基本年収のベイズ信頼区間の描画
# 観測値の散布図の描画
ax1.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax1.plot(X_new, posts_means, color='tomato')
# 予測区間の描画（95%、50%）
ax1.fill_between(X_new, posts_95ci[0], posts_95ci[1], color='tomato', alpha=0.2)
ax1.fill_between(X_new, posts_50ci[0], posts_50ci[1], color='tomato', alpha=0.4)
# 修飾
ax1.set(xlabel='$X$', ylabel='$Y$', title='基本年収のベイズ信頼区間')
ax1.grid(lw=0.5);
## 年収の予測区間の描画
# 観測値の散布図の描画
ax2.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax2.plot(X_new, preds_means, color='tomato')
# 予測区間の描画（95%、50%）
ax2.fill_between(X_new, preds_95ci[0], preds_95ci[1], color='tomato', alpha=0.2)
ax2.fill_between(X_new, preds_50ci[0], preds_50ci[1], color='tomato', alpha=0.4)
# 修飾
ax2.set(xlabel='$X$', ylabel='$Y$', title='年収のベイズ予測区間')
ax2.grid(lw=0.5);

【実行結果】

4.4.12 transformed parametersブロックとgenerated quantitiesブロック

Y の事後分布と予測分布をMCMCで算出します。
まずモデルの定義から。

### モデルの定義 ◆model4-4.stan

with pm.Model() as model1_gen:
    
    ### データ関連定義
    ## coordの定義
    model1_gen.add_coord('data', values=data.index, mutable=True)
    model1_gen.add_coord('dataNew', values=range(len(X_new)), 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')
    # 予測用の説明変数 x_new
    XNew = pm.ConstantData('XNew', value=X_new, dims='dataNew')

    ### 事前分布
    a = pm.Uniform('a', lower=-10000, upper=10000)
    b = pm.Uniform('b', lower=-10000, upper=10000)
    sigma = pm.Uniform('sigma', lower=0, upper=10000)
    mu = pm.Deterministic('mu', a + b * X, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

    ### 計算
    yBaseNew = pm.Deterministic('yBaseNew', a + b * XNew, dims='dataNew')
    yNew = pm.Normal('yNew', mu=yBaseNew, sigma=sigma, dims='dataNew')

モデルを表示します。

### モデルの表示
model1_gen

【実行結果】

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

【実行結果】

MCMCを実行します。

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

【実行結果】省略

収束確認します（収束してます）。

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

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

【実行結果】

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

【実行結果】

### トレースプロットの表示
var_names = ['a', 'b', 'sigma', 'mu', 'yBaseNew', 'yNew']
pm.plot_trace(idata1_gen, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

MCMCでサンプリングした23～60歳における基本年収の事後分布と年収の予測分布を描画します。

### 23～60歳における基本年収のベイズ信頼区間と年収のベイズ予測区間 ◆run-model4-4.R

### 描画用データ加工
# 2.5%,25%,50%,75%,97.5%パーセンタイル点を算出してデータフレーム化する関数の定義
def make_quantile_df(x, y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = [str(s) + '%' for s in probs]
    quantiles = np.percentile(y, probs, axis=1)
    return pd.DataFrame(quantiles.T, columns=columns, index=x)

# 推論データidata_genからyの事後分布yBasenewと予測分布yNewを取り出し
y_base_news = idata1_gen.posterior.yBaseNew.stack(sample=('chain', 'draw')).data
y_news = idata1_gen.posterior.yNew.stack(sample=('chain', 'draw')).data

# 基本年収の事後分布のパーセンタイル点の算出（関数利用）
y_base_df = make_quantile_df(X_new, y_base_news)
print('基本年収の事後分布')
display(y_base_df.head().round(1))
# 年収の予測分布のパーセンタイル点の算出（関数利用）
y_pred_df = make_quantile_df(X_new, y_news)
print('年収の予測分布')
display(y_pred_df.head().round(1))

### 描画
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
## 基本年収のベイズ信頼区間の描画
# 観測値の散布図の描画
ax1.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax1.plot(y_base_df.index, y_base_df['50%'], color='tomato')
# 予測区間の描画（95%、50%）
ax1.fill_between(y_base_df.index, y_base_df['2.5%'], y_base_df['97.5%'],
                  color='tomato', alpha=0.2)
ax1.fill_between(y_base_df.index, y_base_df['25%'], y_base_df['75%'],
                 color='tomato', alpha=0.4)
# 修飾
ax1.set(xlabel='$X$', ylabel='$Y$', title='基本年収のベイズ信頼区間')
ax1.grid(lw=0.5);
## 年収の予測区間の描画
# 観測値の散布図の描画
ax2.scatter(data['X'], data['Y'], s=100, alpha=0.5)
# 予測の平均値の描画
ax2.plot(y_pred_df.index, y_pred_df['50%'], color='tomato')
# 予測区間の描画（95%、50%）
ax2.fill_between(y_pred_df.index, y_pred_df['2.5%'], y_pred_df['97.5%'],
                 color='tomato', alpha=0.2)
ax2.fill_between(y_pred_df.index, y_pred_df['25%'], y_pred_df['75%'],
                 color='tomato', alpha=0.4)
# 修飾
ax2.set(xlabel='$X$', ylabel='$Y$', title='年収のベイズ予測区間')
ax2.grid(lw=0.5);

【実行結果】

4.4 節は以上です。

第４章末の練習問題

### 練習問題で用いるデータの作成
rng = np.random.default_rng(seed=123)
N1 = 30
N2 = 20
Y1 = rng.normal(loc=0, scale=5, size=N1) # 平均0、標準偏差5の正規分布乱数
Y2 = rng.normal(loc=1, scale=4, size=N2) # 平均1、標準偏差4の正規分布乱数

問題(1) ヒストグラムバージョン

### 問題(1)各グループの差をおおよそ把握するための図の描画 ヒストグラム
sns.histplot(Y1, bins=10, stat='density', kde=True, ec='white', alpha=0.3,
             label='Y1')
sns.histplot(Y2, bins=10, stat='density', kde=True, ec='white', alpha=0.3,
             label='Y2')
plt.legend();

【実行結果】

問題(1) 箱ひげ図バージョン

### 問題(1)各グループの差をおおよそ把握するための図の描画 箱ひげ図
sns.boxplot([Y1, Y2], fill=False)
sns.swarmplot([Y1, Y2], s=10, alpha=0.5)
plt.xticks([0, 1], ['Y1', 'Y2'])
plt.grid(lw=0.5);

【実行結果】

問題(2)(3) モデルの作成～MCMC実行～収束確認

### 問題(2)(3)各グループで標準偏差が等しいと仮定してモデルを作成～MCMC実行

with pm.Model() as model_exe1:
    
    ### データ関連定義
    ## coordの定義
    model_exe1.add_coord('data1', values=range(len(Y1)), mutable=True)
    model_exe1.add_coord('data2', values=range(len(Y2)), mutable=True)
    ## dataの定義
    y1 = pm.ConstantData('y1', value=Y1, dims='data1')
    y2 = pm.ConstantData('y2', value=Y2, dims='data2')

    ### 事前分布
    mu1 = pm.Uniform('mu1', lower=-100, upper=100)
    mu2 = pm.Uniform('mu2', lower=-100, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    
    ### 尤度関数
    obs1 = pm.Normal('obs1', mu=mu1, sigma=sigma, observed=y1, dims='data1')
    obs2 = pm.Normal('obs2', mu=mu2, sigma=sigma, observed=y2, dims='data2')

### モデルの表示
model_exe1

【実行結果】

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

【実行結果】

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

【実行結果】省略

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

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

【実行結果】

### 推論データの要約統計情報の表示
pm.summary(idata_exe1, hdi_prob=0.95, round_to=2)

【実行結果】

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

【実行結果】

問題(4)

### 練習問題(4) Prob[μ1<μ2]を計算する
mu1_samples1 = idata_exe1.posterior.mu1.stack(sample=('chain', 'draw')).data
mu2_samples1 = idata_exe1.posterior.mu2.stack(sample=('chain', 'draw')).data
print(f'μ1 < μ2の確率: {(mu1_samples1 < mu2_samples1).mean():.3f}')

【実行結果】
練習問題の解答（サンプルコード）と値が全然違います・・・。

問題(5) モデルの作成～MCMC実行～収束確認～確率計算

### 練習問題(5) グループごとの標準偏差が異なる場合のモデルを作成～MCMC実行

with pm.Model() as model_exe2:
    
    ### データ関連定義
    ## coordの定義
    model_exe2.add_coord('data1', values=range(len(Y1)), mutable=True)
    model_exe2.add_coord('data2', values=range(len(Y2)), mutable=True)
    ## dataの定義
    y1 = pm.ConstantData('y1', value=Y1, dims='data1')
    y2 = pm.ConstantData('y2', value=Y2, dims='data2')

    ### 事前分布
    mu1 = pm.Uniform('mu1', lower=-100, upper=100)
    mu2 = pm.Uniform('mu2', lower=-100, upper=100)
    sigma1 = pm.Uniform('sigma1', lower=0, upper=100)
    sigma2 = pm.Uniform('sigma2', lower=0, upper=100)
    
    ### 尤度関数
    obs1 = pm.Normal('obs1', mu=mu1, sigma=sigma1, observed=y1, dims='data1')
    obs2 = pm.Normal('obs2', mu=mu2, sigma=sigma2, observed=y2, dims='data2')

### モデルの表示
model_exe2

【実行結果】

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

【実行結果】

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

【実行結果】省略

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

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

【実行結果】

### 推論データの要約統計情報の表示
pm.summary(idata_exe2, hdi_prob=0.95, round_to=2)

【実行結果】

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

【実行結果】

確率の計算です。
練習問題の解答（サンプルコード）と値が全然違います・・・。

### Prob[μ1<μ2]を計算する
mu1_samples2 = idata_exe2.posterior.mu1.stack(sample=('chain', 'draw')).data
mu2_samples2 = idata_exe2.posterior.mu2.stack(sample=('chain', 'draw')).data
print(f'μ1 < μ2の確率: {(mu1_samples2 < mu2_samples2).mean():.3f}')

【実行結果】

以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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