「データ解析のための統計モデリング入門」をPythonで写経 Vol.14 ～ ７章「一般化線形混合モデル（GLMM）」①一般化線形混合モデルの骨格～無限混合分布

７章「一般化線形混合モデル（GLMM）」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」７章「一般化線形混合モデル（GLMM）」の Python写経活動記録 です。 

書籍は第７章で 一般化線形混合モデル（GLMM） をやり切ります！
非常に濃い章ですので、記事を分割します。

この記事は１．の 一般化線形混合モデルの骨格 を学びます。
過分散、ランダム切片、無限混合分布 が主要テーマです。
統計モデルの構築は次回以降に行う予定です。

長文になります。
マイペースでお読みくださいな🍀

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

はじめに

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

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

準備

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd
from scipy.special import expit  # ロジスティック関数

# 統計計算
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

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

イントロダクション

一般化線形混合モデル（GLMM）とは

一般化線形混合モデル（Generalized Linear Mixed Model：GLMM）を一言で表してみます。

前回までの記事は GLM：一般化線形モデルに焦点を当てました。
７章では「混合」が加わった GLMM を学びます。
GLM に混合するのは「ランダム効果」です。
上の文章のランダム切片、ランダム傾きがランダム効果を表現しています。

ちなみにランダム効果を除くと、残りは「固定効果」です。
固定効果とランダム効果が混合しているモデルなので一般化線形混合モデルなのです。

７章で取り組む GLMM

７章では、二項分布とロジットリンク関数の GLM に「ランダム切片 $${r_i}$$」を加えた GLMM に取り組みます。
この GLM と GLMM の線形予測子を見ておきます。

📊 二項分布とロジットリンク関数の GLM の線形予測子

📊 上の GLM にランダム切片を追加した GLMM の線形予測子

固定効果とランダム効果を区別するとこんな感じに。

ランダム切片の正体

ランダム切片 $${r_i}$$ をよく見ると $${i}$$ が付いています。
$${\beta_1, \beta_2}$$ には無いインデックス $${i}$$。
これは、$${r_i}$$ がデータ１行ごとに存在することを表しています。
７章の例題データの場合、植物の個体 $${i}$$ ごとに $${r_i}$$ が現れます。

しかもこの $${r_i}$$ はデータ内の変数に含まれていません。
観測されないのです！
つまりGLMM モデルの中にのみ現れる「潜在的な変数」らしいのです。
ますます謎は深まります…

テキストによると、今回のモデルの $${r_i}$$ は、個体ごとの「ばらつき」を表す「個体差」です。
ばらつきと言えば「確率分布」ですね～。
テキストは次のように教えてくれます。

個体差 $${r_i}$$ が正規分布に従う様子を数式表現します。

ちょっと待って…
この GLMM は確か、二項分布とロジットリンク関数の GLM を拡張したものでした。
観測値 $${y_i}$$ は「二項分布」に従っているということですよね…
ということは…

この GLMM には、
・$${y_i}$$ のばらつきを表現する「二項分布」
・個体差のばらつきを表現する「正規分布」
の２つの確率分布が混載されている、ということですよね！

ややこしいですね…（汗）

GLMM のパラメータ推定～さらなる謎へ

いったん落ち着きましょう。

GLM のモデリングで「最尤推定」したのは線形予測子に含まれるパラメータ（係数）でした。
たとえば、$${\beta_1, \beta_2}$$ です。
またひっそりと、ばらつきを表す分散パラメータ（dispersion parameter） の $${\phi}$$ を最尤推定の外で推定しました。

では GLMM で推定するパラメータは何でしょう？
$${\beta_1, \beta_2}$$ に加えて、$${r_i}$$ も？、$${s}$$ も？

実は…
固定効果のパラメータ $${\beta_1, \beta_2}$$ とランダム効果の標準偏差パラメータ $${s}$$（または、分散パラメータ $${s^2}$$）を最尤推定で推定します。
なんと、$${r_i}$$ は最尤推定に含まれないのです。
（じゃあどうするの？は次回記事で！）

ここまでの仮まとめ

■ モデルの概要（部分）
次の線形予測子を持つ GLMM を７章の例題データに当てはめます。

■ ランダム切片
特徴的なのはランダム切片 $${r_i}$$ です。
$${r_i}$$ は植物の個体ごとのばらつき＝個体差を表す潜在変数であり、平均０、分散 $${s^2}$$ の正規分布に従います。

■ パラメータ推定
最尤推定で推定するパラメータは $${\beta_1, \beta_2, s\ (またはs^2)}$$ です。

GLMM のイントロを頭の片隅に置いて、テキストによる GLMM の学びを進めていきます！

過分散が見られるデータ

データの確認

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

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

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

【実行結果】
データの個数（標本サイズ）は 100 です。
植物の個体（個体数 100）に関する仮想の観測データです。
葉数 x の値ごとに 20 の個体を調査しています。

【変数の説明】
植物の個体ごとに採取した調査種子数 N、生存種子数 y（目的変数）、葉数 x です。

◆ ◆ ◆

■ 要約統計量の確認
ざっくりデータの特徴を見ます。

① 要約統計量の確認

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

【実行結果】

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

【実行結果】
生存種子数 y の分散 9.428 が大きいことに注目です。

◆ ◆ ◆

■ 生存種子数 y の可視化
y のヒストグラムを描画します。

# yのヒストグラムの描画

# 描画領域の設定
plt.figure(figsize=(5, 4))
# ヒストグラムの描画
sns.histplot(data=data, x='y', bins=np.arange(-0.5, 8.6), edgecolor='white',
             alpha=0.7)
# 修飾
plt.xlabel('生存種子数 $y_i$', fontsize=12)
plt.ylabel('頻度（個体数）', fontsize=12);

【実行結果】
Ｕ字型の分布になっています。

葉数別にヒストグラムもどきを描画します。

# 葉数別・生存種子数別の個体数

# 葉数別・生存種子数別の個体数をカウント
data_groups = data.groupby(['x', 'y'])['y'].count().rename('freq').reset_index()
# 描画領域の設定
plt.figure(figsize=(5, 4))
# 葉数別・生存種子数別の個体数の折れ線グラフの描画
sns.lineplot(data=data_groups, x='y', y='freq', marker='o', ms=7, lw=2,
             hue='x', palette='Set1')
# 修飾
plt.xlabel('生存種子数 $y_i$', fontsize=12)
plt.ylabel('頻度（個体数）', fontsize=12)
plt.xticks(range(0, 9))
plt.yticks(range(0, 10))
plt.legend(title='葉数 $x_i$', bbox_to_anchor=(1, 1));

【実行結果】
葉数が小さいと生存種子数が小さい傾向、葉数が大きいと生存種子数が大きい傾向が分かります。

◆ ◆ ◆

■ 葉数 x と生存種子数 y の可視化
散布図と生存種子数の平均値プロットを描画します。

