「データ解析のための統計モデリング入門」をPythonで写経 Vol.20 ～ 10章「階層ベイズモデル」①ベイズ統計モデリング～階層ベイズモデル（個体差）

10章「階層ベイズモデル」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」10章「階層ベイズモデル」の Python写経活動記録 です。 

この記事は GLMM の 階層ベイズモデル 化に取り組みます。
二項分布・ロジットリンク関数・ランダム切片（個体差）の GLMM がベイズ統計モデルに！

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

準備

■ 記事の範囲
この記事はテキスト10章の以下の節を取り扱います。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト

Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd
from scipy.special import expit  # ロジスティック関数（シグモイド関数）

# 統計計算
import scipy.stats as stats

# PyMC
import pymc as pm
import arviz as az

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

統計モデリング・サマリー

この記事で扱う統計モデリングの概要です。

■ 統計モデル
GLMM のベイズ統計モデルです。

■ モデリング手続き

データの確認

データを読み込み、データの外観を眺めてから、統計モデルを選択するためのデータの特徴確認を行います。

◼️ データの読み込み
data7a.csv ファイルを pandas データフレームの data に読み込みます。

# データの読み込み
data = pd.read_csv('./data/ch10/data7a.csv')
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 100 です。
100 個体の植物に関する仮想の観測データです。
１個体から８個の種子を調べて生存種子数をカウントしたデータです。
試行回数にあたる調査種子数 $${N}$$ は 8 です。

【変数の説明】
植物の個体 id ごとの生存種子数 y です。

◼️ データの確認
基本的な統計量やチャートでデータを概観します。

① 要約統計量の表示

# 要約統計量
data.describe().round(3)

【実行結果】
生存種子数 y の値は 0 ～ 8 個です。

② 標本分散の表示

# 標本分散
data.var(ddof=1).rename('var').to_frame().T.round(3)

【実行結果】
ｙ の分散はおよそ 9.93 です。

③ ヒストグラムと散布図の描画
seaborn の pairplot() を利用します。
回帰直線を添えます。

# ヒストグラムと回帰直線付き散布図のペアプロットの描画
sns.pairplot(data=data, kind='reg',
             plot_kws={'scatter_kws': {'s': 70, 'alpha': 0.6},
                       'line_kws': {'color': 'tab:red', 'alpha': 0.7}},
             diag_kws={'edgecolor': 'white', 'bins': 5});

【実行結果】
生存種子数 y は Ｕ字型の分布です。
７章で見た「過分散」の分布に似ています。
id と y の間に相関関係は無さそうです。

④ データの分散と二項分布の分散の関係
テキスト p.225 にならって、二項分布における過分散の状況を確認します。

✅ 二項分布の生存確率パラメータ $${q}$$ の最尤推定
テキスト p.173 の生存確率 $${q}$$ の最尤推定量の数式をお借りします。
標本サイズ $${n=100}$$ です。

例題データから生存確率 $${q}$$ の最尤推定量を算出します。

# dataのもとでのqの最尤推定値: q_hat=合計生存種子数/合計調査種子数 p.173の数式を利用
q_hat = data.y.sum() / (8 * len(data))
q_hat

【実行結果】
およそ 0.504 です。

✅ 二項分布の期待値と分散
生存確率 $${q=0.503 \cdots}$$、試行回数 $${N=8}$$ の二項分布の期待値と分散を算出します。

# 生起確率0.504、N=8の二項分布の期待値と分散 p.225
N, p = 8, q_hat
expected_value = N * p
variance = N * p * (1 - p)
print(f'期待値 = {expected_value:.3f}, 分散 = {variance:.3f}')

【実行結果】
分散はおよそ 2.00 です。

✅ 例題データの標本分散
観測値 y の標本分散はさきほど見ました。およそ 9.93 です。
念のため。

# 実際の観測データの分散 ※過分散 p.225
print(f'観測データの分散 = {data.y.var(ddof=1):.3f}')

