StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第８章「練習問題」

ネイピア DS

2024年3月10日 10:53

第８章「階層モデル」

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

この記事は、テキスト第８章「階層モデル」の「練習問題」の 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を動かすまでの準備」章をご覧ください。

8章練習問題

インポート

### インポート

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

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

練習問題(1)

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

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

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

モデル式8-2の再掲です。

### モデルの定義 ◆モデル式8-2 model8-2.stan

with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data1.index, mutable=True)
    model1.add_coord('kaisha', values=sorted(data1.KID.unique()), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data1['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス 0始まりに変換
    KIdx = pm.ConstantData('KIdx', value=data1['KID'].values - 1, dims='data')

    ### 事前分布
    a = pm.Uniform('a', lower=-10000, upper=10000, dims='kaisha')
    b = pm.Uniform('b', lower=-10000, upper=10000, dims='kaisha')
    sigma = pm.Uniform('sigma', lower=0, upper=10000)

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

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

MCMCを実行します。

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

【実行結果】省略

《解答》
出題のベイズ予測区間を描画します。

### 会社ごとの年収のベイズ予測区間の描画 ◆解答案

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
a_samples = idata1.posterior.a.stack(sample=('chain', 'draw')).data
b_samples = idata1.posterior.b.stack(sample=('chain', 'draw')).data
sigma_samples = idata1.posterior.sigma.stack(sample=('chain', 'draw')).data

# x軸の値の設定
x_vals = np.linspace(data1['X'].min(), data1['X'].max(), 101)

## 描画設定
# マーカーの設定
markers = {1: 'o', 2: '^', 3: 'X', 4: 'd'}
# カラーマップの取得
cmap = plt.get_cmap('tab10')

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(2, 2, figsize=(8, 6), sharex=True, sharey=True)
# 会社の数だけ、年収の予測分布を算出＆描画を繰り返し処理
for i in range(4):
    ## 描画用データの作成
    # 年収の予測分布を算出 scipy.statsで正規分布乱数を生成
    y_preds = np.array([stats.norm.rvs(loc=a + b*x_vals, scale=s)
                        for (a, b, s)
                        in zip(a_samples[i], b_samples[i], sigma_samples)]).T
    # 予測分布の95%区間、50%区間を算出
    y_preds_95ci = np.quantile(y_preds, q=[0.025, 0.975], axis=1)
    y_preds_50ci = np.quantile(y_preds, q=[0.250, 0.750], axis=1)
    
    ## 描画設定
    # axesの位置を算出
    pos = divmod(i, 2)
    # 色を算出
    color = cmap(i / 10)
    
    ## 描画
    # 予測分布の中央値の折れ線グラフの描画
    ax[pos].plot(x_vals, np.median(y_preds, axis=1), color='tab:red')
    # 予測分布の95%CIの塗りつぶし描画
    ax[pos].fill_between(x_vals, y_preds_95ci[0], y_preds_95ci[1],
                         color='tomato', alpha=0.1)
    # 予測分布の50%CIの塗りつぶし描画
    ax[pos].fill_between(x_vals, y_preds_50ci[0], y_preds_50ci[1],
                         color='tomato', alpha=0.3)
    # 観測値の散布図の描画
    ax[pos].scatter(data=data1[data1['KID']==i+1], x='X', y='Y',
                    marker=markers[i+1], s=100, color=color, alpha=0.8, label=i+1)
    # 修飾
    ax[pos].set(title=f'会社ID: {i+1}')
    ax[pos].grid(lw=0.5)

# 全体修飾
fig.supxlabel('年齢 $X$ [-23歳]', fontsize=16)
fig.supylabel('年収 $Y$ [万円]', fontsize=16)
fig.suptitle('年収 $Y$ の予測分布: 会社別', fontsize=16)
fig.legend(bbox_to_anchor=(1.1, 0.9), title='会社ID')
plt.tight_layout();

【実行結果】

練習問題(2)

練習問題(1)のデータを用いて、モデル式8-3のモデルを再掲します。

### モデルの定義 ◆モデル式8-3 model8-3.stan

with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data1.index, mutable=True)
    model2.add_coord('kaisha', values=sorted(data1.KID.unique()), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data1['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス 0始まりに変換
    KIdx = pm.ConstantData('KIdx', value=data1['KID'].values - 1, dims='data')

    ### 事前分布
    a0 = pm.Uniform('a0', lower=-10000, upper=10000)
    b0 = pm.Uniform('b0', lower=-10000, upper=10000)
    sigmaA = pm.Uniform('sigmaA', lower=0, upper=1000)
    sigmaB = pm.Uniform('sigmaB', lower=0, upper=1000)
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=10000)
    ak = pm.Normal('ak', mu=0, sigma=sigmaA, dims='kaisha')
    bk = pm.Normal('bk', mu=0, sigma=sigmaB, dims='kaisha')

    ### 線形予測子
    a = pm.Deterministic('a', a0 + ak, dims='kaisha')
    b = pm.Deterministic('b', b0 + bk, dims='kaisha')
    mu = pm.Deterministic('mu', a[KIdx] + b[KIdx] * X, dims='data')

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

MCMCを実行します。

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

【実行結果】省略

《解答》
出題のベイズ予測区間を描画します。

### 会社ごとの年収のベイズ予測区間の描画 ◆解答案

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
a_samples = idata2.posterior.a.stack(sample=('chain', 'draw')).data
b_samples = idata2.posterior.b.stack(sample=('chain', 'draw')).data
sigmaY_samples = idata2.posterior.sigmaY.stack(sample=('chain', 'draw')).data

# x軸の値の設定
x_vals = np.linspace(data1['X'].min(), data1['X'].max(), 101)

## 描画設定
# マーカーの設定
markers = {1: 'o', 2: '^', 3: 'X', 4: 'd'}
# カラーマップの取得
cmap = plt.get_cmap('tab10')

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(2, 2, figsize=(8, 6), sharex=True, sharey=True)
# 会社の数だけ、年収の予測分布を算出＆描画を繰り返し処理
for i in range(4):
    ## 描画用データの作成
    # 年収の予測分布を算出 scipy.statsで正規分布乱数を生成
    y_preds = np.array([stats.norm.rvs(loc=a + b*x_vals, scale=s)
                        for (a, b, s)
                        in zip(a_samples[i], b_samples[i], sigmaY_samples)]).T
    # 予測分布の95%区間、50%区間を算出
    y_preds_95ci = np.quantile(y_preds, q=[0.025, 0.975], axis=1)
    y_preds_50ci = np.quantile(y_preds, q=[0.250, 0.750], axis=1)
    
    ## 描画設定
    # axesの位置を算出
    pos = divmod(i, 2)
    # 色を算出
    color = cmap(i / 10)
    
    ## 描画
    # 予測分布の中央値の折れ線グラフの描画
    ax[pos].plot(x_vals, np.median(y_preds, axis=1), color='tab:red')
    # 予測分布の95%CIの塗りつぶし描画
    ax[pos].fill_between(x_vals, y_preds_95ci[0], y_preds_95ci[1],
                         color='tomato', alpha=0.1)
    # 予測分布の50%CIの塗りつぶし描画
    ax[pos].fill_between(x_vals, y_preds_50ci[0], y_preds_50ci[1],
                         color='tomato', alpha=0.3)
    # 観測値の散布図の描画
    ax[pos].scatter(data=data1[data1['KID']==i+1], x='X', y='Y',
                    marker=markers[i+1], s=100, color=color, alpha=0.8, label=i+1)
    # 修飾
    ax[pos].set(title=f'会社ID: {i+1}')
    ax[pos].grid(lw=0.5)

# 全体修飾
fig.supxlabel('年齢 $X$ [-23歳]', fontsize=16)
fig.supylabel('年収 $Y$ [万円]', fontsize=16)
fig.suptitle('年収 $Y$ の予測分布: 会社別', fontsize=16)
fig.legend(bbox_to_anchor=(1.1, 0.9), title='会社ID')
plt.tight_layout();

【実行結果】

練習問題(3)

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

### データの読み込み ◆データファイル8.2 data-salary-3.txt
# X:年齢-23, Y:年収, KID:勤務会社ID, GID:会社が所属する業界ID

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

【実行結果】

データの前処理を行います。

### データの追加的作成
# 会社インデックス、会社が所属する業界のインデックスの作成（0始まりにする）

# 会社が所属する業界IDを格納するデータフレーム
k2g_df = data2[['KID', 'GID']].drop_duplicates().sort_values(['KID']) - 1

【実行結果】省略

モデル式8-5の再掲です。

### モデルの定義 ◆モデル式8-5 model8-5.stan

with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    model3.add_coord('data', values=data2.index, mutable=True)
    model3.add_coord('kaisha', values=sorted(data2.KID.unique()), mutable=True)
    model3.add_coord('gyokai', values=sorted(data2.GID.unique()), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data2['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data2['X'].values, dims='data')
    # 説明変数 KIdx 会社インデックス 0始まりに変換
    KIdx = pm.ConstantData('KIdx', value=data2['KID'].values - 1, dims='data')
    # 説明変数 K2G 会社が所属する業界ID
    K2G = pm.ConstantData('K2G', value=k2g_df.GID.values, dims='kaisha')

    ### 事前分布
    a0 = pm.Uniform('a0', lower=-10000, upper=10000)
    b0 = pm.Uniform('b0', lower=-10000, upper=10000)
    sigmaAg = pm.Uniform('sigmaAg', lower=0, upper=10000)
    sigmaBg = pm.Uniform('sigmaBg', lower=0, upper=10000)
    sigmaA = pm.Uniform('sigmaA', lower=0, upper=10000)
    sigmaB = pm.Uniform('sigmaB', lower=0, upper=10000)
    sigmaY = pm.Uniform('sigmaY', lower=0, upper=10000)
    ag = pm.Normal('ag', mu=a0, sigma=sigmaAg, dims='gyokai')
    bg = pm.Normal('bg', mu=b0, sigma=sigmaBg, dims='gyokai')
    a = pm.Normal('a', mu=ag[K2G], sigma=sigmaA, dims='kaisha')
    b = pm.Normal('b', mu=bg[K2G], sigma=sigmaB, dims='kaisha')

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

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

MCMCを実行します。

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

【実行結果】省略

《解答》
図8.9を参考にして出題の事後分布を描画します。
まずはパラメータ a から。

### 業界平均のパラメータa[g]の事後分布のヴァイオリンプロットの描画 ◆解答案

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
ag_samples = idata3.posterior.ag.stack(sample=('chain', 'draw')).data
# データフレーム化
params_df = pd.DataFrame({
    'a_業界平均[1]': ag_samples[0],
    'a_業界平均[2]': ag_samples[1],
    'a_業界平均[3]': ag_samples[2],
})

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 5))
# ヴァイオリンプロットの描画
sns.violinplot(data=params_df, orient='h', fill=False, inner=None, linewidth=5,
               alpha=0.5, ax=ax)
# エラーバー付きポイントプロットの描画（ヴァイオリンプロットの中）
# estimatorで中央値、errorbarで95パーセンタイルインターバル(pi)を指定
sns.pointplot(data=params_df, estimator='median', errorbar=('pi', 95),
              orient='h', linestyle='none', color='tab:red', ms=10, ax=ax)
# 修飾
ax.set(xlabel='value', ylabel='Parameter')
plt.grid(lw=0.5);

【実行結果】

続いてパラメータ b 。

### 業界平均のパラメータb[g]の事後分布のヴァイオリンプロットの描画 ◆解答案

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
bg_samples = idata3.posterior.bg.stack(sample=('chain', 'draw')).data
# データフレーム化
params_df = pd.DataFrame({
    'b_業界平均[1]': bg_samples[0],
    'b_業界平均[2]': bg_samples[1],
    'b_業界平均[3]': bg_samples[2],
})

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 5))
# ヴァイオリンプロットの描画
sns.violinplot(data=params_df, orient='h', fill=False, inner=None, linewidth=5,
               alpha=0.5, ax=ax)