# x,yの散布図と平均値プロット

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
# 観測値の平均値の折れ線グラフの描画
sns.lineplot(data=data, x='x', y='y', label='平均値', ax=ax)

# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存種子数 $y_i$', fontsize=12);

【実行結果】

オレンジの帯は生存種子数の平均値の 95% 信頼区間です。
葉数が大きくなるにつれて生存種子数が大きくなる傾向が分かります。
ただし、観測値のばらつき（縦方向）が大きい印象があり、平均値と 95% 信頼区間が各データ点をカバーできないようです。
このばらつきが「個体差」（個体のばらつき）なのでしょうか？

◆ ◆ ◆

■ データの散布図と真の生存確率に基づく予測値の重ね描き
テキスト p.146 の図 7.2 に相当する散布図を描画します。
「真の生存確率」に基づく生存種子数の予測値的な曲線を重ねます。
真のパラメータ値は $${\beta_1=-4,\ \beta_2=1}$$ です。

# GLMではうまくあつかえない生存種子数の例題 p.146 図7.2(B)

## 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# x軸の値
x_val = np.linspace(data.x.min(), data.x.max(), 100)
# 真の生存確率の一例（点線）のy軸の値（ロジスティック関数）
y_val = expit(beta1 + beta2 * x_val) * data.y.max()

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
# 真の生存確率の一例の点線の描画
ax.plot(x_val, y_val, color='tab:red', ls='--', label='真の生存確率の例')
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存種子数 $y_i$', fontsize=12);
ax.set_title('(B) 全100個体の $x_i$ と $y_i$')
ax.legend(loc='upper left');

【実行結果】
真の生存確率に基づく曲線は、上の平均値プロットよりも「傾きの大きい曲線」になっています。
この曲線は個体のばらつきを考慮していると考えられます。
ただ、曲線を見ただけでは、個体のばらつきとの関係を読み取れません。
（モヤッとします…）

ロジスティック回帰が合わない

テキスト p.147 に従って、GLM であるロジスティック回帰を例題データに当てはめます。

■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family で二項分布 Binomial() を指定します。
リンク関数はデフォルトのロジット関数が適用されます。
当てはめ結果を変数 result_binom に格納します。

# 二項分布とロジットリンク関数を組み合わせたモデルの当てはめ p.147

# 設定
family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布

# モデルの当てはめ ※目的変数には生存数yと死滅数N-yを与える
result_binom = smf.glm(formula='y + I(N-y) ~ x', data=data, family=family).fit()

# 結果表示
result_binom.summary()

【実行結果】
特に問題は無さそうです。

■ GLM で説明しきれない観測データ
ロジスティック回帰や二項分布がこのデータに合わないことを、テキスト p.147 図 7.3 をお借りして確認します。

# ロジスティック回帰がうまくいかない例題 p.147 図7.3

## 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# 調査種子数
N = 8
# x軸の値
x_val = np.linspace(data.x.min(), data.x.max(), 100)
# 真の生存確率の一例（点線）のy軸の値（ロジスティック関数）
y_val = expit(beta1 + beta2 * x_val) * N

## 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

## 左の(A)の描画
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax1)
# 真の生存確率の一例の点線の描画
ax1.plot(x_val, y_val, color='tab:red', ls='--')
# GLM(ロジスティック回帰)の予測値の緑実線の描画
ax1.plot(x_val, result_binom.predict(dict(x=x_val)) * N, color='green', lw=2)
# 修飾
ax1.set_xlabel('葉数 $x_i$', fontsize=12)
ax1.set_ylabel('生存種子数 $y_i$', fontsize=12);
ax1.set_title('(A) 全100個体の $x_i$ と $y_i$')

## 右の(B)の描画
# 設定
xi = 4                                         # xiの値=4
num_xi = len(data[data.x==xi])                 # x=4のデータ個数
plot_data = data[data.x==xi]['y'].value_counts().to_frame() # 生存種子数別個体数
k = list(range(N+1))                           # x軸の値
p = result_binom.predict(dict(x=[xi])).values  # x=4のときの生存確率pの推定値:0.47
# x=4のときの生存種子数yの観測された個体数の描画
sns.scatterplot(data=plot_data, x='y', y='count', s=60, alpha=0.7, ax=ax2)
# x=4のときのパラメータpの推定値を用いた二項分布の確率質量関数の描画
ax2.plot(k, stats.binom.pmf(k=k, n=N, p=p) * num_xi, '-o', ms=7, c='green',
         mec='white')
# 修飾
ax2.set_xlabel('生存種子数 $y_i$', fontsize=12)
ax2.set_ylabel('観測された個体数', fontsize=12);
ax2.set_title(f'(B) $x_i=${xi} での種子数分布')
ax2.set_ylim(top=6.2);

【実行結果】

【左のチャートの考察】
緑実線がロジスティック回帰の予測値の曲線（ほぼ直線）です。
真の生存確率に基づく赤点線よりも「葉数に対する傾きの値が小さく」なってしまっています。
ロジスティック回帰はパラメータを適切に推定できなかったようです。

【右のチャートの考察】
葉数＝４のデータに特化して、生存種子数の個体数（度数）の分布を描画しています。
青い点が観測値であり、Ｖ字型になっています。
緑実線がロジスティック回帰で推定したパラメータを用いて、生存確率 $${q_i = \text{logistic}(-2.15 + 0.51 \times 4) \approx 0.473 }$$ の二項分布を描いています。
この二項分布は観測値の分布と合っていません。
生存種子数 y が二項分布に従っていないように見えます。

テキストは「二項分布で期待されるよりも大きなばらつきが原因で、ロジスティック回帰の推定値が不正確になった」旨を記しています。

仮定する確率分布で期待される分散よりもデータの分散が大きいことを 過分散 と呼ぶようです。

過分散を深堀り

テキストにならって、二項分布を仮定すると過分散になる状況を数字で追いかけます。
葉数＝４のデータを対象にして、データの分散と二項分布の分散を比較します。

1️⃣ 準備：葉数４のデータを抽出

# x=4のサブセットd4を作る。あわせてデータの分布を表示する p.148

# 葉数x==4のデータを抽出
data_d4 = data[data['x']==4]

# 葉数==4データのyの分布
data_d4['y'].value_counts().sort_index().to_frame().T

【実行結果】
生存種子数 y は０個～８個までばらついている感じがします。

2️⃣ データの平均・分散を調べる

# d4の平均と分散を調べる p.148~149
mean4 = data_d4['y'].mean()
var4 = data_d4['y'].var()
print(f'd4データの平均: {mean4}\nd4データの分散: {var4}')

【実行結果】
データの平均・分散を使って、次の二項分布の分散を算出します。

3️⃣ 二項分布の分散を調べる
二項分布のパラメータを計算して、分散 $${N p (1-p)}$$ を計算します。
・調査種子数（試行回数）パラメータ：$${N=8}$$
・生存確率（成功確率）パラメータ：$${p=データの平均 / 8}$$