【実行結果】

【分析】
二項分布の分散 2.00 と比較して観測データの分散は 4 倍以上あり、過分散になっています。

テキスト p225 図 10.1 に相当する観測値と二項分布のチャートを描きます。

# p.225 観測データと二項分布の描画 図10.1(B)

# 観測データの個体数の散布図の描画
plt.plot(data.y.value_counts(), 'o', ms=8, alpha=0.7, label='観測値')

# 生存確率0.504の二項分布に従う生存種子数ごとの個体数の描画
x_val = range(9)
plt.plot(x_val, stats.binom.pmf(k=x_val, n=data.y.max(), p=0.504) * len(data),
         '-o', color='tab:blue', mfc='white', ms=8, 
         label=f'Binom(8, {q_hat:.3f})')

# 修飾
plt.xlabel('生存種子数 $y_i$', fontsize=12)
plt.ylabel('観測された個体数', fontsize=12)
plt.ylim(-1, 29)
plt.legend();

【実行結果】
二項分布だけでは生存種子数を表現できそうにないことが分かりました。

◼️ データの特徴まとめ
データの特徴を整理します。

テキストは「個体由来の原因不明な差異など組み込む GLMM 」を提案しつつ、個体差・場所差などの複数の◯◯差を含めるとパラメータの最尤推定が難しくなるという理由で、ベイズ統計モデルの検討へと進みます。

特徴①より「二項分布」「ロジットリンク関数」とし、特徴②③より「個体差（ランダム切片）」を考慮する GLMM をベイズ統計モデル化します。
今回取り組むモデルは 階層ベイズモデル と呼ばれています。

階層ベイズモデル（個体差）

ベイズ統計モデルの理解

◼️ パラメータの事後分布
ベイズ統計モデリングでは「パラメータの事後分布の推定」を中心に置いて動きます。
ベイズ統計モデルの事後分布は尤度と事前分布の積に比例します。

◼️ 今回モデルの概観
今回のベイズ統計モデリングで推定するパラメータは $${\beta, s, \{r_i\}}$$ です（詳細は後ほど！）。
データを $${\bm Y}$$ とし、事後分布を数式化します。
$${n}$$ は標本サイズであり、例題データの場合は 100 です。

◼️ GLMM の構成要素
個体 $${i}$$ の生存種子数 $${y_i}$$ は二項分布 $${p(y_i \mid 8, q_i)}$$ に従うとします。
生存確率パラメータ $${q_i}$$ は線形予測子とロジットリンク関数を用いて $${\text{logit}(q_i) = \beta + r_i}$$ とします。
$${\beta}$$ は全個体共通のパラメータです。
個体差 $${r_i}$$ は平均パラメータ 0、標準偏差パラメータ $${s}$$ の正規分布 $${\text{Normal}(0, s^2)}$$ に従うとします。

🍀🍀🍀

◼️ 尤度
尤度は次の数式で表されます。

になります。
PyMC のコードに寄せて、次のように表現します。

🍀🍀🍀

◼️ 事前分布
事後分布の数式化では「$${\beta}$$ の事前分布 $${p(\beta)}$$」、「$${s}$$ の事前分布 $${p(s)}$$」、「$${r_i \mid s}$$ の事前分布 $${p(r_i \mid s)}$$ の積 $${\prod_{i=1}^n  p(r_i \mid s)}$$」が示されていました。
３つの事前分布を解きほぐします。

1️⃣ $${\beta}$$ の事前分布
無情報事前分布を指定します。
平均パラメータ 0、標準偏差パラメータ 100 の「すごくひらべったい正規分布」です。

2️⃣ $${r_i}$$ の事前分布
平均パラメータ 0、標準偏差パラメータ $${s}$$ の正規分布を指定します。
階層事前分布に該当します（詳細は後ほど）。