# エラーバー付きポイントプロットの描画（ヴァイオリンプロットの中）
# estimatorで中央値、errorbarで95パーセンタイルインターバル(pi)を指定
sns.pointplot(data=params_df, estimator='median', errorbar=('pi', 95),
              orient='h', linestyle='none', color='tab:red', ms=10, ax=ax)
# 修飾
ax.set(xlabel='value', ylabel='Parameter')
plt.grid(lw=0.5);

【実行結果】

練習問題(4)

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

### データの読み込み ◆データファイル8.4 data-attendance-4-1.txt
# PersonID:学生ID, A:バイト好き区分(1:好き), Score:学問の興味の強さ（0~200）

data3_1 = pd.read_csv('./data/data-attendance-4-1.txt')
print('data3_1.shape: ', data3_1.shape)
display(data3_1.head())

【実行結果】

### データの読み込み ◆データファイル8.5 data-attendance-4-2.txt
# PersonID:学生ID, CourseID:科目ID, Weather:天気（A:晴れ、B:曇り、C:雨）
# Y:出欠区分：授業に出席したかどうか（0:欠席、1:出席）

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

【実行結果】

《解答》
出題の２つのヒストグラムを描画します。
学生別出席率のヒストグラムを描画します。

### 学生別出席率のヒストグラムの描画 ◆解答例

