見出し画像

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

ネイピア DS

第５章「基本的な回帰とモデルのチェック」

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

この記事は、テキスト第５章「基本的な回帰とモデルのチェック」の「練習問題」の 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を動かすまでの準備」章をご覧ください。


5章 練習問題

インポート

### インポート

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

5.1 重回帰の関連問題

① 練習問題 (1)

5.1節で利用したサンプルコードのデータを再度、読み込みます。

### データの読み込み ◆データファイル5.1 data-attendance-1.txtの構成
# A:バイト好き区分(1:好き), Score:学問の興味の強さ（0~200）, Y:授業の出席率（1年間）

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

【実行結果】

モデルを定義します。5.1節のモデルと同じです。

### モデルの定義

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

    ### 事前分布
    b1 = pm.Uniform('b1', lower=-10, upper=10)
    b2 = pm.Uniform('b2', lower=-10, upper=10)
    b3 = pm.Uniform('b3', lower=-10, upper=10)
    sigma = pm.Uniform('sigma', lower=0, upper=10)

    ### 線形予測子
    mu = pm.Deterministic('mu', b1 + b2 * A + b3 * Score, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

PythonでMCMCを実行します。

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

【実行結果】省略

練習問題 (1) の解答のコードを描きます。
推論データ idata1 の mu の事後分布と観測値 Y の値を用いて 誤差項 ε の事後分布を算出します。

### 練習問題(1) 解答

# 推論データidataからmuの事後分布サンプルを取り出し shape=(4000, 50)
mu_samples = idata1.posterior.mu.stack(s=('chain', 'draw')).T.data
# εの算出
epsilon = np.array([data1.Y - mu_sample for mu_sample in mu_samples])
# εの概要の表示
print('epsilon.shape: ', epsilon.shape)
print('epsilon[0, :]: \n', epsilon[0, :])

【実行結果】

② 練習問題 (2)

MCMCで ε の事後分布サンプリングを行えるようにモデルを作成します。

### モデルの定義

with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data1.index, mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
    # 説明変数 A
    A = pm.ConstantData('A', value=data1['A'].values, dims='data')
    # 説明変数 Score / 200
    Score = pm.ConstantData('Score', value=data1['Score'].values / 200,
                            dims='data')

    ### 事前分布
    b1 = pm.Uniform('b1', lower=-10, upper=10)
    b2 = pm.Uniform('b2', lower=-10, upper=10)
    b3 = pm.Uniform('b3', lower=-10, upper=10)
    sigma = pm.Uniform('sigma', lower=0, upper=10)

    ### 線形予測子
    mu = pm.Deterministic('mu', b1 + b2 * A + b3 * Score, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

    ### 計算値
    epsilon = pm.Deterministic('epsilon', Y - mu, dims='data')

PythonでMCMCを実行します。

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

【実行結果】省略

練習問題 (2) の解答のコードを描きます。
推論データ idata1 の ε の事後分布サンプリングデータを取り出します。

### epsilonの概要の表示

# 推論データidataからepsilonの事後分布サンプルを取り出し shape=(4000, 50)
epsilon_samples = idata2.posterior.epsilon.stack(s=('chain', 'draw')).data.T
# εの概要の表示
print('epsilon_samples.shape: ', epsilon_samples.shape)
print('epsilon_samples[0, :]: \n', epsilon_samples[0, :])

【実行結果】
練習問題 (1) の ε の値と比べると驚きの事実が！

5.3 ロジスティック回帰の関連問題

③ 練習問題 (3)

5.3節で利用したサンプルコードのデータを再度、読み込みます。

### データの読み込み ◆データファイル5.3 data-attendance-3.txtの構成
# PersonID:学生ID、A:バイト好き区分(1:好き), Weather:天気（A:晴れ、B:曇り、C:雨）
# Y:出欠区分：授業に出席したかどうか（0:欠席、1:出席）

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

【実行結果】

練習問題 (3) の解答のコードを描きます。
pandas の pivot_table() で集計します。

### 練習問題(3) 解答
### アルバイトが好きかどうかごとにYを集計する
data2.pivot_table(index='A', columns='Y', values='PersonID', aggfunc='count')

【実行結果】

④ 練習問題 (4)

以下のコードの全体で練習問題 (4) を解きます。

天気のインデックスを作成します。

### 天気のインデックスの作成
Weather_dict2 = {'A': 0, 'B': 1, 'C': 2}
data2['Weather_idx'] = data2['Weather'].map(Weather_dict2)
data2.head()

【実行結果】

モデルを定義します。5.1節のモデルと同じです。
パラメータ bw を追加しています。

### モデルの定義

with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    model3.add_coord('data', values=data2.index, mutable=True)
    model3.add_coord('beta', values=[1, 2, 3], mutable=True)
    model3.add_coord('weather', values=['晴れ', '曇り', '雨'], mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data2['Y'].values, dims='data')
    # 説明変数 A
    A = pm.ConstantData('A', value=data2['A'].values, dims='data')
    # 説明変数 Score / 200
    Score = pm.ConstantData('Score', value=data2['Score'].values / 200,
                            dims='data')
    # 説明変数 WIdx (Weather_idx)
    WIdx = pm.ConstantData('WIdx', value=data2['Weather_idx'].values, 
                           dims='data')

    ### 事前分布
    ## β0~β2
    b = pm.Uniform('b', lower=-10, upper=10, dims='beta')
    ## 天気の切片 bw
    # 曇りの切片bw2と雨の切片bw3
    bw2 = pm.Uniform('bw2', lower=-10, upper=10)
    bw3 = pm.Uniform('bw3', lower=-10, upper=10)
    # bw （晴れ0, 曇りbw2, 雨bw3
    bw = pt.stack([0, bw2, bw3])

    ### 線形予測子の逆ロジット変換
    q = pm.Deterministic(
        'q',
        pm.invlogit(b[0] + b[1] * A + b[2] * Score + bw[WIdx]),
        dims='data')
    
    ### 尤度関数
    obs = pm.Bernoulli('obs', p=q, observed=Y, dims='data')

モデルの概観を確認します。

### モデルの表示
model3

【実行結果】

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

【実行結果】

MCMCを実行します。

### 事後分布からのサンプリング 1分10秒
with model3:
    idata3 = 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 = idata3        # idata名
threshold = 1.01         # しきい値

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

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

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

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

【実行結果】

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

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

【実行結果】

練習問題 (4) の解答の最後のコードを描きます。
パラメータの要約統計量の表示です。

### パラメータの要約を確認

## 統計量算出関数：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

## 要約統計量の算出・表示
# 事後分布サンプリングデータidataからパラメータbを取り出してデータフレーム化
param_samples = pd.DataFrame(
    idata3.posterior.b.stack(sample=('chain', 'draw')).T.data,
    columns=['b1', 'b2', 'b3'])
param_samples = pd.concat(
    [param_samples, 
    idata3.posterior[['bw2', 'bw3']].to_dataframe().reset_index(drop=True)],
    axis=1)
# 上記データフレームを統計量算出関数に与えて事後統計量データフレームを作成
params_stats_df = make_stats_df(param_samples)
# 事後統計量データフレームの表示
display(params_stats_df.round(2))

【実行結果】

5.4 ポアソン回帰の関連問題

⑤ 練習問題 (5)

5.4節で利用したサンプルコードのデータを再度、読み込みます。

### データの読み込み ◆データファイル5.2 data-attendance-2.txtの構成
# PersonID:学生ID、A:バイト好き区分(1:好き), Score:学問の興味の強さ（0~200）
# M:総授業回数（3か月間）, Y:授業出席回数（3か月間）

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

【実行結果】

5.4 節と同じモデルでMCMCを実行し、事後予測データを取得します。

### モデルの定義

with pm.Model() as model4:
    
    ### データ関連定義
    ## coordの定義
    model4.add_coord('data', values=data3.index, mutable=True)
    model4.add_coord('beta', values=[1, 2, 3], mutable=True)
    ## dataの定義
    # 目的変数 M
    M = pm.ConstantData('M', value=data3['M'].values, dims='data')
    # 説明変数 A
    A = pm.ConstantData('A', value=data3['A'].values, dims='data')
    # 説明変数 Score / 200
    Score = pm.ConstantData('Score', value=data3['Score'].values / 200,
                            dims='data')

    ### 事前分布
    b = pm.Uniform('b', lower=-10, upper=10, dims='beta')

    ### lambda 線形予測子の指数変換
    lam = pm.Deterministic('q', pt.exp(b[0] + b[1] * A + b[2] * Score),
                           dims='data')
    
    ### 尤度関数
    obs = pm.Poisson('obs', mu=lam, observed=M, dims='data')
### 事後分布からのサンプリング 10秒
with model4:
    idata4 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

### Yの事後予測分布のサンプリング
with model4:
    idata4.extend(pm.sample_posterior_predictive(idata4, random_seed=1234))

【実行結果】省略

### ppcプロットの描画
pm.plot_ppc(idata4, num_pp_samples=100);

【実行結果】

練習問題 (5) の解答のコードを描きます。
観測値と予測値のプロットを描画します。

### 観測値と予測値のプロット

## 描画用データの作成 yPredの個人別の中央値と80%区間を算出
# 事後予測サンプリングデータからobsを取り出し
m_pred_samples = (idata4.posterior_predictive.obs
                  .stack(sample=('chain', 'draw')).data)
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
m_pred_df = pd.DataFrame(
    np.quantile(m_pred_samples, q=[0.1, 0.5, 0.9], axis=1).T,
    columns=['10%', 'median', '90%'])
# 観測値dataと予測値m_pred_dfのデータフレームを結合
m_pred_df = pd.concat([data3, m_pred_df], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
m_pred_df['err_lower'] = m_pred_df['median'] - m_pred_df['10%'] 
m_pred_df['err_upper'] = m_pred_df['90%']  - m_pred_df['median']
# アルバイト0とアルバイト1に分離
m_pred_A0 = m_pred_df[m_pred_df['A']==0]
m_pred_A1 = m_pred_df[m_pred_df['A']==1]

## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# アルバイト0の描画（エラーバー付き散布図）
ax.errorbar(m_pred_A0['M'], m_pred_A0['median'],
            yerr=[m_pred_A0['err_lower'], m_pred_A0['err_upper']],
            color='tab:blue', alpha=0.7, marker='o', ms=10, linestyle='none',
            label='0')
# アルバイト1の描画（エラーバー付き散布図）
ax.errorbar(m_pred_A1['M'], m_pred_A1['median'],
            yerr=[m_pred_A1['err_lower'], m_pred_A1['err_upper']],
            color='tab:orange', alpha=0.7, marker='^', ms=10, linestyle='none',
            label='1')
# 赤い対角線の描画
ax.plot([10, 90], [10, 90], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $M$の観測値', ylabel='Predicted: $M$の予測値(中央値)',
       title='$M$ の観測値と予測値（中央値）のプロット 80%区間バー付き')
ax.legend(title='アルバイト好き', bbox_to_anchor=(1, 1))
ax.grid(lw=0.5);

【実行結果】

上記の図はテキストの模範解答図（GitHub）とほぼ同じ結果になりました。
予測値は芳しくない感じがします。

練習問題 (6)

GitHub より練習問題用データを読み込んで、ポアソン回帰を実行します。

データを読み込みます。

### データの読み込み
# y:種子数（目的変数）, x:体サイズ, f:処理の違い 
data4 = pd.read_csv('./data/data3a.csv')
print('data4.shape: ', data4.shape)
display(data4.head())

【実行結果】

変数 f の要素を確認します。

### fの分布確認
data4['f'].value_counts()

【実行結果】
C と T が半分ずつです！

散布図行列でデータを確認します。

### 散布図行列の描画

## 描画領域の指定
fig, ax = plt.subplots(3, 3, figsize=(8, 8))
ax = ax.ravel() # １次元でaxesを指定したいので

## 番地0,0：ヒストグラムの描画
sns.histplot(ax=ax[0], data=data4, x='y', hue='f', bins=10, kde=True,
             ec='white', legend=None)
ax[0].set(xlabel=None, ylabel='y')
ax[0].grid(lw=0.5)

## 番地1,0：散布図の描画
sns.scatterplot(ax=ax[3], data=data4, x='y', y='x', hue='f',
                size='f', style='f', markers=['o', '^'], sizes=(80, 80),
                alpha=0.5, legend=None)
ax[3].set(xlabel=None)
ax[3].grid(lw=0.5)

## 番地1,1：ヒストグラムの描画
sns.histplot(ax=ax[4], data=data4, x='x', hue='f', bins=10, kde=True,
             ec='white', legend=None)
ax[4].set(xlabel=None, ylabel=None)
ax[4].grid(lw=0.5)

## 番地2,0：箱ひげ図＋ストリッププロットの描画
sns.boxplot(ax=ax[6], x=data4.y, y=data4.f, hue=data4.f, fill=False,
            legend=None)
sns.stripplot(ax=ax[6], x=data4.y, y=data4.f, hue=data4.f, size=8, alpha=0.4,
              legend=None)
ax[6].grid(lw=0.5)

## 番地2,1：箱ひげ図＋ストリッププロットの描画
sns.boxplot(ax=ax[7], x=data4.x, y=data4.f, hue=data4.f, fill=False,
            legend=None)
sns.stripplot(ax=ax[7], x=data4.x, y=data4.f, hue=data4.f, size=8, alpha=0.4,
              legend=None)
ax[7].set(ylabel=None)
ax[7].grid(lw=0.5)

## 番地2,2：ヒストグラムの描画（棒グラフを使用）
bar_f = data4.f.value_counts().sort_index()
sns.barplot(ax=ax[8], x=bar_f.index, y=bar_f, hue=bar_f.index, palette='tab10',
            alpha=0.5, ec='white', legend=None)
ax[8].set(xlabel='f', ylabel=None)
ax[8].grid(lw=0.5)

## スピアマンの順位相関係数を上三角のaxesに表示
# 列名をリスト化
cols = data4.columns
# 列名の組み合わせ 行i,列j ごとにテキスト表示を繰り返す
for i, col1 in enumerate(cols):
    for j, col2 in enumerate(cols):
        # 上三角の位置は 行i < 列j のとき
        if i < j:
            # axesの番号を取得
            pos = i * len(cols) + j
            # スピアマンの順位相関係数を算出
            corr, pval = stats.spearmanr(data4[col1], data4[col2])
            # 枠線等を削除
            ax[pos].set_axis_off()
            # テキスト表示：中央表示に関連する引数: x,y,va,ha,transform
            ax[pos].text(x=0.5, y=0.5, s=round(corr * 100), fontsize=30,
                         va='center', ha='center', transform=ax[pos].transAxes)

plt.tight_layout();

【実行結果】

モデルを定義します。

### モデルの定義

# fのインデックス化 C=0, T=1
f_dict = {'C': 0, 'T': 1}
f_idx = data4['f'].map(f_dict)

# モデルの定義
with pm.Model() as model5:
    
    ### データ関連定義
    ## coordの定義
    # データのインデックス
    model5.add_coord('data', values=data4.index, mutable=True)
    # パラメータβの添字番号
    model5.add_coord('beta', values=[1, 2, 3], mutable=True)
    # 変数fのカテゴリ値(C, T)
    model5.add_coord('f', values=f_dict.keys(), mutable=True)
    ## dataの定義
    # 目的変数 y
    y = pm.ConstantData('y', value=data4['y'].values, dims='data')
    # 説明変数 x
    x = pm.ConstantData('x', value=data4['x'].values, dims='data')
    # 説明変数 fIdx
    fIdx = pm.ConstantData('fIdx', value=f_idx, dims='data')

    ### 事前分布
    b = pm.Uniform('b', lower=-10, upper=10, dims='beta')

    ### lambda 線形予測子の指数変換
    lam = pm.Deterministic('q', pt.exp(b[0] + b[1] * x + b[2] * fIdx),
                           dims='data')
    
    ### 尤度関数
    obs = pm.Poisson('obs', mu=lam, 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)

【実行結果】省略

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

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

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

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

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

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

【実行結果】

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

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

【実行結果】

目的変数 Y の事後予測サンプリングを実行します。

### Yの事後予測分布のサンプリング
with model5:
    idata5.extend(pm.sample_posterior_predictive(idata5, random_seed=1234))

【実行結果】省略

ppcプロットを描画します。

### ppcプロットの描画
pm.plot_ppc(idata5, num_pp_samples=100);

【実行結果】

パラメータの事後統計量を表示します。

### パラメータの要約を確認

## 統計量算出関数：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

## 要約統計量の算出・表示
# 事後分布サンプリングデータidataからパラメータbを取り出してデータフレーム化
param_samples = pd.DataFrame(
    idata5.posterior.b.stack(sample=('chain', 'draw')).T.data,
    columns=['b1', 'b2', 'b3'])
# 上記データフレームを統計量算出関数に与えて事後統計量データフレームを作成
params_stats_df = make_stats_df(param_samples)
# 事後統計量データフレームの表示
display(params_stats_df.round(2))

【実行結果】

目的変数 Y について、観測値と事後予測を描画します。

### 観測値と予測値のプロット

## 描画用データの作成 yPredの個人別の中央値と80%区間を算出
# 事後予測サンプリングデータからobsを取り出し
y_pred_samples = (idata5.posterior_predictive.obs
                  .stack(sample=('chain', 'draw')).data)
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df = pd.DataFrame(
    np.quantile(y_pred_samples, q=[0.1, 0.5, 0.9], axis=1).T,
    columns=['10%', 'median', '90%'])
# 観測値data4と予測値y_pred_dfのデータフレームを結合
y_pred_df = pd.concat([data4, y_pred_df], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
y_pred_df['err_lower'] = y_pred_df['median'] - y_pred_df['10%'] 
y_pred_df['err_upper'] = y_pred_df['90%']  - y_pred_df['median']
# fのCとTでy_pred_dfを分離
y_pred_fC = y_pred_df[y_pred_df['f']=='C']
y_pred_fT = y_pred_df[y_pred_df['f']=='T']

## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# f=Cのy_predの描画（エラーバー付き散布図）
ax.errorbar(y_pred_fC['y'], y_pred_fC['median'],
            yerr=[y_pred_fC['err_lower'], y_pred_fC['err_upper']],
            color='tab:blue', alpha=0.7, marker='o', ms=10, linestyle='none',
            label='C')
# f=Tのy_predの描画（エラーバー付き散布図）
ax.errorbar(y_pred_fT['y'], y_pred_fT['median'],
            yerr=[y_pred_fT['err_lower'], y_pred_fT['err_upper']],
            color='tab:orange', alpha=0.7, marker='^', ms=10, linestyle='none',
            label='T')
# 赤い対角線の描画
ax.plot([0, 20], [0, 20], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $y$の観測値', ylabel='Predicted: $y$の予測値(中央値)',
       title='$y$ の観測値と予測値（中央値）のプロット 80%区間バー付き')
ax.legend(title='f', bbox_to_anchor=(1, 1))
ax.grid(lw=0.5);

【実行結果】
いまいちな予測になったようです。
この記事のモデリングに問題があるのか、この図が妥当なのかが、正解の図が公開されていないようなので、確認できません（汗）
確認できました。本モデリングの場合、この予測値になるようです。

MCMCサンプリングデータの散布図行列を描画します。

### MCMCサンプルの散布図行列の描画

## 描画用データの作成
# MCMCサンプリングデータからq1, q100を取り出し
q1_samples = (idata5.posterior['q'].to_dataframe().reset_index()
              .query('data==0').rename({'q': 'q1'}, axis=1))
q100_samples = (idata5.posterior['q'].to_dataframe().reset_index()
               .query('data==99').rename({'q': 'q100'}, axis=1))
# 描画対象パラメータをデータフレーム化
plot_df = pd.concat([param_samples,
                     q1_samples.reset_index(drop=True)['q1'],
                     q100_samples.reset_index(drop=True)['q100']], axis=1)

## 描画処理
# 相関行列プロットの描画
g = sns.pairplot(plot_df, diag_kws={'kde': True, 'ec': 'white'})
# スピアマンの順位相関係数の表示のためのaxフラット化
ax = g.axes.ravel()

## スピアマンの順位相関係数を上三角のaxesに表示
# 列名をリスト化
cols = plot_df.columns
# 列名の組み合わせ行i, 列j ごとにテキスト表示を繰り返す
for i, col1 in enumerate(cols):
    for j, col2 in enumerate(cols):
        # 上三角の位置は 行i < 列j のとき
        if i < j:
            # axesの番号を取得
            pos = i * len(cols) + j
            # スピアマンの順位相関係数を算出
            corr, pval = stats.spearmanr(plot_df[col1], plot_df[col2])
            # テキスト表示：中央表示に関連する引数: x,y,va,ha,transform
            ax[pos].text(x=0.5, y=0.5, s=round(corr * 100), fontsize=30,
                         va='center', ha='center', transform=ax[pos].transAxes,
                         bbox=dict(boxstyle='round', facecolor='white'))

【実行結果】

練習問題 (7)

GitHub より練習問題用データを読み込んで、二項ロジスティック回帰を実行します。

データを読み込みます。

### データの読み込み
# N:種子数, y:生存種子数（目的変数）, x:体サイズ, f:処理の違い 
data5 = pd.read_csv('./data/data4a.csv')
print('data5.shape: ', data5.shape)
display(data5.head())

【実行結果】

変数 N の要素を確認します。

### Nの要素を確認
data5['N'].value_counts()

【実行結果】
種子数 N はオール８（総種子数が８つ）です。

散布図行列でデータを確認します。

### 散布図行列の描画

## 描画領域の指定
fig, ax = plt.subplots(4, 4, figsize=(10, 10))
ax = ax.ravel() # １次元でaxesを指定したいので

## 番地0,0：ヒストグラムの描画
sns.histplot(ax=ax[0], data=data5, x='N', hue='f', bins=10, ec='white',
             legend=None)
ax[0].set(xlabel=None, ylabel='N')
ax[0].grid(lw=0.5)

## 番地1,0：散布図の描画
sns.scatterplot(ax=ax[4], data=data5, x='N', y='y', hue='f',
                size='f', style='f', markers=['o', '^'], sizes=(80, 80),
                alpha=0.5, legend=None)
ax[4].set(xlabel=None)
ax[4].grid(lw=0.5)

## 番地1,1：ヒストグラムの描画
sns.histplot(ax=ax[5], data=data5, x='y', hue='f', bins=10, kde=True,
             ec='white', legend=None)
ax[5].set(xlabel=None, ylabel=None)
ax[5].grid(lw=0.5)

## 番地2,0：散布図の描画
sns.scatterplot(ax=ax[8], data=data5, x='N', y='x', hue='f',
                size='f', style='f', markers=['o', '^'], sizes=(80, 80),
                alpha=0.5, legend=None)
ax[8].set(xlabel=None)
ax[8].grid(lw=0.5)

## 番地2,1：散布図の描画
sns.scatterplot(ax=ax[9], data=data5, x='y', y='x', hue='f',
                size='f', style='f', markers=['o', '^'], sizes=(80, 80),
                alpha=0.5, legend=None)
ax[9].set(xlabel=None)
ax[9].grid(lw=0.5)

## 番地2,2：ヒストグラムの描画
sns.histplot(ax=ax[10], data=data5, x='x', hue='f', bins=10, kde=True,
             ec='white', legend=None)
ax[10].set(xlabel=None, ylabel=None)
ax[10].grid(lw=0.5)

## 番地3,0：箱ひげ図＋ストリッププロットの描画
sns.boxplot(ax=ax[12], x=data5.N, y=data5.f, hue=data5.f, fill=False,
            legend=None)
sns.stripplot(ax=ax[12], x=data5.N, y=data5.f, hue=data5.f, size=8, alpha=0.4,
              legend=None)
ax[12].grid(lw=0.5)

## 番地3,1：箱ひげ図＋ストリッププロットの描画
sns.boxplot(ax=ax[13], x=data5.y, y=data5.f, hue=data5.f, fill=False,
            legend=None)
sns.stripplot(ax=ax[13], x=data5.y, y=data5.f, hue=data5.f, size=8, alpha=0.4,
              legend=None)
ax[13].grid(lw=0.5)

## 番地3,2：箱ひげ図＋ストリッププロットの描画
sns.boxplot(ax=ax[14], x=data5.x, y=data5.f, hue=data5.f, fill=False,
            legend=None)
sns.stripplot(ax=ax[14], x=data5.x, y=data5.f, hue=data5.f, size=8, alpha=0.4,
              legend=None)
ax[14].grid(lw=0.5)

## 番地3,3：ヒストグラムの描画（棒グラフを使用）
bar_f = data5.f.value_counts().sort_index()
sns.barplot(ax=ax[15], x=bar_f.index, y=bar_f, hue=bar_f.index, palette='tab10',
            alpha=0.5, ec='white', legend=None)
ax[15].set(xlabel='f', ylabel=None)
ax[15].grid(lw=0.5)

## スピアマンの順位相関係数を上三角のaxesに表示
# 列名をリスト化
cols = data5.columns
# 列名の組み合わせ 行i,列j ごとにテキスト表示を繰り返す
for i, col1 in enumerate(cols):
    for j, col2 in enumerate(cols):
        # 上三角の位置は 行i < 列j のとき
        if i < j:
            # axesの番号を取得
            pos = i * len(cols) + j
            # スピアマンの順位相関係数を算出
            corr, pval = stats.spearmanr(data5[col1], data5[col2])
            # 枠線等を削除
            ax[pos].set_axis_off()
            # テキスト表示：中央表示に関連する引数: x,y,va,ha,transform
            if np.isnan(corr):
                s = '-'
            else:
                s = round(corr * 100)
            ax[pos].text(x=0.5, y=0.5, s=s, fontsize=30,
                         va='center', ha='center', transform=ax[pos].transAxes)

plt.tight_layout();

【実行結果】
x と Y は正の相関がありそうです。

モデルを定義します。

### モデルの定義

# fのインデックス化 C=0, T=1
f_dict = {'C': 0, 'T': 1}
f_idx = data5['f'].map(f_dict)

# モデルの定義
with pm.Model() as model6:
    
    ### データ関連定義
    ## coordの定義
    # データのインデックス
    model6.add_coord('data', values=data5.index, mutable=True)
    # パラメータβの添字番号
    model6.add_coord('beta', values=[1, 2, 3], mutable=True)
    # 変数fのカテゴリ値(C, T)
    model6.add_coord('f', values=f_dict.keys(), mutable=True)
    ## dataの定義
    # 目的変数 y
    y = pm.ConstantData('y', value=data5['y'].values, dims='data')
    # 総種子数 N=8
    N = pm.ConstantData('N', value=data5['N'].values, dims='data')
    # 説明変数 x
    x = pm.ConstantData('x', value=data5['x'].values, dims='data')
    # 説明変数 fIdx
    fIdx = pm.ConstantData('fIdx', value=f_idx, dims='data')

    ### 事前分布
    b = pm.Uniform('b', lower=-50, upper=50, dims='beta')

    ### 確率q 線形予測子の指数変換
    q = pm.Deterministic('q', pm.invlogit(b[0] + b[1] * x + b[2] * fIdx),
                         dims='data')
    
    ### 尤度関数
    obs = pm.Binomial('obs', n=N, p=q, observed=y, dims='data')

モデルの概観を確認します。

### モデルの表示
model6

【実行結果】

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

【実行結果】

MCMCを実行します。

### 事後分布からのサンプリング 15秒
with model6:
    idata6 = 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 = idata6        # idata名
threshold = 1.01         # しきい値

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

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

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

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

【実行結果】

トレースプロットを描画します。
q は一部のみの表示です。

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

【実行結果】

目的変数 Y の事後予測サンプリングを実行します。

### Yの事後予測分布のサンプリング
with model6:
    idata6.extend(pm.sample_posterior_predictive(idata6, random_seed=1234))

【実行結果】省略

ppcプロットを描画します。

### ppcプロットの描画
pm.plot_ppc(idata6, num_pp_samples=100);

【実行結果】

パラメータの事後統計量を表示します。

### パラメータの要約を確認

## 統計量算出関数：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

## 要約統計量の算出・表示
# 事後分布サンプリングデータidataからパラメータbを取り出してデータフレーム化
param_samples = pd.DataFrame(
    idata6.posterior.b.stack(sample=('chain', 'draw')).T.data,
    columns=['b1', 'b2', 'b3'])
# 上記データフレームを統計量算出関数に与えて事後統計量データフレームを作成
params_stats_df = make_stats_df(param_samples)
# 事後統計量データフレームの表示
display(params_stats_df.round(2))

【実行結果】

目的変数 Y について、観測値と事後予測を描画します。

### 観測値と予測値のプロット

## 描画用データの作成 yPredの個人別の中央値と80%区間を算出
# 事後予測サンプリングデータからobsを取り出し
y_pred_samples = (idata6.posterior_predictive.obs
                  .stack(sample=('chain', 'draw')).data)
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df = pd.DataFrame(
    np.quantile(y_pred_samples, q=[0.1, 0.5, 0.9], axis=1).T,
    columns=['10%', 'median', '90%'])
# 観測値data4と予測値y_pred_dfのデータフレームを結合
y_pred_df = pd.concat([data5, y_pred_df], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
y_pred_df['err_lower'] = y_pred_df['median'] - y_pred_df['10%'] 
y_pred_df['err_upper'] = y_pred_df['90%']  - y_pred_df['median']
# fのCとTでy_pred_dfを分離
y_pred_fC = y_pred_df[y_pred_df['f']=='C']
y_pred_fT = y_pred_df[y_pred_df['f']=='T']

## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# f=Cのy_predの描画（エラーバー付き散布図）
ax.errorbar(y_pred_fC['y'], y_pred_fC['median'],
            yerr=[y_pred_fC['err_lower'], y_pred_fC['err_upper']],
            color='tab:blue', alpha=0.7, marker='o', ms=10, linestyle='none',
            label='C')
# f=Tのy_predの描画（エラーバー付き散布図）
ax.errorbar(y_pred_fT['y'], y_pred_fT['median'],
            yerr=[y_pred_fT['err_lower'], y_pred_fT['err_upper']],
            color='tab:orange', alpha=0.7, marker='^', ms=10, linestyle='none',
            label='T')
# 赤い対角線の描画
ax.plot([-0.5, 8.5], [-0.5, 8.5], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $y$の観測値', ylabel='Predicted: $y$の予測値(中央値)',
       title='$y$ の観測値と予測値（中央値）のプロット 80%区間バー付き')
ax.legend(title='f', bbox_to_anchor=(1, 1))
ax.grid(lw=0.5);

【実行結果】
こちらも微妙な予測値になった感じがいたします・・・
この記事のモデリングに問題があるのか、この図が妥当なのかが、正解の図が公開されていないようなので、確認できません（汗）
確認できました。本モデリングの場合、この予測値になるようです。

上の作図用に作成した一時データフレーム p_pred_df の中身を覗いてみましょう。

### 上の図の作表に利用したy_pred_dfの内容を表示
display(y_pred_df)

【実行結果】
観測値 Y と予測値 median を比べると、結構外している感じです（汗）

MCMCサンプリングデータの散布図行列を描画します。

### MCMCサンプルの散布図行列の描画

## 描画用データの作成
# MCMCサンプリングデータからq1, q100を取り出し
q1_samples = (idata6.posterior['q'].to_dataframe().reset_index()
              .query('data==0').rename({'q': 'q1'}, axis=1))
q100_samples = (idata6.posterior['q'].to_dataframe().reset_index()
               .query('data==99').rename({'q': 'q100'}, axis=1))
# 描画対象パラメータをデータフレーム化
plot_df = pd.concat([param_samples,
                     q1_samples.reset_index(drop=True)['q1'],
                     q100_samples.reset_index(drop=True)['q100']], axis=1)

## 描画処理
# 相関行列プロットの描画
g = sns.pairplot(plot_df, diag_kws={'kde': True, 'ec': 'white'})
# スピアマンの順位相関係数の表示のためのaxフラット化
ax = g.axes.ravel()

## スピアマンの順位相関係数を上三角のaxesに表示
# 列名をリスト化
cols = plot_df.columns
# 列名の組み合わせ行i, 列j ごとにテキスト表示を繰り返す
for i, col1 in enumerate(cols):
    for j, col2 in enumerate(cols):
        # 上三角の位置は 行i < 列j のとき
        if i < j:
            # axesの番号を取得
            pos = i * len(cols) + j
            # スピアマンの順位相関係数を算出
            corr, pval = stats.spearmanr(plot_df[col1], plot_df[col2])
            # テキスト表示：中央表示に関連する引数: x,y,va,ha,transform
            ax[pos].text(x=0.5, y=0.5, s=round(corr * 100), fontsize=30,
                         va='center', ha='center', transform=ax[pos].transAxes,
                         bbox=dict(boxstyle='round', facecolor='white'))

【実行結果】

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

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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

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

お金について考える

36,077件

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