3️⃣ $${s}$$ の事前分布
無情報事前分布を設定します。
ひらべったい正規分布ではなく、幅が十分に広い連続一様分布です。
先出しすると $${s}$$ の MCMC サンプルは区間 $${(1.9, 4.7)}$$ 程度です。

🍀🍀🍀

◼️ 階層事前分布
$${r_i}$$ と $${s}$$ の事前分布の関係を見つめます。

◆ 階層事前分布と階層ベイズモデル
個体差 $${r_i}$$ の事前分布に「平均パラメータ 0、標準偏差パラメータ $${s}$$ の正規分布 $${\text{Normal}(0, s^2)}$$」を指定しています。
そして標準偏差パラメータ $${s}$$ に事前分布を指定しています。
$${r_i}$$ と $${s}$$ の事前分布が階層的に設定されることになります。

事前分布の階層のうち、親（今回の $${p(r_i \mid s)}$$）の方をテキストは 階層事前分布 と呼びます。
階層事前分布に含まれる子のパラメータ $${s}$$ を ハイパーパラメータ（超パラメータ）、子の事前分布 $${p(s)}$$ を ハイパー事前分布（超事前分布）と呼ぶ場合があるそうです。

階層事前分布を含むベイズ統計モデルを 階層ベイズモデル と呼びます。

◆ 事前分布の種類
テキストは 10.4 節で３つの事前分布の種類を説明しています。
テキスト p.233 図 10.5 をお借りして、３つの事前分布のイメージを可視化します。

# ベイズ統計モデルでよく使われる３種類の事前分布の例 p.233 図10.4

## 設定と準備
x_val = np.linspace(0, 1, 201)            # x軸の値
x_min, x_max = x_val.min(), x_val.max()   # xの値の最小値、最大値
scale = x_max - x_min                     # 一様分布のscaleのパラメータ値
s_max = 0.05                              # (C)の正規分布の標準偏差の基準値
titles = ['(A) 主観的な事前分布', '(B) 無情報事前分布', '(C) 階層事前分布']

## 描画
# 描画領域の設定
fig, ax = plt.subplots(1, 3, figsize=(10, 3), sharey=True)
# (A)主観的な事前分布の正規分布の確率密度関数の描画
ax[0].fill_between(x_val, 0, stats.norm.pdf(x_val, loc=0.3, scale=0.07),
                   alpha=0.3)
# (B)無情報事前分布の一様分布の確率密度関数の描画
ax[1].fill_between(x_val, 0, stats.uniform.pdf(x_val, loc=x_min, scale=scale),
                   alpha=0.3)
# (C)階層事前分布の描画
for i in [4, 2, 1]:
    ax[2].fill_between(x_val, 0, stats.norm.pdf(x_val, loc=0.5, scale=s_max * i),
                       color='tab:blue', alpha=0.1)
# (A)～(C)のテキストの表示
ax[0].text(x=0.5, y=4, s='信じる！', fontsize=12)
ax[1].text(x=0.5, y=4, s='わからない？', horizontalalignment='center', fontsize=12)
ax[2].text(x=0.7, y=4, s='$s$ によって\n変わる…', fontsize=12)
# 修飾
for i in range(3):
    # 修飾
    ax[i].set(xticks=np.arange(0, 1.1, 0.2), ylim=(0, 8.5), title=titles[i])
    # 枠線の消去
    ax[i].spines['left'].set_visible(False)
    ax[i].spines['right'].set_visible(False)
    ax[i].spines['top'].set_visible(False)
    # 目盛り・目盛りラベルの消去
    ax[i].tick_params(left=False, labelleft=False)

plt.show()

【実行結果】

(A) 主観的な事前分布は、データ分析者が持つ「先行情報」や「信念」などに基づいて主観的に設定する事前分布です。
テキストは主観的な事前分布を取り扱わないことにしているそうです。

(B) 無情報事前分布は前回記事でも扱ったひらべったい正規分布や一様分布などの「情報をもたせない」事前分布です。
$${\beta}$$ と $${s}$$ に無情報事前分布を設定しています。