# 学生別出席率の算出
ratio_p = data3_2.groupby(['PersonID'])['Y'].mean()
# binsの設定
binwidth = (ratio_p.max() - ratio_p.min()) / 30
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# ヒストグラムの描画
sns.histplot(ratio_p, binwidth=binwidth, ec='white', alpha=0.7, kde=True, ax=ax)
# 修飾
ax.set(title='学生ごとの出席率のヒストグラム', xlabel='出席率')
ax.grid(lw=0.5);

【実行結果】

科目別の出席率のヒストグラムを描画します。

### 科目別出席率のヒストグラムの描画 ◆解答例

# 科目別出席率の算出
ratio_c = data3_2.groupby(['CourseID'])['Y'].mean()
# binsの設定
binwidth = (ratio_c.max() - ratio_c.min()) / 30
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# ヒストグラムの描画
sns.histplot(ratio_c, binwidth=binwidth, ec='white', alpha=0.7, kde=True, ax=ax)
# 修飾
ax.set(title='科目ごとの出席率のヒストグラム', xlabel='出席率')
ax.grid(lw=0.5);

【実行結果】

練習問題(5)

既に読み込み済みのデータを使います。
データの前処理を行います。

### データの前処理：天気の重み係数を設定
Weather_dict2 = {'A': 0, 'B': 0.2, 'C': 1}
data3_2['Weather_w'] = data3_2['Weather'].map(Weather_dict2)
data3_2.head()