# d4の生存割合（生存確率）と二項分布の分散 N*p*(1-p) の計算 p.149
p4 = mean4 / 8
var4_binom = 8 * p4 * (1 - p4)
print(f'd4データの生存確率の平均: {p4}\nd4データの二項分布の分散: {var4_binom}')

【実行結果】

データの分散 $${8.37}$$は二項分布の分散 $${2.00}$$ の４倍以上です。
テキストは次のように過分散の状況をまとめています。

◆ ◆ ◆

⏰️ アディショナルタイム！ ⏰️
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。

せっかくなので葉数４以外のデータも確かめましょう。
データの分散と二項分布の分散を比較する関数を作ります。

# データの分散と二項分布の分散を比較する関数の定義
# 引数 data: データ(array-like), N: 二項分布の試行回数パラメータ(int)

def binom_var_checker(data, N):

    ## 準備
    # データをNumpy配列化
    data = np.array(data)
    
    ## データの分散を算出 ※不偏分散
    var_data = data.var(ddof=1)

    ## 二項分布の分散の算出
    # データを二項分布に当てはめてパラメータ推定する ※bounds: N, p
    result = stats.fit(stats.binom, data, bounds=[(N, N), (0, 1)])
    # 推定した成功確率pを取得
    p = result.params.p
    # 二項分布の分散を算出 ※var = N * p * (1-p)
    var_binom = stats.binom.var(n=N, p=p)
    
    ## 結果を比較表示
    print(f'データの分散　: {var_data}\n二項分布の分散: {var_binom}')

【実行結果】なし

最初に全データの分散を確かめます。

# 全データの分散の比較
N = 8
print('全データ:')
binom_var_checker(data['y'], N)

【実行結果】
データの分散は二項分布の分散よりも大きい過分散になっています！

続いて葉数ごとの分散を確かめます。

# 各葉数ごとの分散の比較
for i in range(data.x.min(), data.x.max()+1):
    print(f'葉数x={i}:')
    binom_var_checker(data[data['x'] == i]['y'], N)
    print()

【実行結果】
すべての葉数で、データの分散が二項分布の分散よりも大きい過分散になっています！

今回のデータは、二項分布の想定よりもばらつきが大きく、過分散になっていることが分かりました。

いよいよ真打ち・個体差「$${r_i}$$」の登場です！
$${r_i}$$ はランダム切片です！

GLMMへ橋渡し：ランダム切片の導入

ランダム切片・再考

■ ランダム切片と個体差

まずは、ChatGPTとGeminiがまとめてくれたランダム切片をご覧ください。

「グループ」は、テキストが扱う個体ごとだったり、ある群単位だったりします。
テキストは生物の生態学をベースにしていて、グループごとの差の例に「個体差」「場所差」「ブロック差」「植木鉢差」を挙げています。

個体差に関して、テキスト p.150～ は「生物の個体差をもたらしうる要因は少なくとも２種類ある」とし、次の２つを説明しています。

• 生物的な要因
  例：個体の遺伝子、年齢、経験履歴の相違

• 非生物的な要因
  例：栄養・水分・光などの生育環境の相違

■ ランダム切片の例
生物の個体差以外の「グループの差」の例をChatGPTに教えてもらいました。

個体差の導入

■ ランダム切片：個体差の導入
テキストは、個体差や場所差の影響を「原因不明のまま取り込む」統計モデルが必要だとして、GLMM に進みます。

📊 個体差（ランダム切片）を追加した GLMM の線形予測子
生存確率 $${q_i}$$ に関する GLMM の線形予測子をあらためて掲載します。

■ 個体差で生存確率が変わることを可視化で確認
個体差 $${r_i}$$ の値によって、生存確率 $${q_i}$$ が変わることを可視化します。
テキスト p.152 図 7.5 に相当します。

真のパラメータ $${\beta_1 = -4,\ \beta_2 = 1}$$ を用いて、$${r_i = \{ -1.5,\ 0,\ 1.5 \}}$$ のケースを見ます。
テキストによると $${r_i = 0}$$ は平均的な個体に相当し、GLM の線形予測子、つまり個体差を考慮しない線形予測子に該当します。

# 個体差 r_i と生存確率 q_i p.153 図7.5
# logit(q_i) = beta1 + beta2 x_i + r_i

## 設定 ※beta1, beta2, rsの設定内容はテキストと異なる可能性があります
beta1, beta2 = -4, 1                           # 線形予測子の係数
rs = [1.5, 0, -1.5]                            # 線形予測子のr
colors = ['tab:green', 'tab:blue', 'tab:red']  # 線の色
labels = ['$r_i>0$', '$r_i=0$', '$r_i<0$']     # ラベル
x_val = np.linspace(1, 7, 100)                 # x軸の値

## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# rの値ごとに線形予測子のロジスティック関数の曲線を描画
for r, color, label in zip(rs, colors, labels):
    # 曲線を描画
    ax.plot(x_val, expit(beta1 + beta2 * x_val + r), color=color, label=label)
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存確率 $q_i$', fontsize=12)
ax.set(xlim=(2, 6), ylim=(0, 1))
ax.legend();

【実行結果】

個体差 $${r_i}$$ によって生存確率 $${q_i}$$ が変わることが分かりました。
生存確率 $${q_i}$$ が変わると、もちろん、生存種子数 $${y_i}$$ も変わってきます！

個体差のばらつきと過分散

■ ランダム切片：個体差 $${r_i}$$ が従う確率分布

個体差 $${r_i}$$ が平均０、標準偏差 $${s}$$（分散 $${s^2}$$）の正規分布に従うと仮定します。

📊 個体差 $${r_i}$$ のばらつきを表す確率分布（再掲）

また、個体差 $${r_i}$$ は個体間で独立（独立同分布 $${\text{i.i.d.}}$$）と仮定すると、確率密度関数 $${f(r_i \mid \mu=0,\ \sigma=s)}$$ は次のようになります。

■ 個体差 $${r_i}$$ が従う正規分布の可視化
平均０、標準偏差 $${s}$$ の正規分布の確率密度関数 $${f(r_i \mid \mu=0,\ \sigma=s)}$$ を可視化します。
テキスト p.153 図 7.6 に相当します。
標準偏差 $${s}$$ のパターンはテキストの $${1.0,\ 1.5,\ 3.0}$$ に、図 7.7 に現れる $${0.5}$$ を加えます。

# 個体差をあらわす平均ゼロの正規分布の確率密度関数 f(r_i|s) p.153 図7.6

## 設定
sigmas = [0.5, 1, 1.5, 3]        # 標準偏差s
x_val = np.linspace(-8, 8, 100)  # x軸の値

## 描画
# sigmaの値ごとに正規分布の確率密度関数の描画を繰り返し処理
for s in sigmas:
    # μ=0, σ=sの正規分布の確率密度関数 f(r_i|μ=0, σ=s) の曲線を描画
    plt.plot(x_val, stats.norm.pdf(x_val, loc=0, scale=s), label=f'$s=${s:.1f}')