(C) 階層事前分布はハイパーパラメータ（今回モデルの $${s}$$）の推定値によって分布の形状が変わります。
$${r_i}$$ に階層事前分布を設定しています。

テキストには、無情報事前分布と階層事前分布のどっちを使うのがいいかのガイダンスが掲載されています。
表 10.1 をお借りします。

$${\beta}$$ はデータ全体を説明する大域的なパラメータですので、無情報事前分布を設定しているようです。

$${r_i}$$ は 個体ごと（y の１個ごと）という、データ全体のうちのごく一部を説明する局所的なパラメータです。
この$${r_i}$$ たちは「ある範囲に分布」すると仮定されます。
この２点を踏まえて、$${r_i}$$ には共通のばらつきパラメータ $${s}$$ を持つ階層事前分布を設定します。
これは $${r_i}$$ のひとつひとつに異なる事前分布を設定しないことを意味します。
データ全体で推定する $${s}$$ を通じて、各 $${r_i}$$ にデータ全体の情報を共有し、かつ、$${r_i}$$ が自由すぎる値を取らないように階層事前分布で拘束して、「ある範囲に分布」する状況が生み出されるようです。

ということで、$${s}$$ は大域的なパラメータであり、無情報事前分布を指定します。

🍀🍀🍀

今回のベイズ統計モデルの数式をまとめます。
$${q_i}$$ は PyMC の設定を考慮して、ロジットリンク関数で包んでいます。

$${\sim}$$ は左側の確率変数が右側の確率分布に従うことを意味します。
$${=}$$ で表現された変数は、等号で結ばれた数式どおりに特定の値が決定される「決定論的変数」です。

パラメータの事後分布の推定

PyMC ライブラリでベイズ統計モデリングを実装します。
テキストの WinBUGS の設定を解釈しつつ、PyMC コードに書き換えていきます。

◼️ 共通設定
調査種子数（試行回数）$${N}$$ と 標本サイズ $${n}$$ を定義します。

# 変数の設定

# 調査種子数（試行回数）
N = 8
# 標本サイズ = 100
n = len(data)

◼️ モデルの定義

# モデルの定義

# coordsの設定
coords = {'id': data.id.values}   # 座標ラベルの設定：データ行の識別子

# モデリング
with pm.Model(coords=coords) as model:
    
    # dataの定義: 目的変数=生存種子数Y
    Y = pm.Data('Y', value=data['y'].values, dims='id')

    # 事前分布
    beta = pm.Normal('beta', mu=0, sigma=100)     # β: 無情報事前分布 N(0,100)
    s = pm.Uniform('s', lower=0, upper=10**4)     # s: 無情報事前分布 U(0, 10^4)
    r = pm.Normal('r', mu=0, sigma=s, dims='id')  # r: 階層事前分布   N(0, s)

    # 線形予測子: 生存確率のロジットlogit_qを用いる
    logit_q = pm.Deterministic('logit_q', beta + r, dims='id')

    # 尤度関数: 試行回数n=8, 生存確率p=qの二項分布
    obs = pm.Binomial('obs', n=8, logit_p=logit_q, observed=Y, dims='id')

【実行結果】なし

🍀🍀🍀

◼️ モデルの確認
モデルの数式と有向グラフを可視化します。
数式を表示します。

# モデルの表示
model

【実行結果】

モデルの有向グラフを描画します。

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

【実行結果】
$${r_i}$$ は個体 id 単位で推定されます。

🍀🍀🍀

◼️ MCMC の実行
MCMC サンプリングを行います。
NUTS サンプラーに nutpie を使います。

%%time
# MCMCサンプリング
# chain=4, draws=1000, tune=1000, thinなし, パラメータ初期値なし, NUTSサンプラー

with model:
    idata = pm.sample(
        random_seed=42,
        nuts_sampler='nutpie',  # nutpieを使わない場合はこの行を削除
    )