【実行結果】

出題のモデルを定義します。

### モデルの定義

with pm.Model() as model4:
    
    ### データ関連定義
    ## coordの定義
    model4.add_coord('data', values=data3_2.index, mutable=True)
    model4.add_coord('person', values=data3_1['PersonID'].values, mutable=True)
    model4.add_coord('course', values=sorted(data3_2['CourseID'].unique()),
                    mutable=True)
    model4.add_coord('coef', values=[1, 2, 3, 4], mutable=True)

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data3_2['Y'].values, dims='data')
    # 説明変数 A
    A = pm.ConstantData('A', value=data3_1['A'].values, dims='person')
    # 説明変数 Score
    Score = pm.ConstantData('Score', value=data3_1['Score'].values / 200,
                            dims='person')
    # 説明変数 Weather
    Weather = pm.ConstantData('Weather', value=data3_2['Weather_w'].values,
                             dims='data')
    # インデックス PersonIdx data2のPersonIDを0始まりに変換
    PersonIdx = pm.ConstantData('PersonIdx',
                                value=data3_2['PersonID'].values - 1,
                                dims='data')
    # インデックス CourseIdx data2のCourseIDを0始まりに変換
    CourseIdx = pm.ConstantData('CourseIdx',
                                value=data3_2['CourseID'].values - 1,
                                dims='data')
    
    ### 事前分布
    # 回帰係数
    b = pm.Uniform('b', lower=-10, upper=10, dims='coef')
    # 学生に依存する項
    bP = pm.Uniform('bP', lower=-10, upper=10, dims='person')
    xP = pm.Deterministic('xP', b[1]*A + b[2]*Score + bP, dims='person')
    # 科目に依存する項
    bC = pm.Uniform('bC', lower=-10, upper=10, dims='course')
    xC = pm.Deterministic('xC', bC, dims='course')
    # 授業に依存する項
    xJ = pm.Deterministic('xJ', b[3]*Weather, dims='data')

    ### 線形予測子 x  ※切片はこちらに集約
    x = pm.Deterministic('x', b[0] + xP[PersonIdx] + xC[CourseIdx] + xJ,
                         dims='data')
    ### 確率パラメータq ※線形予測子を逆ロジット変換
    q = pm.Deterministic('q', pm.invlogit(x), dims='data')

    ### 尤度関数
    obs = pm.Bernoulli('obs', p=q, observed=Y, dims='data')