# x=0の垂直線の描画
plt.axvline(0, color='black', ls='--')
# 修飾
plt.xlabel('個体差 $r_i$', fontsize=12)
plt.ylabel('確率密度 $f\ (r_i \mid s)$', fontsize=12)
plt.xlim(-7.9, 7.9)
plt.legend();

【実行結果】

平均０の正規分布の確率密度関数なので、平均０の確率密度が大きくなっています。
また標準偏差 $${s}$$ の値が大きくなるにつれて、峰が潰れて裾が厚い形状になります。
テキストは確率密度と標準偏差 $${s}$$ を次のように説明しています。

◆ ◆ ◆

■ 標準偏差 $${s}$$ と過分散の関係
テキスト p.153 図 7.7 の図を描画して、「個体差 $${r_i}$$ のばらつきの大小」と「過分散による二項分布からの乖離」の関係を直感的に確認します。
例題データを使わずに、正規分布乱数⇒二項分布乱数を生成しています。
生存確率は $${q_i = \text{logistic}(r_i)}$$ で算出しています。

# 個体差r_iのばらつきの大きさsと過分散の関係 p.153 図7.7
# 線形予測子はbeta1=beta2=0を仮定し、ばらつきの大きさr_iのみとなる
# テキストと乱数生成が異なるため、結果はテキストと異なります

## 設定
n = 50                                         # サンプルサイズ
sigmas = [0.5, 3]                              # 個体差のばらつきの大きさ s
titles1 = ['A', 'B']                           # グラフタイトルのアルファベット
titles2 = ['小さい', '大きい']                  # グラフタイトルの大小
x1_val = np.linspace(-8, 8, 100)               # １行目の確率密度関数のx軸の値
x2_val = list(range(9))                        # ２行目の生存種子数のx軸の値
rng = np.random.default_rng(seed=19)           # 乱数生成器 0 17 19 33

## 描画
# 描画領域の設定
fig, ax = plt.subplots(2, 2, figsize=(10, 8), tight_layout=True)

# ばらつきsごとにのグラフ描画を繰り返し処理
for i, (s, title1, title2) in enumerate(zip(sigmas, titles1, titles2)):

    ## １行目の設定と準備
    # p(r_i|s)が生成した50個体分の{r_i}（正規分布乱数）
    rs = stats.norm.rvs(loc=0, scale=s, size=n, random_state=rng)
    
    ## １行目のr_iの分布の描画
    # 平均0、標準偏差sの正規分布の確率密度関数の描画
    ax[0, i].plot(x1_val, stats.norm.pdf(x=x1_val, loc=0, scale=s),
                  label=f'$s=${s:.1f}')
    # r_iのラグ（赤い線）の描画
    sns.rugplot(x=rs, ax=ax[0, i], height=0.07, color='tomato', label='観測値')
    # x=0の垂直線の描画
    ax[0, i].axvline(0, color='black', ls='--', lw=1)
    # 修飾
    ax[0, i].set_xlabel('個体差 $r_i$', fontsize=12)
    ax[0, i].set_ylim(-0.2, 16.2)
    ax[0, i].set(xlim=(-7.9, 7.9), ylim=(0, 0.85),
                 title=f'({title1}) 個体差のばらつきが{title2}場合')
    ax[0, i].legend() 

    ## ２行目の設定と準備
    # 確率 q_i=1/(1+exp(-r_i))の二項乱数で生存種子数y_iを生成
    y = stats.binom.rvs(size=n, n=max(x2_val), p=1 / (1 + np.exp(-rs)),
                        random_state=rng)
    # 生存種子数ごとの個体数のカウント
    num_y, counts = np.unique(y, return_counts=True)
    
    ## ２行目の散布図と二項分布の描画
    # 観察された個体数の散布図の描画 ※標本分散には不偏分散を記載
    ax[1, i].plot(num_y, counts, 'o', ms=10, color='tomato', mec='red',
                  alpha=0.5, label=f'個体数:標本分散 {np.var(y, ddof=1):.1f}')
    # 試行回数8, 生存確率0.5の二項分布の確率質量関数の描画
    ax[1, i].plot(x2_val, stats.binom.pmf(k=x2_val, n=max(x2_val), p=0.5) * n,
                  '-o', ms=8, color='tab:blue', mec='white', alpha=0.7,
                  label=f'二項分布\np=0.5, 分散{8*0.5*0.5:.1f}')
    # 修飾
    ax[1, i].set_xlabel('生存種子数 $y_i$', fontsize=12)
    ax[1, i].set_ylim(-0.2, 16.2)
    ax[1, i].legend()

ax[0, 0].set_ylabel('確率密度', fontsize=12);
ax[1, 0].set_ylabel('観察された個体数', fontsize=12);

【実行結果】
標準偏差 $${s}$$ が大きい＝個体差 $${r_i}$$ のばらつきが大きいとデータの分散が二項分布の分散より大きくなる「過分散」、そして、過分散によってデータの分布が二項分布と乖離する様子を実に巧みに表現しています。
著者の先生、凄い！

個体差のばらつきが大きいと過分散になり、データの分布と二項分布が乖離していく様子が分かりました。

二項分布で表せない生存種子数 $${y_i}$$。
ではどうすれば、ばらつきのある生存種子数 $${y_i}$$ を表現できるのでしょう？

テキストは 二項分布と正規分布を混ぜ合わせた分布 で表現する一般化線形混合モデルへと進みます！

GLMMと無限混合分布

条件付き分布と無限混合分布

生存種子数 $${y_i}$$ はどんな確率分布に従うのでしょう？
ChatGPT が教えてくれたストーリーで確率分布を探しに行きます！

◆ ◆ ◆ 

個体差 $${r_i}$$ で条件付けた生存種子数 $${y_i}$$ は二項分布に従う、からリスタートです。

1️⃣ 条件付き分布の世界
個体差 $${r_i}$$ が与えられているとき、生存種子数 $${y_i}$$ は二項分布に従います。

ここで、

つまり、特定の $${r_i}$$ のもとでは単純な二項分布モデルです。
この段階では $${r_i}$$ は固定された条件として扱います。

2️⃣ ランダム効果の分布
実際には $${r_i}$$ は観測されない潜在変数で、平均 0、分散 $${s^2}$$ の正規分布に従うと仮定します。

3️⃣ 最尤推定の課題
もし $${r_i}$$ をパラメータとして同時推定すると、個体数分の未知パラメータが増え、通常の最尤法では不安定になります（パラメータ爆発問題）。

4️⃣ 周辺化でランダム効果を消す
そこで、$${r_i}$$ はパラメータではなく確率変数として扱い、その分布に基づき積分で消去します。

こうして尤度は、固定効果と分散パラメータだけの関数になります。

5️⃣ 無限混合分布
この積分は「$${r_i}$$ の取り得るすべての値について二項分布を混ぜ合わせる」ことを意味します。