【実行結果】
Divergences は０個です。

パラメータ推定値の確認

ざっくり収束の確認などを行います。

◼️ $${\widehat{R}}$$ の確認
$${\widehat{R}}$$ が 1.01 以下になっていることを確認します。
全パラメータの「$${\widehat{R}}$$ ＞1.01」の個数が０になればOKです。

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

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

【実行結果】
すべてのパラメータの $${\widehat{R}}$$ は1.01 以下です。

◼️ 有効サンプルサイズ（ESS）
「実質的に独立なサンプル数」であるESS が 400 以上になっていることを確認します。
全パラメータの「ESS ＜ 400」の個数が０になればOKです。

# 有効サンプルサイズ しきい値を400とした
print((az.ess(idata) < 400).sum())

【実行結果】
すべてのパラメータの ESS は400 以上です。

◼️ トレースプロットの確認

# トレースプロットの表示
var_names = ['beta', 's', 'r', 'logit_q']
pm.plot_trace(idata, var_names=var_names, figsize=(8, 7),
              backend_kwargs={'tight_layout': True});

【実行結果】
右のトレースプロットは、各チェーンが「毛糸玉」のようにゲジゲジと混ざり合い、ドリフトがないのでＯＫ（収束ＯＫ）としましょう。

🍀🍀🍀

◼️ 事後分布の要約統計量
パラメータの事後分布の推定値を確認します。

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

【実行結果】

$${\beta}$$ の事後平均は 0.029、95% HDI は -0.677 ～ 0.715 です。
$${s}$$ の事後平均は 3.039、95% HDI は 2.330 ～ 3.809 です。

ロジット単位の $${\beta}$$ を生存確率 $${q}$$ の単位に変換します。

# βを生存確率qの単位に変換
print(f'logistic(mean)     = {expit(summary.loc['beta', 'mean']):.3f}')
print(f'logistic(hdi 2.5%) = {expit(summary.loc['beta', 'hdi_2.5%']):.3f}')
print(f'logistic(hdi97.5%) = {expit(summary.loc['beta', 'hdi_97.5%']):.3f}')

【実行結果】
生存確率 $${q_i}$$ が $${\beta}$$ だけで構成される場合、平均 0.507、95% HDI 0.337 ～ 0.772 と解釈できます。

🍀🍀🍀

◼️ 事後分布の確認（テキストのチャートの前座）
テキスト p.231 図 10.3 のパラメータの事後分布（密度）を横目に見つつ、PyMC＆ArviZ の標準機能で対応できることをやってみます。
$${\beta, s}$$ の事後分布です。

# 事後分布プロット
pm.plot_posterior(idata, var_names=['beta', 's'], hdi_prob=0.95, figsize=(8, 3),
                  backend_kwargs={'tight_layout': True, 'sharey': True});

【実行結果】

続いて $${r_i}$$ から。先頭の $${r_0, r_1, r_2}$$ の事後分布です。

# 事後分布プロット r の最初の３つ
fig, ax = plt.subplots(1, 3, figsize=(10, 2.5), sharey=True, tight_layout=True)
pm.plot_posterior(
    idata.posterior.r.head({'id': 3}), hdi_prob=0.95, round_to=3, textsize=10,
    ax=ax)
for ax_ in ax.flat:
    ax_.set_xlim(-10, 10)

【実行結果】

◼️ 事後分布の確認（テキストの図 10.3）
テキスト p.231 図 10.3 の形式でパラメータの事後分布を可視化します。
$${\beta, s}$$ の事後分布です。

# MCMCサンプリングによって推定されたβとsの事後分布 p.231 図10.3(A)(B)