モデルを確認します。

### モデルの表示
model4

【実行結果】

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

【実行結果】

MCMCを実行します。

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

【実行結果】省略

収束確認を行います。

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

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

【実行結果】
４本のマルコフ連鎖について、すべてのパラメータが$${\hat{R}<1.1}$$ですので、収束しているとみなします。

要約統計量を表示します。

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

【実行結果】

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

### トレースプロットの表示
pm.plot_trace(idata4, compact=False, var_names=['b'])
plt.tight_layout();

【実行結果】
b[1], b[2], b[3] は４本のマルコフ連鎖がバラバラであり、-10～10まで平坦な感じがします。

《解答》
問題文には具体的な解答形式が記載されていないので、統計量と可視化を掲載してみます。

事後統計量を算出します。

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

### 要約統計量の算出・表示
b_samples_df = pd.DataFrame(
            idata4.posterior.b.stack(sample=('chain', 'draw')).data.T,
            columns=[f'b[{i+1}]' for i in range(4)])
display(make_stats_df(b_samples_df).round(2))

【実行結果】

### 要約統計量の算出・表示
bP_samples_df = pd.DataFrame(
            idata4.posterior.bP.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bC[{i+1}]' for i in range(50)])
display(make_stats_df(bP_samples_df).round(2))

【実行結果】

### 要約統計量の算出・表示
bC_samples_df = pd.DataFrame(
            idata4.posterior.bC.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bC[{i+1}]' for i in range(10)])
display(make_stats_df(bC_samples_df).round(2))

【実行結果】

b_P（学生差）とb_C（科目差）のKDE曲線を描画します。

### パラメータb_Pの事後分布のKDEプロットの描画

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_P_samples = idata4.posterior.bP.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab20')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_C[n]のKDEプロットの描画
for i, b_C in enumerate(b_P_samples):
    # 色の設定
    color = cmap(i % 20)
    # x軸の値を設定
    x_ticks = np.linspace(b_C.min(), b_C.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_C, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.01), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'学生差のパラメータ $\beta_P$ の事後分布：KDEプロット')
# plt.ylim(-3, 45)
plt.grid(lw=0.5)

【実行結果】

### パラメータb_Cの事後分布のKDEプロットの描画

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_C_samples = idata4.posterior.bC.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab10')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_C[n]のKDEプロットの描画
for i, b_C in enumerate(b_C_samples):
    # 色の設定
    color = cmap(i / 10)
    # x軸の値を設定
    x_ticks = np.linspace(b_C.min(), b_C.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_C, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.01), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'科目差のパラメータ $\beta_C$ の事後分布：KDEプロット')
# plt.ylim(-3, 45)
plt.grid(lw=0.5)

【実行結果】

練習問題(6)

データを読み込みます。

### データの読み込み
# id:個体番号, y:調査した８個中、生存していた種子数（目的変数）
data4 = pd.read_csv('./data/data7a.csv')
print('data4.shape: ', data4.shape)
display(data4.head())

【実行結果】

Y のヒストグラムを描画します。