$${r_i}$$ の値は連続的に無限にあるため、これは 無限混合分布 です。

6️⃣ まとめ

• 条件付き分布：$${r_i}$$ が既知なら二項分布

• 推定方法：$${r_i}$$ を確率変数として周辺化

• 結果：二項分布と正規分布の 無限混合分布 が生存種子数の分布になる

生存種子数は「二項分布と正規分布の無限混合分布」に従うことが ChatGPT から明かされました。
この無限混合分布のパラメータは $${\beta_1, \beta_2, s}$$ であり、個体差 $${r_i}$$ は消えてしまいました…

無限混合分布の尤度

無限混合分布をもとにした「個体 $${i}$$ ごとの尤度関数」はテキスト p.156 に掲載された尤度 $${L_i}$$ になります。

個体ごとの尤度からも個体差 $${r_i}$$ が消えます。

全体の尤度関数 $${L}$$ は個体ごとの尤度関数 $${L_i}$$ の積です。

二項分布と正規分布の無限混合分布

生存種子数が従う無限混合分布も、尤度関数も、二項分布と正規分布を混ぜ合わせています。
テキストによると「二項分布の確率を正規分布の確率密度 $${f(r_i \mid 0, s^2)}$$ で重み付けして」２つの確率分布を混合しています。

「二項分布の確率 $${\times}$$ 正規分布の確率密度」をテキストに準拠して可視化します。
テキスト p.157 図 7.8 に相当します。

この図例は以下の条件で作成されます。

個体差 $${r=\{ -2.20,\ -0.60,\ 1.00,\ 2.60 \}}$$ をチョイスして、次のチャートを描きます。

# 「分布を混ぜる」という考え方Ⅰ 二項分布と正規分布の無限混合分布の例 p.157 図7.8

## 設定
N = 8                                          # 生存種子数の最大値
sigma = 3                                      # rの分散s
rs = [-2.20, -0.60, 1.00, 2.60]                # rのリスト
y_val_bin = list(range(N + 1))                 # 二項分布の横軸のyの値
r_val = np.linspace(-sigma, sigma, 100) * 2.7  # 正規分布の横軸のrの値
norm_dist = stats.norm(loc=0, scale=sigma)     # 平均0,標準偏差3の正規分布の作成

## 描画
# 描画領域の設定
fig, ax = plt.subplots(4, 2, figsize=(6, 10), tight_layout=True)
# rの値ごとに二項分布（左）と正規分布（右）の描画を繰り返し処理
for i, r in enumerate(rs):
    ## 設定
    q = expit(r)                # 二項分布の確率の算出 q=1/(1+exp(-r))
    pdf_r = norm_dist.pdf(x=r)  # r=rのときの正規分布の確率密度関数の算出
    
    ## 個体差rごとに異なる二項分布の確率質量関数の描画
    # 二項分布の確率質量関数の描画
    ax[i, 0].plot(y_val_bin, stats.binom.pmf(k=y_val_bin, n=N, p=q), '-o', ms=8, 
                  mec='white')
    # 修飾
    ax[i, 0].set(ylim=(-0.03, 0.6), xlabel='$y$',
                 title=f'$r=${r:.2f},  $q=${q:.2f}')
    
    ## 集団内のrの分布：重みf(r|s)の正規分布の確率密度関数の描画
    # rの垂直点線の描画
    ax[i, 1].vlines(r, -0.03, pdf_r, color='black', ls='--', lw=1)
    # 正規分布の確率密度関数の描画
    ax[i, 1].plot(r_val, norm_dist.pdf(x=r_val))
    # f(r|s)の白抜き点の描画
    ax[i, 1].plot(r, pdf_r, 'o', ms=8, mec='tab:blue', mfc='white')
    # 修飾
    ax[i, 1].set(ylim=(0, 0.15), xlabel='$r$', title=f'$f(r)=${pdf_r:.2f}')

【実行結果】
「左の二項分布の確率 $${\times}$$ 右の正規分布の確率密度」は、無限混合分布の中の特定の $${r}$$ をピックアップしたものになります。

続いて、すべての個体差 $${r}$$ に関する積分結果の可視化です。
テキストでは「集団をあらわす混合された分布」と紹介されています。

## （続き）描画２：集団全体をあらわす混合された分布

## 設定と準備
# 正規分布に従うrの200個のパーセンタイル点の取得
rs2 = norm_dist.ppf(q=np.linspace(0.001, 0.999, 200))
# 上記rで生存確率logistic(r)の二項分布に従う生存種子数yの確率の取得 shape=(200, 9)
y_mix_bin = np.array(
    [stats.binom.pmf(k=y_val_bin, n=N, p=expit(r)) for r in rs2]
)

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
# 混合された分布の描画 ※y_mix_bin.mean()は等確率分位点の等重み平均
ax.plot(y_val_bin, y_mix_bin.mean(axis=0), '-o', ms=8, mfc='white')
# 修飾
ax.set_xlabel('$y$', fontsize=12)
ax.set_ylabel('確率', fontsize=12)
ax.set(ylim=(0, 0.25));

【実行結果】
$${y}$$ の確率がＵ字形になっています。
データの確認時に描画した生存種子数のヒストグラムもＵ字形でした。
二項分布と正規分布の無限混合分布が例題データを表現できそうな気がします！

この２つの図の解説をChatGPTにまとめてもらいました。

このコードがやっていることはシンプルで「条件付きの二項分布を、正規分布で“混ぜて”周辺分布（無限混合分布）を作る」デモです。
ポイントを段階的にほどきます。

📈 前半：条件付きの世界を見せる（r を固定した二項分布）

rs = [-2.20, -0.60, 1.00, 2.60]
q = expit(r)                # p = logistic(r)
stats.binom.pmf(k=y_val_bin, n=N, p=q)

• ここでは個体差 r を固定（-2.20 など）して、$${q=\text{logistic}(r)}$$ の 二項分布の 確率質量（pmf） を描いています（左列のグラフ）。

• 右列はその r の位置における正規分布の密度 $${f(r\mid 0,\sigma^2)}$$ を示して、「混ぜるときの重みは正規密度だよ」と可視化しています。

📈 後半：無限混合分布（周辺 pmf）を近似で作る

# 正規の 0.1%〜99.9% の 200 個のパーセンタイル点
rs2 = norm_dist.ppf(q=np.linspace(0.001, 0.999, 200))
# 各 r での二項 pmf を行列に（200 × 9）
y_mix_bin = np.array([stats.binom.pmf(k=y_val_bin, n=N, p=expit(r)) for r in rs2])
# 列ごと平均をとる → 混合（周辺 pmf）の近似
ax.plot(y_val, y_mix_bin.mean(axis=0), ...)

■ この図は何を示しているの？
本来の周辺 pmf は積分で：

これを 数値的に近似しています。
ここでは 正規分布の分位点（ppf） を 200 個とってきて、

• 各 $${r}$$ ごとに条件付き pmf を計算