## 設定: ２つのグラフの個別設定
vars = ['beta', 's']                            # パラメータ名
xlabels = ['$\\beta$', '$s$']                   # xlabel
heads = ['A', 'B']                              # グラフタイトルのアルファベット
xlims = [(-11, 11), (-1, 7)]                    # x軸の描画範囲
xtickses = [[-10, -5, 0, 5, 10], [0, 2, 4, 6]]  # x軸の軸目盛り

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(7, 3), sharey=True, tight_layout=True)
# for文で使う変数の設定
params = zip(vars, heads, xlabels, xlims, xtickses, axes.flat)
# グラフごとにKDEプロットの描画を繰り返し処理
for i, (var, head, xlabel, xlim, xticks, ax) in enumerate(params):
    # KDEプロットの描画
    sns.kdeplot(idata.posterior[var].stack(sample=('chain', 'draw')), ax=ax)
    # 修飾
    ax.set(xlim=xlim, xticks=xticks,
           title=f'({head}) {xlabel} の事後分布\n全個体共通')
    ax.set_xlabel(xlabel, fontsize=12);

【実行結果】

続いて $${r_i}$$ から。先頭の $${r_1, r_2, r_3}$$ の事後分布です。

# MCMCサンプリングによって推定されたrの先頭3つの事後分布  p.231 図10.3(C)～(E)

## 設定と準備
heads = ['C', 'D', 'E']                      # グラフタイトルのアルファベット
s_median = idata.posterior.s.median().data   # sの事後分布の中央値
x_val = np.linspace(-10, 10)                 # x軸の値
r_prior = stats.norm.pdf(x=x_val, loc=0, scale=s_median)  # 全r_i共通の事前分布

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 2.5), sharey=True, tight_layout=True)
# グラフごとにKDEプロットのと全r_i共通の事前分布の描画を繰り返し処理
for i, (head, ax) in enumerate(zip(heads, axes.flat)):
    # KDEプロットの描画
    sns.kdeplot(idata.posterior.r.stack(sample=('chain', 'draw'))[i], ax=ax)
    # 全r_i共通の事前分布の描画
    ax.fill_between(x_val, 0, r_prior, color='gray', alpha=0.2)
    # 修飾
    ax.set(xlim=(-11, 11), title=f'({head}) $r_{{{i + 1}}}$ の事後分布')
# 修飾2
axes[0].text(x=4, y=0.07, s='全 $r_i$ 共通の\n事前分布')
plt.show()

【実行結果】

グレイの分布は、$${s}$$ の値が事後分布の中央値の場合の事前分布 $${p(r_i \mid s)}$$ です（３チャート共通）。

生存種子数の予測

◼️ 生存種子数 y の予測分布 $${p(y \mid \beta, s)}$$
テキスト p.231 図 10.4 に相当する「生存種子数 $${y}$$ の予測チャート」を描画します。
テキストによると、生存種子数 y の確率分布は、以下の数式で示される「二項分布 $${p(y \mid \beta, r)}$$と正規分布 $${p(r \mid s)}$$ の無限混合分布」です。

ラスボス「積分」が出現しました…
ここでは MCMC サンプルとモンテカルロ積分を用いて、$${r_i}$$ の乱数生成⇒$${q_i}$$の算出⇒ y の確率の推定値 $${\hat{p}}$$ を算出します。
その後、y の平均予測（の個体数）を $${\hat{p} \times n}$$ で算出、y の事後予測（の個体数）を多項分布乱数で生成します。

%%time
# 生存種子数yの予測分布p(y|β,s)と観測データの可視化 p.231 図10.4
# 無限混合分布の積分にモンテカルロ積分を使用

# 描画データの作成

## 設定
x_val = np.arange(N + 1)              # 生存種子数の要素(0～8の整数)
m = 1000                              # モンテカルロ法のサンプルサイズ
rng = np.random.default_rng(seed=42)  # 乱数生成器

## データの準備
# MCMCサンプルからパラメータβとsを取り出し shape=(4000,)
betas = idata.posterior.beta.stack(sample=('chain', 'draw')).data
sigmas = idata.posterior.s.stack(sample=('chain', 'draw')).data