### yのヒストグラムの描画
sns.histplot(data=data4, x='y', bins=9, ec='white')
plt.grid(lw=0.5);

【実行結果】

PyMCのモデルを定義します。

### モデルの定義

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

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data4['y'].values, dims='data')
    # インデックス idIdx idを0始まりに変換
    idIdx = pm.ConstantData('idIdx', value=data4['id'].values - 1, dims='data')
    
    ### 事前分布
    # 回帰係数
    b = pm.Uniform('b', lower=-100, upper=100)
    sigmaI = pm.Uniform('sigmaI', lower=0, upper=100)
    bI = pm.Normal('bI', mu=0, sigma=sigmaI, dims='data')

    ### 線形予測子
    x = pm.Deterministic('x', b + bI[idIdx], dims='data')
    ### 確率パラメータq ※線形予測子を逆ロジット変換
    q = pm.Deterministic('q', pm.invlogit(x), dims='data')

    ### 尤度関数
    obs = pm.Binomial('obs', p=q, n=8, observed=Y, dims='data')

モデルの定義内容を確認します。

### モデルの表示
model5

【実行結果】

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

【実行結果】

MCMCを実行します。

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

【実行結果】省略

収束確認を行います。

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

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

【実行結果】
４本のマルコフ連鎖について、すべてのパラメータが$${\hat{R}<1.1}$$ですので、収束しているとみなします。

要約統計量を表示します。

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

【実行結果】

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

### トレースプロットの表示
var_names = ['b', 'sigmaI', 'bI', 'x', 'q']
pm.plot_trace(idata5, 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

### 要約統計量の算出・表示
vars = ['b', 'sigmaI']
param_samples = idata5.posterior[vars].to_dataframe().reset_index(drop=True)
bI_samples_df = pd.DataFrame(
            idata5.posterior.bI.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bI[{i+1}]' for i in range(data4.shape[0])])
param_samples = pd.concat([param_samples, bI_samples_df], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】

個体差の大きさ・パラメータ b_I の事後分布のKDEプロットを描画します。

### パラメータb_Iの事後分布のKDEプロットの描画

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_I_samples = idata5.posterior.bI.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab20')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_I[n]のKDEプロットの描画
for i, b_I in enumerate(b_I_samples):
    # 色の設定
    color = cmap(i % 20)
    # x軸の値を設定
    x_ticks = np.linspace(b_I.min(), b_I.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_I, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.02), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'個体差のパラメータ $\beta_I$ の事後分布：KDEプロット')
plt.xlim(-10, 10)
plt.grid(lw=0.5)

【実行結果】

練習問題(7)

データを読み込みます。

### データの読み込み
# id:個体番号, pot:植木鉢番号, f:処理の違い, y:種子数（目的変数）

data5 = pd.read_csv('./data/d1.csv')
print('data5.shape: ', data5.shape)
display(data5.head())

【実行結果】

pot と f の要素を確認します。

### potとfの値および要素数を把握
print(data5['pot'].value_counts())
print(data5['f'].value_counts())

【実行結果】

Y の箱ひげ図を描画します。

### 箱ひげ図の描画
sns.boxplot(data=data5, x='pot', y='y', hue='f', fill=False)
sns.swarmplot(data=data5, x='pot', y='y', hue='f', alpha=0.5, legend=None)
plt.grid(lw=0.5)
plt.legend(loc='upper left', title='f:処理');

【実行結果】

PyMCのモデルを定義します。

### モデルの定義

## インデックス化
# pot：0始まりの整数値に変換
pot_cat = sorted(data5['pot'].unique())
pot_idx = data5['pot'].apply(lambda x: pot_cat.index(x))
# f：0始まりの整数値に変換
f_cat = ['C', 'T']
f_idx = data5['f'].apply(lambda x: f_cat.index(x))

## F: 処理fの辞書
F = {'C': 0, 'T': 1}

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

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data5['y'].values, dims='data')
    # インデックス fIdx fをC:0, T:1の数値に変換
    fIdx = pm.ConstantData('fIdx', value=f_idx, dims='data')
    # インデックス potIdx potを0始まりの数値に変換
    potIdx = pm.ConstantData('potIdx', value=pot_idx, dims='data')
    
    ### 事前分布
    # 切片
    b1 = pm.Uniform('b1', lower=-100, upper=100)
    # 処理差
    b2 = pm.Uniform('b2', lower=-100, upper=100)
    # 個体差
    sigmaI = pm.Uniform('sigmaI', lower=0, upper=100)
    bI = pm.Normal('bI', mu=0, sigma=sigmaI, dims='data')
    # 植木鉢差
    sigmaP = pm.Uniform('sigmaP', lower=0, upper=100)
    bP = pm.Normal('bP', mu=0, sigma=sigmaP, dims='pot')

    ### 線形予測子
    loglam = pm.Deterministic('loglam', b1 + b2*fIdx + bI + bP[potIdx],
                              dims='data')
    lam = pm.Deterministic('lam', pt.exp(loglam), dims='data')

    ### 尤度関数
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='data')