• それらを 等重み平均（$${\mathtt{mean(axis=0)}}$$）
  → 無限混合（積分）を有限個の代表点の平均で置き換えるという発想

■ なんで等重みでいいの？

• $${\mathtt{np.linspace(0.001, 0.999, 200)}}$$ は 確率空間で等間隔の点（= 各区間の確率がほぼ同じ）をとっています。
  つまり 「等確率の層」から1点ずつサンプルして平均しているイメージで、
  層化モンテカルロのような効果で安定した近似になります。

• 厳密な「数値積分としての重み付き和」にしたいなら Gauss–Hermite を使います（下に補足）。

📈 数式でつなぐ（コードがやっていることの芯）

• 条件付き：$${Y\mid r \sim \text{Binomial}(N, \text{logistic}(r))}$$

• ランダム効果：$${r\sim \text{Normal}(0,\sigma^2)}$$

• 周辺（無限混合）：

ここで $${r_k}$$ は正規分布の分位点からとった代表値（このコードの $${\mathtt{rs2}}$$）。
$${\mathtt{y\_mix.mean(axis=0)}}$$ はこの右辺の 列方向の平均に相当します。

📈 （参考）Gauss–Hermite
近似精度を少し上げたい場合は $${\mathtt{rs2}}$$ の点数を増やす（例：300〜500）か、Gauss–Hermite を使います。

Gauss–Hermiteの参考コード

nodes, weights = np.polynomial.hermite.hermgauss(K)  # e^{-x^2} 用
rj = np.sqrt(2) * sigma * nodes
wj = weights / np.sqrt(np.pi)  # 正規(0,σ^2) に合わせる重み
pmf_mix = sum(w * stats.binom.pmf(y_val, N, expit(r)) for r, w in zip(rj, wj))

これは 無限混合の積分を「重み付き有限和」で近似する厳密派。
書籍のコードは「等確率分位点の等重み平均」という直感派。
どちらも正しい近似です。

📈 まとめ（このコードの本質）

1. 左の列：固定 r での二項分布 pmf（条件付きの世界）

2. 右の列：混ぜる時の重み（正規分布の pdf）

3. 最後の図：代表 r をたくさんとって 等重み平均 → 周辺 pmf（無限混合分布の近似）

これで、「条件付きの束を正規で混ぜると“太った分布”になる（過分散）」が 視覚と計算の両面で納得できるはずです。

うおーーーっ（なぜか叫びたくなるこの感じ…）
二項分布と正規分布を $${\int}$$（積分）するところが難しいですが、Pythonコード（テキストは R コード）で積分テクニックを駆使することで、無限混合分布を近似できることが分かりました。

ポアソン分布と正規分布の無限混合分布

ここから話題が変わります
テキストは続けて「ポアソン分布 $${\times}$$ 正規分布」の無限混合分布のケースに移ります。

テキストに準拠して「ポアソン分布 $${\times}$$ 正規分布」の無限混合分布を可視化します。
テキスト p.158 図 7.9 に相当します。

この図例は以下の条件で作成されます。

個体差 $${r=\{ -1.10,\ -0.30,\ 0.50,\ 1.30 \}}$$ をチョイスして、次のチャートを描きます。

# 「分布を混ぜる」という考え方Ⅱ ポアソン分布と正規分布の無限混合分布の例 p.158 図7.9

## 設定
N = 10                                         # 生存種子数の最大値
sigma = 1                                      # rのばらつきの大きさs
rs = [-1.10, -0.30, 0.50, 1.30]                # rのリスト
y_val_poi = list(range(N + 1))                 # ポアソン分布の横軸のyの値
r_val = np.linspace(-sigma, sigma, 100) * 2.7  # 正規分布の横軸のrの値
norm_dist = stats.norm(loc=0, scale=sigma)     # 平均0,標準偏差sigmaの正規分布の作成

## 描画
# 描画領域の設定
fig, ax = plt.subplots(4, 2, figsize=(6, 10), tight_layout=True)
# rの値ごとに二項分布（左）と正規分布（右）の描画を繰り返し処理
for i, r in enumerate(rs):
    ## 設定
    lam = np.exp(0.5 + r)       # ポアソン分布の平均の算出 λ=exp(beta1 + r)
    pdf_r = norm_dist.pdf(x=r)  # r=rのときの正規分布の確率密度関数の算出
    
    ## 個体差rごとに異なるポアソン分布の確率質量関数の描画
    # ポアソン分布の確率質量関数の描画
    ax[i, 0].plot(y_val_poi, stats.poisson.pmf(k=y_val_poi, mu=lam), '-o', ms=8, 
                  mec='white')
    # 修飾
    ax[i, 0].set(ylim=(-0.03, 0.62), xlabel='$y$',
                 title=f'$r=${r:.2f},  $\\lambda=${lam:.2f}')
    
    ## 集団内のrの分布：重みf(r|s)の正規分布の確率密度関数の描画
    # rの垂直点線の描画
    ax[i, 1].vlines(r, -0.03, pdf_r, color='black', ls='--', lw=1)
    # 正規分布の確率密度関数の描画
    ax[i, 1].plot(r_val, norm_dist.pdf(x=r_val))
    # f(r|s)の白抜き点の描画
    ax[i, 1].plot(r, pdf_r, 'o', ms=8, mec='tab:blue', mfc='white')
    # 修飾
    ax[i, 1].set(ylim=(0, 0.43), xlabel='$r$', title=f'$f(r)=${pdf_r:.2f}')

【実行結果】

続いて、すべての個体差 $${r}$$ に関する積分結果の可視化です。
テキストでは「集団をあらわす混合された分布」と紹介されています。

## （続き）描画２：集団全体をあらわす混合された分布

## 設定と準備
# 正規分布に従うrの200個のパーセンタイル点の取得
rs2 = norm_dist.ppf(q=np.linspace(0.001, 0.999, 200))
# 上記rで平均exp(0.5 + r)のポアソン分布に従う生存種子数yの確率の取得 shape=(200, 9)
y_mix_poi = np.array(
    [stats.poisson.pmf(k=y_val_poi, mu=np.exp(0.5 + r)) for r in rs2]
)

## 描画
# 混合された分布の描画
fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(y_val_poi, y_mix_poi.mean(axis=0), '-o', ms=8, mfc='white')
# 修飾
ax.set_xlabel('$y$', fontsize=12)
ax.set_ylabel('確率', fontsize=12)
ax.set(ylim=(0, 0.31));

【実行結果】

⏰️ アディショナルタイム ⏰️

（注意）
こちらは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

Gauss–Hermite 積分とモンテカルロ法

二項分布と正規分布の無限混合分布の可視化の際、ChatGPTが積分の方法「Gauss–Hermite 積分」を提言してくれました。
また記事外でChatGPTがそっと「モンテカルロ法」を教えてくれました。
この２つの方法で、無限混合分布に従う生存種子数 $${y_i}$$ の確率を可視化します！

◆ ◆ ◆