## 生存種子数yの確率 p_hat(y=k) = mean(Binom(k | 8, logistic(β + r)))
# 個体差rの推定: r ~ Normal(0, s), shape=(m, 4000)
rs = stats.norm.rvs(
    loc=0, scale=sigmas, size=(m, len(sigmas)), random_state=rng)
# 生存確率qの推定: logistic(β + r), shape=(4000, m)
qs = expit(betas + rs).T
# yの確率p(y)の推定: p_hat(y=k) ~ mean(Binom(k | 8, q)), shape=(4000, 9)
y_prob = stats.binom.pmf(
    k=x_val[None, None, :], n=N, p=qs[:, :, None]).mean(axis=1)

## 生存種子数yの平均予測
# yの平均予測（個体数）の中央値y_mean_pred_medの算出　⇒描画へ
# p_hat(y) * n, shape=(9,) ※テキストp.230脚注*11参照
y_mean_pred_med = np.median(y_prob * n, axis=0)

## 生存種子数yの事後予測
# yの事後予測（個体数）y_pred_countの算出 shape=(4000, 9)
y_pred_count = np.array(
    [stats.multinomial.rvs(n=n, p=y_prob[i], random_state=rng)
     for i in range(len(y_prob))]
)
# y_pred_countの95%予測区間（2.5%点、97.5%点）の算出 shape=(2, 9)　⇒描画へ
y_pred_count_quantile = np.quantile(y_pred_count, q=[0.025, 0.975], axis=0)

【実行結果】
モンテカルロ法の繰り返し数を 1000 にして、以下の処理時間がかかりました。

算出結果を使ってチャートを描画します。

# 描画

# 描画領域の設定
fig, ax = plt.subplots()

# 観測値の散布図（青い点）の描画 ※Rスクリプトのplot.dataに相当
ax.plot(data.y.value_counts().sort_index(), 'o', ms=8, alpha=0.7, label='観測値')

# yの平均予測の中央値の折れ線グラフ（赤色）の描画
ax.plot(y_mean_pred_med, '-o', ms=10, c='tab:red', mec='white',
        label='予測値（中央値）')

# yの事後予測（個体数）の95%予測区間の塗りつぶし描画
ax.fill_between(x_val, y_pred_count_quantile[0], y_pred_count_quantile[1],
                color='salmon', alpha=0.15, label='予測値（95%区間）')

# 修飾
ax.set_xlabel('生存種子数 $y$', fontsize=12)
ax.set_ylabel('個体数', fontsize=12)
ax.set(ylim=(-1, 28))
ax.legend(bbox_to_anchor=(1.4, 1));

【実行結果】
テキストのチャートによく似ています！

【チャートの補足説明】
赤い折れ線グラフは y の平均予測の中央値です。
薄赤色の塗りつぶしは y の事後予測（ばらつき反映）の 95% 信用区間です。
青い点の y の観測値と平均予測の中央値はまずまず傾向が似ています。
また観測値は事後予測の 95% 信用区間に含まれています。

パラメータの事後分布はうまく推定されているように感じます！

まとめ

今回のベイズ統計モデルをまとめます。

🔷 尤度
観測データは二項分布に従います。

🔷 リンク関数と線形予測子
二項分布のパラメータ $${\text{logit(p)}}$$ は線形予測子のロジット関数です。

🔷 事前分布
線形予測子の切片 $${\beta}$$ の事前分布は、平均０、標準誤差 100 の「すごくひらべったい正規分布」です。

個体差（ランダム切片）$${r_i}$$ の事前分布は、平均０、標準誤差 $${s}$$ 正規分布であり、階層事前分布です。

個体差のばらつきパラメータ $${s}$$ は階層事前分布に関するハイパーパラメータであり、事前分布は $${0}$$ から $${10^4}$$ までの一様分布です。

今回のブログは以上です。

次回は 個体差＋場所差の GLMM をベイズ統計モデル化 します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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