モデルの定義内容を確認します。

### モデルの表示
model6

【実行結果】

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

【実行結果】

MCMCを実行します。

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

【実行結果】省略

収束確認を行います。

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

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

【実行結果】
４本のマルコフ連鎖について、すべてのパラメータが$${\hat{R}<1.1}$$ですので、収束しているとみなします。

要約統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['b1', 'b2', 'sigmaP', 'sigmaI', 'bP', 'bI']
pm.summary(idata6, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

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

### トレースプロットの表示
var_names = ['b1', 'b2', 'sigmaP', 'sigmaI', 'bP', 'bI', 'lam']
pm.plot_trace(idata6, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

《解答》
問題文には具体的な解答形式が記載されていないので、統計量と可視化を掲載してみます。

事後統計量を算出します。

### 要約統計量の算出・表示
vars = ['b1', 'b2', 'sigmaP', 'sigmaI']
param_samples = idata6.posterior[vars].to_dataframe().reset_index(drop=True)
bP_samples_df = pd.DataFrame(
            idata6.posterior.bP.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bP[{i+1}]' for i in range(len(pot_cat))])
param_samples = pd.concat([param_samples, bP_samples_df], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】

### 要約統計量の算出・表示 ※図8.9左のパラメータ
bI_samples_df = pd.DataFrame(
            idata6.posterior.bI.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bI[{i+1}]' for i in range(data5.shape[0])])
display(make_stats_df(bI_samples_df).round(2))

【実行結果】

個体差の大きさ・パラメータ b_I の事後分布のKDEプロットを描画します。

### パラメータb_Iの事後分布のKDEプロットの描画

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_I_samples = idata6.posterior.bI.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab20')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_I[n]のKDEプロットの描画
for i, b_I in enumerate(b_I_samples):
    # 色の設定
    color = cmap(i % 20)
    # x軸の値を設定
    x_ticks = np.linspace(b_I.min(), b_I.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_I, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.02), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'個体差のパラメータ $\beta_I$ の事後分布：KDEプロット')
plt.xlim(-5, 5)
plt.grid(lw=0.5)

【実行結果】

植木鉢差の大きさ・パラメータ b_P の事後分布のKDEプロットを描画します

### パラメータb_Pの事後分布のKDEプロットの描画

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_P_samples = idata6.posterior.bP.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab10')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_P[pot]のKDEプロットの描画
for i, b_P in enumerate(b_P_samples):
    # 色の設定
    color = cmap(i / 10)
    # x軸の値を設定
    x_ticks = np.linspace(b_P.min(), b_P.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_P, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.02), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'植木鉢差のパラメータ $\beta_P$ の事後分布：KDEプロット')
plt.grid(lw=0.5)

【実行結果】

第８章練習問題は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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件

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第８章「練習問題」

第８章「階層モデル」

はじめに

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

引用表記

PyMC環境の準備

8章 練習問題

インポート

練習問題(1)

練習問題(2)

練習問題(3)

練習問題(4)

練習問題(5)

練習問題(6)

練習問題(7)

シリーズの記事

ブログの紹介

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

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

8章練習問題