🔷 Gauss–Hermite 積分
ChatGPT の軽快な説明をご覧ください。

（１）全体像：３コマでスッキリ理解

① 条件付き分布の世界
個体差 $${r}$$ が与えられたら：

② 無限混合分布の世界（周辺化）
$${r \sim \mathcal N(0,s^2)}$$ を“混ぜる”ので

⇒ 二項×正規の無限混合分布（これが“周辺 pmf”）

③ Gauss–Hermiteで有限混合分布に近似
上の積分を、有限個の代表点 $${r_j}$$ と 重み $${w_j}$$ の和で置き換える：

⇒ 無限混合分布 ≈ 有限混合分布（スマートな点選び）

（２）数式：Gauss–Hermiteの“ちょい足し”だけ覚えればOK

$${\mathtt{numpy.polynomial.hermite.hermgauss(K)}}$$ の出力 $${\mathtt{nodes,\ weights}}$$ は、$${\int e^{-x^2} f(x)dx \approx \sum w_j f(x_j)}$$ 用になっている。
正規分布に合わせるには次のスケーリング処理が必要。

スケーリング処理の $${r_j,\ \widetilde w_j }$$ を用いて、有限混合分布の最終形を得る。

（３）コード例：Gauss–Hermite で周辺 pmf を算出

import numpy as np
from scipy.stats import binom
from scipy.special import expit

def binom_marginal_pmf_GH(y, n, beta1, beta2, x, sigma, K=20):
    # Gauss–Hermiteの実行 ※ノードと重みは「e^{-x^2} 用」
    nodes, weights = np.polynomial.hermite.hermgauss(K)
    # 正規分布(0, sigma^2) に合わせてノードrjと重みwjをスケーリング
    rj = np.sqrt(2.0) * sigma * nodes
    wj = weights / np.sqrt(np.pi)
    # 各ノードrjで二項分布pmfを計算して重み付き平均をとる
    qj = expit(beta1 + beta2 * x + rj)      # 各ノードの成功確率qj
    pmf_nodes = binom.pmf(y, n, pj)         # 二項分布pmf, shape=(K,)
    # 戻り値：周辺分布pmf p(Y=y)の近似値
    return np.sum(wj * pmf_nodes)

◆ ◆ ◆

🔷 モンテカルロ法
引き続き ChatGPT の軽快な説明をご覧ください。

やっていることはシンプルで「積分をサンプルの平均で近似する方法」です。

① 今回の周辺 pmf の場合、周辺 pmf は積分で書けます。

② モンテカルロ法 ではこのように近似します。

1. 個体差 $${r}$$ の $${R}$$ 個のサンプル $${r^{(1)}, r^{(2)}, \dots, r^{(R)} }$$ を $${\text{Normal}(0, \sigma^2)}$$ からランダムに生成

2. 各サンプルで $${P(Y=y \mid r^{(k)})}$$ を計算

3. その平均を取る：

③ ポイント

• Gauss–Hermite 積分 は “決まった位置と重み” で積分を近似

• モンテカルロ法 は “ランダムに引いた値” で積分を近似

• モンテカルロ法 は概念的に簡単ですが、精度を上げるにはサンプル数 $${R}$$ を多くする必要があります（遅くなりがち）

◆ ◆ ◆

🔷 実装
上記 ChatGPT の解説に寄り添ったコードを ChatGPT が書いてくれました。
先ほど描画したテキストの無限混合分布の「生存種子数 $${y}$$ ごとの確率」のチャートを Gauss–Hermite 積分とモンテカルロ法で再現します。

① 関数定義
２つの方法で無限混合分布の確率の近似値を計算する関数を定義します。

# インポート
import numpy as np
from scipy import stats
from scipy.special import expit, logsumexp

# ---------------------------------------------------------
# ユーティリティ関数：
# Gauss–Hermite 積分用のノード（積分点）と重みを、
# 正規分布 N(0, sigma^2) に合わせてスケーリングする
# ---------------------------------------------------------
def _gh_nodes_weights(sigma, K):
    # Gauss–Hermite 積分のノードと重み（元は e^{-x^2} 用）
    nodes, weights = np.polynomial.hermite.hermgauss(K)
    # ノードを N(0, σ^2) のスケールに変換
    rj = np.sqrt(2.0) * sigma * nodes
    # 重みも正規分布の形に合わせる
    wj = weights / np.sqrt(np.pi)
    return rj, wj  # rj: ノード値, wj: 重み（合計 ≈ 1）


# ---------------------------------------------------------
# 無限混合分布の周辺 PMF を Gauss–Hermite で近似
# 対応分布：二項分布（binomial） / ポアソン分布（poisson）
# ---------------------------------------------------------
def marginal_pmf_GH(dist, y, *, N=None, beta1=0.0, beta2=0.0, x=0.0,
                    sigma=1.0, K=80):
    """
    dist  : 'binomial' または 'poisson'
    y     : 計算したい観測値（スカラー or 1次元配列）
    N     : 二項分布の試行回数（dist='binomial' のとき必須）
    beta1 : 切片
    beta2 : 傾き（説明変数の係数）
    x     : 説明変数の値（スカラー）
    sigma : ランダム効果 r ~ N(0, sigma^2) の標準偏差
    K     : Gauss–Hermite の積分点数（大きいほど精度↑）

    戻り値: y と同じ形の周辺 pmf（numpy.ndarray）
    """
    # 分布タイプのチェック
    if dist not in {'binomial', 'poisson'}:
        raise ValueError("dist は 'binomial' か 'poisson' を指定してください。")
    if dist == 'binomial' and N is None:
        raise ValueError("dist='binomial' のときは N を指定してください。")

    # y を配列化（ベクトル化計算対応）
    y = np.atleast_1d(y)

    # Gauss–Hermite のノード rj と重み wj を取得
    rj, wj = _gh_nodes_weights(sigma, K)  # 形=(K,)

    # 各ノードでの線形予測子 η_j を計算
    eta_j = beta1 + beta2 * x + rj

    if dist == 'binomial':
        # 二項分布：η_j をロジスティック関数で成功確率 p_j に変換
        pj = expit(eta_j)  # 形=(K,)
        # 各 y に対する pmf(y | p_j) を計算（形=(len(y), K)）
        pmf_nodes = stats.binom.pmf(y[:, None], N, pj[None, :])
    else:
        # ポアソン分布：η_j を exp() で平均 λ_j に変換
        lam_j = np.exp(eta_j)
        pmf_nodes = stats.poisson.pmf(y[:, None], lam_j[None, :])

    # log-sum-exp で安定的に Σ wj * pmf(y | rj) を計算
    eps = 1e-300  # log(0) 回避用の極小値
    log_terms = np.log(wj)[None, :] + np.log(pmf_nodes + eps)
    pmf = np.exp(logsumexp(log_terms, axis=1))  # 形=(len(y),)

    # 戻り値：周辺 pmf ※入力がスカラーならスカラーで返す
    return pmf if y.ndim > 0 else pmf.item()


# ---------------------------------------------------------
# 無限混合分布の周辺 PMF を モンテカルロ法 で近似
# 対応分布：二項分布（binomial） / ポアソン分布（poisson）
# ---------------------------------------------------------
def marginal_pmf_MC(dist, y, *, N=None, beta1=0.0, beta2=0.0, x=0.0,
                    sigma=1.0, R=200_000, seed=0):
    """
    dist  : 'binomial' または 'poisson'
    y     : 計算したい観測値（スカラー or 1次元配列）
    N     : 二項分布の試行回数（dist='binomial' のとき必須）
    beta1 : 切片
    beta2 : 傾き
    x     : 説明変数の値（スカラー）
    sigma : ランダム効果 r ~ N(0, sigma^2) の標準偏差
    R     : モンテカルロサンプル数（大きいほど精度↑）
    seed  : 乱数シード

    戻り値: y と同じ形の周辺 pmf（numpy.ndarray）
    """
    # 分布タイプのチェック
    if dist not in {'binomial', 'poisson'}:
        raise ValueError("dist は 'binomial' か 'poisson' を指定してください。")
    if dist == 'binomial' and N is None:
        raise ValueError("dist='binomial' のときは N を指定してください。")

    # y を配列化
    y = np.atleast_1d(y)

    # 乱数生成器を初期化
    rng = np.random.default_rng(seed)

    # ランダム効果 r ~ N(0, σ^2) を R 個サンプリング
    r = rng.normal(0.0, sigma, size=R)

    # 線形予測子を計算
    eta = beta1 + beta2 * x + r

    if dist == 'binomial':
        # 二項分布の成功確率 p を計算
        p = expit(eta)  # 形=(R,)
        pmf_nodes = stats.binom.pmf(y[:, None], N, p[None, :])
    else:
        # ポアソン分布の平均 λ を計算
        lam = np.exp(eta)
        pmf_nodes = stats.poisson.pmf(y[:, None], lam[None, :])

    # 条件付き確率の平均（周辺 pmf のMC近似）
    pmf = pmf_nodes.mean(axis=1)

    # 戻り値：周辺 pmf ※入力がスカラーならスカラーで返す
    return pmf if y.ndim > 0 else pmf.item()

【実行結果】なし

② 二項分布＋正規分布の無限混合分布の近似値の算出と可視化
書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法の確率を比べます。

# 二項分布＋正規分布の無限混合分布の確率を算出

# Gauss–Hermite 積分
N = 8
ys = np.arange(N+1)
pmf_binom_gh = marginal_pmf_GH('binomial', ys, N=N,
                               beta1=-4, beta2=1, x=4, sigma=3, K=100)

# モンテカルロ法
pmf_binom_mc = marginal_pmf_MC('binomial', ys, N=N,
                               beta1=-4, beta2=1, x=4, sigma=3, R=500_000)

# 結果の表示（書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法）
y_mix_bin.mean(axis=0), pmf_binom_gh, pmf_binom_mc

【実行結果】
上から、書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法の順で、$${y=\{0, 1, 2, \cdots, 8\}}$$ のときの確率の近似値を表示しています。
３つの方法による確率はかなり近い感じがします。

続いて可視化です。

# 二項分布＋正規分布の無限混合分布の確率を可視化

# 書籍の計算方法による確率の描画
plt.plot(y_val_bin, y_mix_bin.mean(axis=0), '-o', color='gray',
         ms=8, mfc='white', label='書籍の計算方法')
# Gauss–Hermite 積分による確率の描画
plt.plot(pmf_binom_gh, color='tab:blue', label='Gauss-Hermite')
# モンテカルロ法による確率の描画
plt.plot(pmf_binom_mc, color='tab:red', ls='--', label='モンテカルロ近似')
# 修飾
plt.xlabel('生存種子数 $y$', fontsize=12)
plt.ylabel('確率', fontsize=12)
plt.ylim(0, 0.25)
plt.legend();

【実行結果】
３つの方法による確率の曲線は重なっています。
Gauss–Hermite 積分、モンテカルロ法でほどよく近似できたことが分かります。

③ ポアソン分布＋正規分布の無限混合分布の近似値の算出と可視化
書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法の確率を比べます。

# ポアソン分布＋正規分布の無限混合分布の確率を算出

# Gauss–Hermite 積分
ys_pois = np.arange(11)
pmf_pois_gh = marginal_pmf_GH('poisson', ys_pois,
                              beta1=0.5, beta2=0, x=0, sigma=1.0, K=80)
# モンテカルロ法
pmf_pois_mc = marginal_pmf_MC('poisson', ys_pois,
                               beta1=0.5, beta2=0, x=0, sigma=1.0, R=500000)

# 結果の表示（書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法）
y_mix_poi.mean(axis=0), pmf_pois_gh, pmf_pois_mc

【実行結果】
上から、書籍の計算方法、Gauss–Hermite 積分、モンテカルロ法の順で、$${y= \{0, 1, 2, \cdots, 10 \}}$$ のときの確率の近似値を表示しています。
３つの方法による確率はかなり近い感じがします。

続いて可視化です。

# ポアソン分布＋正規分布の無限混合分布の確率を可視化

# 書籍の計算方法による確率の描画
plt.plot(y_val_poi, y_mix_poi.mean(axis=0), '-o', color='gray',
         ms=8, mfc='white', label='書籍の計算方法')
# Gauss–Hermite 積分による確率の描画
plt.plot(pmf_pois_gh, color='tab:blue', label='Gauss-Hermite')
# モンテカルロ法による確率の描画
plt.plot(pmf_pois_mc, color='tab:red', ls='--', label='モンテカルロ近似')
# 修飾
plt.xlabel('生存種子数 $y$', fontsize=12)
plt.ylabel('確率', fontsize=12)
ax.set(ylim=(0, 0.31))
plt.legend();

【実行結果】
３つの方法による確率の曲線は重なっています。
Gauss–Hermite 積分、モンテカルロ法でほどよく近似できたことが分かります。

まとめ

一般化線形混合モデル GLMM および無限混合分布で「過分散のデータを解析できること」を学びました。

■ 過分散
仮定する確率分布で期待される分散よりもデータの分散が大きいことを 過分散 と呼びます。

■ 個体差・グループ差
観測されない（データ化できていない）個体差やグループ差が存在するが、過分散の原因になることがあります。

■ 一般化線形混合モデル
一般化線形混合モデルは、GLMに個体差・グループ差を組み込んだモデルです。
今回は個体差 $${r_i}$$ をランダム切片として追加します。
個体差は平均０、分散 $${s^2}$$ の正規分布に従うと仮定します。

■ 無限混合分布
パラメータを最尤推定できるように、二項分布と正規分布を「混ぜて」個体差 $${r_i}$$ を消します。
二項分布と正規分布を混ぜた分布を無限混合分布と呼びます。

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

次回は一般化線形混合モデルのパラメータ推定を実践します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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