「データ解析のための統計モデリング入門」をPythonで写経 Vol.22 ～ 11章「空間構造のある階層ベイズモデル」①ベイズ統計モデリング～一次元空間の階層ベイズモデル

11章「空間構造のある階層ベイズモデル」

書籍の著者　久保拓弥 先生

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

この記事は 一次元空間の空間相関を考慮した階層ベイズモデル に取り組みます。
隣り合う場所が互いに相関し合うモデルに取り入れます。
場所差の事前分布を条件付き自己回帰の一種 Intrinsic CAR モデル（ICAR）で表現します。

アディショナルでは「空間相関の可視化」と「ICARのお気持ち」に取り組みます。

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

はじめに

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

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

準備

準備

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

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

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計計算
import scipy.stats as stats

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# Rデータセットの読み込み
import rdata

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

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

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

■ 統計モデル
一次元空間の空間相関を考慮するベイズ統計モデルです。
ランダム切片「場所差」の空間相関を $${\text{ICAR}}$$ で表現します。

■ モデリング手続き

データの確認

通常のデータの確認

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

◼️ データの読み込み
Y.RData ファイルを rdata ライブラリで変換して、pandas データフレームの data に読み込みます。

# 例題：一次元空間上の個体数分布のデータの読み込み p.242~

# RDataファイルの読み込み
data = rdata.read_rda('./data/ch11/Y.RData', default_encoding='ASCII')
data = pd.DataFrame(data)

# csvファイルの出力
# data.to_csv('./data/ch11/data.csv')

# データフレームの表示
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 50 です。
１本の直線上に等間隔に配置された「50 区画」を観察して記録した「個体数」に関する仮想の観測データです。
区画（あるいは位置）は インデックス $${j}$$ で示されます。

【変数の説明】
区画 $${j}$$ の個体数 y です。

50 の区画で構成される一次元空間のイメージ図です。

🍀🍀🍀

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

① 要約統計量の表示

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

【実行結果】
個体数 y の値は 0 ～ 19 です。

② 標本分散の表示

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

【実行結果】

y の分散は 27.373 です。
y の分布にポアソン分布を仮定すると標本分散は標本平均と近い値になるはずですが、分散は平均 10.880 の３倍近くの値であり、過分散が発生することになります。

③ y のヒストグラムの描画
個体数 y のヒストグラムで分布の形状を確認します。
データから推定したポアソン分布の確率質量関数もあわせて描画します。

# ヒストグラム, KDE, ポアソン分布の確率質量の描画

# yのヒストグラム・KDE曲線の描画
sns.histplot(x=data.Y, bins=8, stat='density', kde=True, ec='white',
             label='観測値');

# ポアソン分布の平均パラメータを推定して確率質量関数を描画
# scipy.stats の fit() で平均パラメータを推定
x_val = np.arange(0, data.Y.max()+1)
params = stats.fit(stats.poisson, data.Y, bounds=[(0, 100)]).params # λの推定
plt.plot(x_val, stats.poisson.pmf(x_val, *params), color='tab:red',
         label=f'$\\lambda=${params[0]:.1f}のポアソン分布')

# 修飾
plt.legend();

【実行結果】
y のヒストグラムとKDE曲線は平たくて分散が大きい感じがします。
一方でデータから推定した $${\lambda \approx 10.9}$$ のポアソン分布の確率質量関数は峰が高くて裾が重い（両端の確率が小さい）形状になっています。
過分散の現れなのかもしれません。

④ 区画と個体数の関係の可視化
テキスト p.244 図 11.2 に相当する区画と個体数 y・局所密度 m の関係を可視化します。

# 例題の一次元空間上の架空データの可視化 p.244 図11.2

# 区画 j と個体数 y の散布図の描画
plt.plot(data.Y, 'o', ms=7, alpha=0.7, label='観測値')
# （通常は分からない）区画 j の局所密度 m の赤点線の描画
plt.plot(data.m, color='tab:red', ls='--', label='平均値')
# 修飾
plt.xlabel('位置 $j$', fontsize=12)
plt.ylabel('個体数 $y_j$', fontsize=12)
plt.legend();

【実行結果】

位置は区画のことです。
個体数は各区間でばらついていますが、近くの区画、隣り合う区画では値が似ていることを確認できます。
局所密度（赤点線）は滑らかな曲線であり、やはり、近くの区画、隣り合う区画では値が似ていることを確認できます。

🍀🍀🍀

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

特徴①②より「ポアソン分布」「対数リンク関数」とし、特徴③より「場所差（ランダム切片）」を取り入れて、特徴④より「近傍区間の場所差の影響を考慮」する空間統計モデルをベイズ統計モデル化します。
今回取り組むモデルは 階層ベイズモデル です。

（注意）
アディショナルは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

アディショナル：空間相関の可視化

近傍区間に相関関係がある空間データを扱うということで…
ChatGPTに 空間相関の可視化手法を訊いてみました。
亜種含めて 10 種類教えてもらいました（うち２つは既知）。

優しさ度の高い３種の可視化で例題データの空間相関を満喫します。
すべて１次元空間データを想定してます。

描画用（およびモデリング用）に共通変数を作成します。

# 描画用・モデリング用データの作成

# 個体数 y
y = data.Y.values
# 標本サイズ
n = len(y)
# yの標本平均
y_mean = y.mean()
# y の標本標準偏差（分母はn）
y_std = y.std(ddof=0)
# 標準化したy
z = (y - y_mean) / y_std

【実行結果】なし

🍀🍀🍀

◼️ 類似度ヒートマップ（個体数の差）
２区画の「個体数の差（絶対値）」を２次元のヒートマップで可視化します。
シンプルならではの「直感的な理解」に期待します。

２つの区画の個体数の差が小さいセルは赤く、差が大きくなるにつれて薄赤⇒薄青⇒青へとセルの色が変わります。
近くの区画の個体数の差が小さい場合、左上から右下にかけての対角線付近が「赤く」なる予定です。

# 1-1) 類似度ヒートマップ：2地点の個体数の差のヒートマップ

# 2地点の個体数の差（絶対値）の算出
diff_abs = [[abs(y[i] - y[j]) for j in range(n)] for i in range(n)]

# 描画
plt.figure(figsize=(6, 5))
sns.heatmap(diff_abs, cmap='coolwarm_r', square=True,
            cbar_kws={'label':f'2地点の個体数の差（絶対値）'})
plt.title('類似度ヒートマップ：個体数の差')
plt.xlabel('位置 index j')
plt.ylabel('位置 index i')
plt.tight_layout()
plt.show()

【実行結果】

【チャートの見方】
行の区間 $${i}$$ と 列の区間 $${j}$$ の個体数の差（絶対値）が…
０：赤 ⇒ 薄い赤 ⇒ 10：白 ⇒ 薄い青 ⇒18：青 になっています。
近い区間は左上から右下への対角線に近いセルです。

【観察】
対角線付近のセル、つまり、近い区間同士は差が小さい＝似たような個体数ということが分かります。

さらには、対角線から遠いセルも赤が多い印象です。
これは個体数の分布が山なりになっているので、上り坂と下り坂で似たような値に遭遇するからでしょう。
例えば、区画 10 は 区画 20、30、40 付近にも同じような個体数になっているので、遠い区間とも相関しているように見えるのです。

【お知らせ】
PyMC モデリングで突如、「隣接行列 W」が出現します。
１つ隣の区間に１、その他の区間に０を設定した行列なのですが、ヒートマップで可視化すると、対角線上に「値あり」の色が付きます。
今回の類似度ヒートマップを用いて「近い区間は対角線上に現れる」ことに見慣れておくのが良さげです。

🍀🍀🍀

◼️ 自己相関・偏自己相関のコレログラム
時系列分析でおなじみの自己相関・偏自己相関です！
時系列分析は「時間という１次元の系列データ」を扱っているので、今回の「一次元空間の系列データ」にも転用できます（と思います）。

【一次元空間データのラグ】
区画の距離（離れ具合）をラグと呼びます。
横一列に等間隔で区画が並んでいるとすると…
　・左右の１つ隣の区画を「ラグ１」
　・左右の２つ隣までの区画を「ラグ２」
となります。

【自己相関と偏自己相関の違い】
自己相関コレログラムは、ラグごとの相関を可視化します。
隣はさらに隣の影響（自己相関）があり、隣の隣はさらにその隣の影響（自己相関）があるといった、影響のじわじわ伝達感を残したまま計算するのが自己相関です。
一方で、偏自己相関コレログラムは、影響のじわじわ伝達感を断ち切って、１つ隣だけの偏自己相関（ラグ１）、２つ隣だけの偏自己相関（ラグ２）を個別に計算します。

statsmodels の可視化関数 plot_acf（自己相関コレログラム）、plot_pacf（偏自己相関コレログラム）を利用します。

# 追加インポート
import statsmodels.tsa.api as tsa  # 自己相関・偏自己相関プロット

# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=True, tight_layout=True)
# 自己相関・偏自己相関のコレログラムの描画
tsa.graphics.plot_acf(y, lags=10, auto_ylims=True, title='自己相関', ax=ax[0])
tsa.graphics.plot_pacf(y, lags=10, auto_ylims=True, title='偏自己相関', ax=ax[1])
# 修飾
for ax_ in ax.flat:
    ax_.set_xlabel('ラグ', fontsize=12)
ax[0].set_ylabel('自己・偏自己相関係数', fontsize=12);

【実行結果】

【チャートの見方】
横軸はラグ、縦軸は係数（$${-1}$$ ～ $${1}$$）です。
青塗り領域は 95% 信頼区間です。

【観察】
自己相関コレログラムは、青塗り領域の外まで線が伸びているラグ２くらい、つまり、２つ隣まで自己相関がありそうなことを語っています。
偏自己相関コレログラムは、青塗り領域の外まで線が伸びているラグ１くらい、つまり、１つ隣まで偏自己相関がありそうなことを語っています。

ここでは偏自己相関を信じて「例題データの各区画の個体数は１つ隣の区画（両隣）の影響を受けているようだ」とまとめます。

🍀🍀🍀

◼️ モランの散布図
モランの散布図は初見です！
モランの $${I}$$ という統計量も関係するようです。

モランの散布図の概要を ChatGPT に訊きました。

ひとまず散布図を見てみましょう。
コードは長いです（ChatGPT製）。

# ===========================================================
# 1) （Moran's scatterplot）
#    数式との対応：
#      ・標準化値 z_i = (y_i - y_bar) / sigma_y
#      ・空間ラグ (Wz)_i = sum_j w_{ij} z_j
#      ・モランの I：
#         I = (n / S0) 
#             * [ sum_i sum_j w_{ij}(y_i - y_bar)(y_j - y_bar) ]
#             / [ sum_i (y_i - y_bar)^2 ]
#      ・行標準化（各行の重み和が1）を採用
# ===========================================================

def spatial_weights_1d_adjacent(n: int):
    '''
    一次元・等間隔・最近接（距離1）だけを「近い」とみなす隣接行列 W を作る。
    端点は隣が1つ、内部は隣が2つ。
    返り値：
        W : (n, n) の隣接行列（非標準化：隣接に1、それ以外0）
    '''
    W = np.zeros((n, n), dtype=float)
    for i in range(n):
        if i - 1 >= 0:
            W[i, i-1] = 1.0
        if i + 1 < n:
            W[i, i+1] = 1.0
    # 戻り値：隣接行列
    return W

def row_standardize(W: np.ndarray):
    '''
    行標準化：隣接行列の各行の和で割って、各行の和を1にそろえる。
    ※ブールインデックスを使わず、ブロードキャスト除算で安全に実装。
    '''
    # 行合計の算出 形状 (n, 1)
    row_sums = W.sum(axis=1, keepdims=True) 
    # 標準化したWを算出
    # ゼロ除算を避けるため、where を使って row_sums != 0 のところだけ割る
    W_std = np.divide(W, row_sums, out=np.zeros_like(W), where=row_sums != 0)
    # 戻り値：標準化した隣接行列
    return W_std

def morans_I(y: np.ndarray, W: np.ndarray):
    '''
    モランの I を定義式にしたがって計算する。
    ここでは y を「元の単位」のままで受け取り、教科書と同じ形で計算する。
      I = (n / S0) * [ Σ_i Σ_j w_ij (y_i - y_bar)(y_j - y_bar) ]
           / [ Σ_i (y_i - y_bar)^2 ]
    返り値：
        I : モランの I
    '''
    # 設定と準備
    y_bar = y.mean()    # 標本平均
    yc = y - y_bar      # yの中心化
    
    # モランのIの算出
    n = len(y)                                   # 標本サイズ
    S0 = W.sum()                                 # 重み総和
    num = (yc[:, None] * W * yc[None, :]).sum()  # 分子
    den = (yc ** 2).sum()                        # 分母 偏差平方和
    I = (n / S0) * (num / den)                   # モランのI
    
    # 戻り値：モランのI
    return I


def moran_scatterplot(z: np.ndarray, W_std: np.ndarray, title=None):
    '''
    モランの散布図を描く。
      横軸：z_i（標準化値）
      縦軸：(Wz)_i（近傍の平均＝空間ラグ）
    原点を通る回帰直線（最小二乗）の傾きは、行標準化 W のときモランの I に一致する。
    '''

    # 空間ラグの算出
    Wz = W_std @ z

    # 原点通過の回帰直線の傾きの算出（最小二乗）： slope = (z @ Wz) / (z @ z)
    slope = (z @ Wz) / (z @ z)

    # 描画領域の設定
    plt.figure(figsize=(6, 5))

    # 標準化したyと空間ラグの散布図の描画
    plt.scatter(z, Wz, s=40, alpha=0.7)
    
    # 回帰直線の描画
    xx = np.linspace(z.min(), z.max(), 200)
    plt.plot(xx, slope * xx, color='tab:red', lw=2)
    
    # x=0, y=0の点線の描画
    plt.axhline(0, color='gray', lw=1, ls='--')
    plt.axvline(0, color='gray', lw=1, ls='--')
    
    # グラフタイトルの表示
    title_text = 'モランの散布図（行標準化・隣接）'
    if title is not None:
        title_text += f'\n{title}'
    plt.title(f'{title_text}\n回帰直線の傾き（= モランの I）: {slope:.3f}')
    
    # 修飾
    plt.xlabel('$y_i$ の標準化値 $z_i$', fontsize=12)
    plt.ylabel('空間ラグ $(Wz)_i$：近傍の平均値', fontsize=12)
    plt.tight_layout()
    plt.show()


# ---- モラン散布図の実行例 ----

## 行標準化した隣接行列の算出
# 隣接行列の作成
W = spatial_weights_1d_adjacent(n)
# 隣接行列の（行）標準化
W_std = row_standardize(W)

# モランのIの算出 ※教科書の定義式どおりに計算
I_def = morans_I(y, W_std)

# モランの散布図の描画
# ※散布図の回帰直線の傾き（moran_scatterplot 内で計算）と I_def は一致する
moran_scatterplot(z, W_std, title=f'定義式でのモランの I ： {I_def:.3f}')

【実行結果】
直感的に「データには正の相関がありそう！」と思えます。

【チャートの見方】
横軸は自分の区画の個体数です。標準化しています。
縦軸は１つ隣の区画の個体数の平均値です。標準化した値を使っています。
ざっくり「自分の区画の個体数と、両隣の区画の個体数の平均値の散布図」なのです。

【観察】
自分の区画の個体数が大きくなるにつれて、両隣の区画の個体数も大きくなるのだなぁ、と読みました。
つまり、自分の区画の個体数と隣の区画の個体数は似ている、と読めます。

赤実線は、このデータ点の単回帰分析にて得られた回帰直線です。
回帰直線の傾きが「モランの $${I}$$」と呼ばれる空間自己相関の指標なのです。
ChatGPT によると指標の目安として、モランの $${I}$$ が 0.5 を超える場合、強い正の相関があるとのこと。
例題データは、１つ隣の空間自己相関は強い正の相関だ と言えそうです。

【読み解きヒント】
縦軸の「空間ラグ」は定義式などを調べると「？？？」ってなりますが…
安心して下さい！
今回のコードでは単に両隣の個体数の平均を計算しているだけです。
仕組みとしては、行標準化した隣接行列という難しい名前のデータを活用しています。
この行列を可視化しましょう。

# 行標準化した隣接行列 W_std の可視化・ヒートマップ

# 描画領域の設定
plt.figure(figsize=(6, 5))
# ヒートマップの描画
sns.heatmap(W_std, cmap='Greens', square=True)
# 修飾
plt.title('標準化した隣接行列')
plt.xlabel('位置 index')
plt.xlabel('位置 index');

【実行結果】

対角線上の緑色の値は 0.5 です。
行ごとに見ていくと、両隣だけ 0.5 で残りは 0 が設定されています。
この行に 50 区画の個体数を掛け算して結果を足すのです。
この計算は、両隣の個体数 $${\times 0.5}$$ を足し合わせているだけなのです。
つまり、両隣の個体数の平均値を計算しています。

ところで０行目と49行目の濃緑色の値は 1 です。
区画の両端は隣が１つしかありません。
反対側には区画がないからです。
ですので両端に限っては、片隣の１区画の個体数が両隣の個体数の平均値のフリをして、混入しています。

以上でアディショナル：空間相関の可視化は終了です。
面白かったですね！

🍀🍀🍀

◼️ （参考）ChatGPTと一緒に考えたチャート一覧
末尾★★★は記事で紹介しました。

• 相関ヒートマップ

  • 類似度ヒートマップ

    • 類似度：個体数の差の絶対値★★★

    • 類似度：標準化した個体数の外積

  • 局所類似度ヒートマップ

    • 類似度：小さな窓内の平均外積

  • ラグ相関ヒートマップ

    • 自己相関をスプライン補間でスムージング

    • ガウス過程カーネル相関

  • 局所相関ヒートマップ（小さな窓内の相関）

    • ピアソン相関

    • スピアマン相関

• 自己相関・偏自己相関のコレログラム★★★

• モランの散布図★★★

• セミバリオグラム

階層ベイズモデル（一次元空間・ICAR）

ベイズ統計モデルの理解

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

◼️ 今回モデルの概観
今回のベイズ統計モデルを「ICAR モデル」と呼びましょう。
ICAR は 条件付き自己回帰（conditional auto regressive, CAR）モデルの一種であり、intrinsic gaussian CAR モデルと分類されます。

推定するパラメータは $${\beta, \{r_j\}, s}$$ です（詳細は後ほど！）。
データを $${\bm Y}$$ とし、事後分布を数式化します。
$${n}$$ は標本サイズであり、例題データの場合は 50 です。

◼️ GLMM の構成要素
区画 $${j}$$ の個体数 $${y_j}$$ はポアソン分布 $${p(y_j \mid \lambda_j)}$$ に従うとします。
平均パラメータ $${\lambda_j}$$ は線形予測子と対数リンク関数を用いて $${\log \lambda_j = \beta + r_j}$$ とします。
$${\beta}$$ は全区画共通のパラメータ（大域的なパラメータ）です。
場所差 $${r_j}$$ は局所的なパラメータであり、標準偏差パラメータ $${s}$$ の ICAR モデル $${\text{ICAR}(s)}$$ で表現できるとします。

🍀🍀🍀

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

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

🍀🍀🍀

◼️ 事前分布
事後分布の数式化に現れた３つの事前分布を具体化します。

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

2️⃣ $${r_j}$$ の事前分布
標準偏差パラメータ $${s}$$ の ICAR を指定します。
階層事前分布に該当します。
ICAR の深掘りは後ほどアディショナルにて！

$${j \sim j'}$$ は区画 $${j}$$ と近傍するすべての区画 $${j'}$$ の組み合わせ $${\{j, j'\}}$$ です。
今回は１つ隣を近傍と扱うので、両隣がある区画の場合は $${\{j-1, j, j+1\}}$$、片隣しかない区画の場合は $${\{j, j+1\}}$$ か $${\{j-1, j\}}$$ です。

3️⃣ $${s}$$ の事前分布
無情報事前分布を設定します。
幅が十分に広い連続一様分布です。
テキストは上端を $${10^4}$$ としていますが、この記事のモデルでは $${10^3}$$ とします。

🍀🍀🍀

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

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

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

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

◼️ 隣接行列 W の作成
隣接行列 W は ICAR のパラメータです。
２つの区画が近傍かどうかを表す行列です。
今回のモデルは１つ隣を近傍とするので、区画 $${j}$$ の１つ隣の区画に１、その他の区画に０を設定します。

# 隣接行列Wの作成
W = np.array([[1 if abs(i - j)==1 else 0 for j in range(n)] for i in range(n)])
print('W.shape:', W.shape)
W

【実行結果】
行数・列数は例題データの標本サイズ 50 です。
２行目の区画 1 を見ましょう。
１が立っている区画は両隣に相当する 0 と 2 です。

隣接行列をヒートマップで可視化します。

# 隣接行列Wの可視化・ヒートマップ

# 描画領域の設定
plt.figure(figsize=(6, 5))
# ヒートマップの描画
sns.heatmap(W, cmap='Greens', square=True)
# 修飾
plt.title('隣接行列\n１つ隣の位置は１、その他の位置は０')
plt.xlabel('位置 index')
plt.ylabel('位置 index');

【実行結果】
対角線の両隣に濃緑色＝１が設定されています。
今回のモデルは、この隣接行列で「１つ隣を近傍として扱うこと」を指定しています。

テキスト p.246 に掲載の３点の「仮定」とこの隣接行列を関連付けてみます。

• 区画の場所差は近傍区画の場所差にしか影響されない

• 区画 $${j}$$ の近傍の個数 $${n_j}$$ は有限個であり、どの区画が近傍であるかはモデル設計者が指定する

• 近傍の直接の影響はどれも等しく $${1 / n_j}$$

隣接行列で指定した近傍の情報を用いて、近傍区画が決まります。
今回は１つ隣を近傍とするので、１つ隣の区画に 1 、その他の区画に 0 を（モデル設計者が）指定します。
近傍の個数 $${n_j}$$ は、両隣がある区画は 2 個、区画の両端は 1 個です。
近傍の直接の影響は、両隣がある区画は両隣の「場所差」から 1/2 ずつ影響を受け、区画の両端は隣の「場所差」だけから（ つまり 1 ）の影響を受けます。
この影響とは、テキスト p.246 の「$${r_j}$$ の条件付き事前分布に指定する正規分布の平均パラメータ $${\mu_j}$$」に相当します。

🍀🍀🍀

◼️ モデルの定義

# モデルの定義
# ※非中心パラメータ化：
# https://www.pymc.io/projects/docs/en/stable/api/distributions
# /generated/pymc.ICAR.html

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

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

    # 事前分布
    # 切片β: 無情報事前分布 N(0, 100^2)
    beta = pm.Normal('beta', mu=0, sigma=100)
    # 個体差rのICARモデルの標準偏差s: 無情報事前分布 U(0, 100) ※テキストはU(0, 1000)
    s = pm.Uniform('s', lower=0, upper=100)
    # 個体差r: 階層事前分布 ICAR(W)
    r_raw = pm.ICAR('r_raw', W=W, sigma=1, dims='data')  # ※ sigma=1
    r = pm.Deterministic('r', r_raw * s, dims='data')    # ※ r_raw * s

    # 線形予測子: 指数関数でlog(λ)を平均λに変換
    lam = pm.Deterministic('lam', pt.exp(beta + r), dims='data')

    # 尤度関数: 平均λのポアソン分布
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='data')

【実行結果】なし

【コードの補足：r を非中心パラメータ化】
MCMC サンプリングの効率向上などの目的で、階層ベイズモデルなどでは「非中心パラメータ化」と呼ばれるテクニックが使われます。

# 中心パラメータ化
# 標準偏差 = s で ICAR を設定
r = pm.ICAR('r', W=W, sigma=s, dims='data') 

# 非中心パラメータ化

# ICAR を 標準偏差=1 とし、場所差 r の素になる r_raw を設定
r_raw = pm.ICAR('r_raw', W=W, sigma=1, dims='data')

# 標準偏差パラメータ s と r_raw を掛け算して r を設定
r = pm.Deterministic('r', r_raw * s, dims='data')

🍀🍀🍀

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

# モデルの表示
model

【実行結果】

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

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

【実行結果】
個体差 $${\texttt{r}}$$ は区間 data 単位で推定されます。

🍀🍀🍀

◼️ MCMC の実行
MCMC サンプリングを行います。
NUTS サンプラーに nutpie を使います。
合計 4000 個の MCMC サンプルを得ます。

%%time
# MCMCサンプリング
# chain=4, draws=1000, tune=1000, thinなし, パラメータ初期値なし, NUTSサンプラー利用
# テキストのパラメータ: draws=10000, tune=100, chains=3

with model:
    idata = pm.sample(
        draws=1000, tune=1000, chains=4, random_seed=123,
        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', 'lam', 'r']
pm.plot_trace(idata, var_names=var_names, backend_kwargs={'tight_layout': True});

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

🍀🍀🍀

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

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

【実行結果】

テキストのパラメータ推定値と比較できるように、中央値、95% 信用区間を算出します。

# パラメータ beta, s の中央値と95%CIの算出
pd.DataFrame({
    'beta': np.quantile(az.extract(idata).beta.data, q=[0.025, 0.500, 0.975]),
    's': np.quantile(az.extract(idata).s.data, q=[0.025, 0.500, 0.975])
}, index=['CI 2.5%', 'median', 'CI 97.5%']).T.round(3)

【実行結果】
テキストのパラメータ推定値とほぼ同じになりました。

【テキストのパラメータ推定値】

🍀🍀🍀

◼️ 事後分布の確認
$${\beta, s}$$の事後分布を可視化します。

# 事後分布プロット β, s
pm.plot_posterior(
    idata, var_names=['beta', 's'], hdi_prob=0.95, round_to=3,
    ref_val=0, ref_val_color='tab:red',
    backend_kwargs={'tight_layout': True}, figsize=(10, 3));

【実行結果】

後々の MCMC サンプルを用いた分析に使えるよう、MCMC サンプルを numpy 配列化します。

# MCMCサンプルからデータの取り出し

def stack_idata(xa):
    return xa.stack(sample=('chain','draw')).values

betas = stack_idata(idata.posterior.beta)  # β
rs = stack_idata(idata.posterior.r)        # r
ss = stack_idata(idata.posterior.s)        # s

【実行結果】なし

切片 $${\beta}$$ を観測値スケールで確認します。

# =========================================
# 1) 切片 β （対数スケール、観測値スケール）
# =========================================

# betaが0以上の確率、betaの中央値・95%HDIの算出
pr_beta_pos = float((betas > 0).mean())
beta_med = np.quantile(betas, 0.5)
beta_lo, beta_hi = az.hdi(betas, hdi_prob=0.95)

# 観測値スケールでの中央値・95%HDIの算出
exp_beta = np.exp(betas)
exp_beta_med = np.quantile(exp_beta, 0.5)
exp_beta_lo, exp_beta_hi = az.hdi(exp_beta, hdi_prob=0.95)

# 結果の表示
print('【切片 β（対数スケール）】')
print(f'  Pr(β > 0) = {pr_beta_pos:.3f}')
print(f'  β: median={beta_med:.3f}, 95%HDI=[{beta_lo:.3f}, {beta_hi:.3f}]')
print('\n【観測値スケール】')
print(f'  rate ratio: median={exp_beta_med:.3f}, '
      f'95%HDI=[{exp_beta_lo:.3f}, {exp_beta_hi:.3f}]')

# （お好み）観測値スケールのヒストグラム
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(exp_beta, bins=40, alpha=0.85, edgecolor='white')
ax.set_xlabel('exp(β)', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
ax.set_title('$\\beta_1$ の事後分布（観測値スケール）')
plt.show()

【実行結果】
観測値スケールは各区画の個体数です。
中央値 9.7 個、95%HDI は 8.9 ～ 10.7 個です。
観測値スケールの事後分布はきれいなベル型です。

場所差のばらつき $${s}$$ を眺めます。

# =========================================
# 2) ばらつきの強さ：場所差の s
#  ・s は「log(λ) 尺度の標準偏差」
#  ・観測値スケール（λ）の“1SDシフト倍率”は exp(s)
# =========================================

# sの中央値、95%HDIを算出
s_med = np.quantile(ss, 0.5)
s_lo,  s_hi  = az.hdi(ss, hdi_prob=0.95)

# “1SD シフト時の倍率”の直感用（観測値スケール）：exp(sd)
s_factor_med  = float(np.exp(s_med))
print('【個体差 s（対数スケールのSD）】')
print(f'  s: median={s_med:.3f}, 95%HDI=[{s_lo:.3f}, {s_hi:.3f}]')
print(f'  1SDシフト倍率（λの直感）: exp(median s) ≈ ×{s_factor_med:.3f}')

# （お好み）s のMCMCサンプルのヒストグラム
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(ss,  bins=40, alpha=0.85, edgecolor='white')
ax.set_title('s の事後分布')
ax.set_xlabel('$s$', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
plt.tight_layout()
plt.show()

【実行結果】
1SDシフトの倍率は観測値スケールです。
ばらつきパラメータ $${s}$$ が 1シグマ 動くと、個体数は 1.26 倍になる、と読めます。

🍀🍀🍀

◼️ 事後予測チェック
y の事後予測値と観測値を可視化して比べます。

# 事後予測チェック

# 事後予測
with model:
    idata_pp = pm.sample_posterior_predictive(idata)

# 事後予測チェックプロット
pm.plot_ppc(idata_pp, num_pp_samples=100, random_seed=123);

【実行結果】
オレンジ点線は事後予測値の平均、薄青色線は MCMCサンプルのうち 100 個分の事後予測値、黒実線は観測値です。
事後予測値のオレンジ点線は右裾が長いものの、大まかに見て観測値の黒実線とまあまあ似ていると思っておきます。

種子数の予測

 $${\beta, r_j}$$ のMCMC サンプルを用いて、平均個体数 $${\lambda}$$ ベースの予測を可視化します。
テキスト p.250 図 11.4 の予測チャートに相当します。
※ $${\lambda}$$ の MCMC サンプルを用いても大丈夫です。

# 図11.2にモデルによる予測を追加 p.249~250 図11.4 ※テキストと若干異なる

## 設定と準備
# λの予測値の算出
lam_pred = np.exp(betas + rs)
# λの予測値の10%,50%,90%パーセンタイル点の算出
lam_pred_quantiles = np.quantile(lam_pred, q=[0.1, 0.5, 0.9], axis=1)

## 描画
# 観測値の散布図の描画
plt.plot(data.Y, 'o', ms=7, alpha=0.7, label='観測値')
# 観測値のベースになった平均λの折れ線グラフ（緑）の描画
plt.plot(data.m, color='tab:green', ls='--', label='観測値の平均λ')
# λの予測値（中央値）の折れ線グラフ（赤）の描画
plt.plot(lam_pred_quantiles[1, :], color='tab:red', label='λの予測値:中央値')
# λの予測値の80%区間の塗りつぶし描画
plt.fill_between(data.index, lam_pred_quantiles[0, :], lam_pred_quantiles[2, :],
                 color='lightpink', alpha=0.3, label='λの予測値:80%区間')
# 修飾
plt.xlabel('位置 $j$', fontsize=12)
plt.ylabel('個体数 $y_j$', fontsize=12)
plt.legend(bbox_to_anchor=(1, 1));

【実行結果】
青い点が y の観測値、緑点線が局所密度（観測値の平均 $${\lambda}$$）、赤実点が平均予測の中央値、薄赤色が平均予測の 80% 信用区間です。
テキストのチャートよりも 80% 信用区間の幅が狭い印象です。

【考察】
隣り合う区画の相関を捉えて、いい感じの予測になっていると思います。
このモデルが適しているように感じます！

相互作用する確率変数と確率場

テキストによると場所差 $${r_j}$$ は確率場に該当するそうです。

テキストの「ばらつきパラメータ $${s}$$ の大小が確率場 $${\{r_j\}}$$ に与える影響を図示」するためのシミュレーションを追体験します。

$${\beta=2.27}$$ に固定し、さらに $${\ s = \{0.0316, 0.224, 10.0\}}$$ の３つの固定パラメータを用いて、３つのモデルの MCMC サンプルを得て、$${\lambda_j}$$ の予測値を見ます。

３モデルの MCMC サンプリングを行います。

%%time
# ばらつきパラメータsの大小が確率場{r_i}に与える影響の可視化 １．MCMC p.251 図11.5

# 設定
beta_fixed = 2.27               # βの固定値
s_vals = [0.0316, 0.224, 10.0]  # sの3つの値
mcmc_samples = []               # 3つのMCMCサンプルを格納するリスト
coords = {'data': data.index}   # 座標ラベルの設定：データ行の識別子

# ベイズ統計モデリング：sの値ごとにモデリングとMCMCサンプリングを繰り返し処理
for s_val in s_vals:
    with pm.Model(coords=coords) as model2:
        
        ## モデリング
        # dataの定義: 目的変数=生存種子数Y
        Y = pm.Data('Y', value=data.Y.values, dims='data')

        # パラメータ（固定値）の設定 ※推定しない
        # 切片β
        beta = beta_fixed
        # 個体差rのICARモデルの標準偏差s
        s = s_val

        # 事前分布
        # 個体差r: 階層事前分布 ICAR(W) ※非中心:標準偏差sigma=1⇒r_raw * s
        r_raw = pm.ICAR('r_raw', W=W, sigma=1, dims='data')
        r = pm.Deterministic('r', r_raw * s, dims='data')
        
        # 線形予測子: 指数関数でlog(λ)を平均λに変換
        lam = pm.Deterministic('lam', pt.exp(beta + r), dims='data')

        # 尤度関数: 平均λのポアソン分布
        obs = pm.Poisson('obs', mu=lam, observed=Y, dims='data')


        ## MCMCサンプリング
        idata2 = pm.sample(
            draws=1000, tune=1000, chains=4, random_seed=123,
            nuts_sampler='nutpie',  # nutpieを使わない場合はこの行を削除
        )
    
    # MCMCサンプルをリストに格納
    mcmc_samples.append(idata2)

【実行結果】
３つの MCMC サンプルの Divergences は０個です。

🍀🍀🍀

３つのモデルの 平均個体数 $${\lambda_j}$$ の予測値を３本生成して、描画します。

# ばらつきパラメータsの大小が確率場{r_i}に与える影響の可視化 ２．描画 p.251 図11.5

## 設定
titles = ['A', 'B', 'C']
rng = np.random.default_rng(7)

## 描画
# 描画領域の設定
fig, axes = plt.subplots(3, 1, figsize=(6, 10) ,tight_layout=True)

# sの3つの値ごとにMCMCサンプルからrの取り出しと描画を繰り返し処理
for mcmc_sample, s_val, title, ax in zip(mcmc_samples, s_vals, titles, axes.flat):

    ## MCMCサンプルからrの値を取り出し
    r_samples = stack_idata(mcmc_sample.posterior.r_raw)

    ## 観測値の散布図の描画
    ax.plot(data.Y, 'o', ms=7, alpha=0.7)

    ## 3つのλの折れ線グラフの描画
    # rから3つのサンプルをランダムに抽出
    sample3 = rng.choice(a=len(r_samples), size=3)
    # βの固定値・抽出したrのサンプル・s_valを用いて、3つのλを算出
    lam_plot = np.exp(beta_fixed + r_samples[:, sample3] * s_val)
    # 3つのλの折れ線グラフの描画
    ax.plot(lam_plot, lw=3, alpha=0.5)

    # 修飾
    ax.set_ylabel('個体数 $y_j$', fontsize=12)
    ax.set(ylim=(-2, 28), title=f'({title}) $s={{{s_val}}}$')

# 修飾
axes[-1].set_xlabel('位置 $j$', fontsize=12);

【実行結果】
赤、黄色、緑の３本の線が生成した予測値です。

【観察】
ばらつきパラメータ $${s}$$ の値が大きいほど、予測値のばらつきが大きくなっています。

$${s = 0.0316}$$ のケースは、近くの区画とほぼ同じ値を予測しています。
テキストは「両隣の平均と似ている傾向が強くなり、$${r_j}$$ 全体のばらつきは小さくなります」と分析しています。

$${s = 0.24}$$ のケースは、観測値と近い予測しています。
ICAR モデルによる $${s}$$ の推定値とほぼ同じ値なので、観測値と近いのでしょう。

$${s = 10.0}$$ のケースは、近くの区画と値が大きく異なっています。
テキストは「$${s}$$ が大きいときには各 $${r_j}$$ は隣とは無関係に値を選べるようになり、各区画 $${k}$$ ごとにデータに合わせようとするので、ばらつきの大きい確率場になります」と分析しています。

大域的なパラメータ $${s}$$ が場所差 $${r_j}$$ たちの推定を「コントロール」しているようです。

（注意）
アディショナルは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

アディショナル：ICAR をちょっと深堀り

ChatGPT に ICAR の教科書 を作ってもらいました。
PyMC の ICAR クラスのページと接続できるように変数を
$${r_j \rightarrow \phi}$$、$${s \rightarrow \sigma}$$ に変えています。
また、「区画」を「地点」に読み替えてください。

長文をじっくりお読みいただき、ICAR のお気持ちに触れてみましょう。

1. ICAR は何をするための道具か

• 目的：「隣り合う場所は似た値になりやすい」 という空間の性質を、事前分布としてモデルに与えるための仕組みです。

• 対象：地図の地区、格子、一次元の位置列（今回の設定）など、どの地点がどの地点と隣かを定義できるデータ。

• 役割：平均構造にのる空間効果（地点 $${i}$$ ごとの“平滑なずれ”）を表すベクトル

に対して、「隣同士の差が大きくなりすぎる形」を嫌う“罰則”を与えます。

🚗🚗🚗

2. 隣接を行列で表す（W, D, Q）

2.1 隣接行列 $${W}$$

• $${W\in\mathbb{R}^{N\times N}}$$ は 隣接行列 です。

• $${W_{ij}=1}$$：地点 $${i}$$ と地点 $${j}$$ が 隣 を表す。

• $${W_{ij}=0}$$：隣ではない。

• 一次元空間（1,2,…,N）なら、$${W_{i,i-1}=W_{i,i+1}=1}$$、端点（1, N）は片側だけ 1 とする。

• 通常は対称行列になり、対角（自区画）には 0 を設定。

例：5 地点の一次元鎖の隣接行列

2.2 次数行列 $${D}$$

• $${D}$$ は対角行列。

• 対角成分 $${D_{ii}}$$ は「地点 $${i}$$ の隣接数」（一次元の端は 1、内部は 2）。

例：5 地点の一次元鎖の次数行列

2.3 ラプラシアン行列 $${Q=D-W}$$

• これが ICAR の“心臓部”。

• 一次元 5 点の例（鎖）：

• 重要な性質：任意の $${\boldsymbol{\phi}}$$ に対して

ここで $${\sum_{i\sim j}}$$ は「隣接するペア $${i,j}$$ についての和」。
つまり、隣同士の差の二乗の総和がこの二次形式で書けます。

🚗🚗🚗

3. ICAR 事前分布の形（確率密度）

ICAR 事前分布は次の“指数型”の形をしています：

• $${\sigma>0}$$ は スケール（標準偏差）。

• $${\sum_{i\sim j}(\phi_i-\phi_j)^2=\boldsymbol{\phi}^\top Q\boldsymbol{\phi}}$$ は罰則。

• 指数の前に マイナス が付いているので、差が大きい形ほど確率が下がる（起こりにくい）。

PyMCのICARクラスの「対数確率密度」の数式と比べます。

🚗🚗🚗

4. 「罰則が大きい／小さい」と「滑らかさ」のメカニズム

ここがとても重要な肝です。

4.1 罰則の正体
ICARの事前分布では、隣同士の差の二乗を合計した項

が入っています。

• これは「隣同士が離れているほど値が大きくなる」関数です。

• この項が指数関数の中でマイナス符号と一緒に登場するので、

  という形になります。

4.2 確率分布の観点から
上式の意味はこうです：

• 隣同士の差が大きいと $${(\phi_i - \phi_j)^2}$$ が大きくなる。

• すると指数の中が大きなマイナスになる。

• その結果、確率密度 $${p(\boldsymbol\phi)}$$ が小さくなる。

つまり「隣同士で値が離れすぎる形は確率的に起こりにくい」という制約がかかります。

4.3 σ² の役割
この「どのくらい強く嫌うか」を決めるのが $${\sigma^2}$$ です。

• $${\sigma^2}$$ が小さいと → 分母が小さい → 同じ差でも罰則が大きくなる → 隣はほぼ同じ値に押さえ込まれる。

• $${\sigma^2}$$ が大きいと → 分母が大きい → 差を許容しやすい → ギザギザ（起伏）が出やすい。

4.4 直感的な比喩
罰則は「隣同士をゴムひもで結んでいる」と考えるとわかりやすいです。

• 隣同士の差が大きい または $${\sigma^2}$$ が小さい = ゴムが強くて固い → 離れると強く引っ張られる → 隣同士がほぼ同じ高さになる → 曲線がなめらかに。

• 隣同士の差が小さい または $${\sigma^2}$$ が大きい = ゴムがゆるい → 引っ張りが弱い → 隣同士が自由に動ける → ギザギザの形も許される。

4.5 「滑らかさ」の意味
「滑らかになる」とは、隣の地点と値がほとんど変わらないので、グラフに描くと緩やかな曲線になることを指しています。
逆に「起伏を許す」とは、隣同士で急に変わることも認めるので、グラフがギザギザになります。

4.6 ✅ まとめると：

• 罰則が大きい（差が大 または σ²小） → 隣同士が強く結ばれ → 滑らか

• 罰則が小さい（差が小 または σ²大） → 隣同士がゆるく結ばれ → ギザギザもOK

🚗🚗🚗

5. ICAR が「内在的（intrinsic）」と呼ばれる理由（ゼロ和の話）

• $${Q=D-W}$$ は定数ベクトル $${\mathbf{1}=(1,\ldots,1)}$$ に対して $${Q\mathbf{1}=\mathbf{0}}$$。
  つまり $${\boldsymbol{\phi}^\top Q\boldsymbol{\phi}}$$ は、$${\boldsymbol{\phi}}$$ に定数を足しても変わらない。

• そのため、このままでは 正規化できない（不適切＝improper） 事前分布になります。

• 実務では、“ゼロ和” の弱い制約（$${\sum_i\phi_i\approx0}$$）を入れて、基準点を定めます。
  PyMC の $${\texttt{ICAR}}$$ でも、

という弱い正規化を対数密度に足して安定化しています（デフォルトはとても弱い制約）。

🚗🚗🚗

6. 一次元データでの具体像（今回の想定に直結）

• 隣接：隣り合う位置だけ（$${i\leftrightarrow i\pm1}$$）を 1 にした 0/1 の対称行列 $${W}$$。

• ラプラシアン $${Q=D-W}$$：内部は $${[,\ldots,-1,2,-1,\ldots]}$$、端は $${[1,-1]}$$ の三重対角。

• 罰則：$${\sum_{i=1}^{N-1}(\phi_i-\phi_{i+1})^2}$$。
  ⇒ 隣との段差を小さく保つように働く（=なめらかな線が好まれる）。

🚗🚗🚗

7. 観測モデルへの入れ方（概念図）

カウントデータ（種子数）なら、例えば：

• 観測値（ポアソン分布）：

• 線形予測子（対数リンク関数）：

ここで $${\alpha}$$ は切片、$${\mathbf{x}_i}$$ は説明変数、$${\beta}$$ は係数、$${\phi_i}$$ が ICAR の空間効果。

• 事前分布（ICAR）：

🚗🚗🚗

8. PyMC の `ICAR` クラスがしてくれること（実装の勘所）

• $${\texttt{W}}$$：隣接行列（0/1 が基本、重み付きも可）。対称・対角 0 を守る。

• $${\texttt{sigma}}$$：平滑化の強さ（標準偏差）。

  • 小さいほど強く平滑、大きいほど起伏を許す。

  • ふつうは $${\texttt{pm.HalfNormal}}$$ や $${\texttt{pm.HalfStudentT}}$$ などのハイパー事前分布を置いて学習させます。

• $${\texttt{zero\_sum\_stdev}}$$：ゼロ和制約の弱さ（実際に使われる標準偏差は $${N \times \texttt{zero\_sum\_stdev}}$$）。

  • デフォルトはとても弱い固定。特に理由がなければそのままでよい。

🚗🚗🚗

9. 「罰則の大きさ」と「滑らかさ」の関係をもう一度（腹落ち版）

• 罰則は $${\sum_{i\sim j}(\phi_i-\phi_j)^2}$$。

• 確率は $${\exp\bigl(-\frac{1}{2\sigma^2}\times\text{罰則}\bigr)}$$。

• 同じ差でも $${\sigma^2}$$ が小さければ 指数の減衰が速い ⇒ その形は起こりにくい ⇒ 隣は揃えられる。

• $${\sigma^2}$$ が大きければ 減衰が遅い ⇒ 差があっても受け入れる ⇒ ギザギザも許す。

• 比喩：硬いバネ（$${\sigma}$$ 小）＝すぐ戻る＝滑らか。ゆるいバネ（$${\sigma}$$ 大）＝自由＝起伏。

🚗🚗🚗

10. どこまでが ICAR？（近距離・局所のモデルだという自覚）

• ICAR は「近接差」へのペナルティで局所的に滑らかさを与えるモデルです。

• ガウス過程（カーネルで $${\rho(h)}$$ を指定）と違い、距離の連続関数は直接は使っていません。

• その分、隣接行列さえ定まれば軽量で扱いやすく、“隣は似る” という直感とダイレクトに対応します。

🚗🚗🚗

11. 仕上げの実務ポイント（軽くチェック）

1. 切片と $${\sum\phi_i=0}$$ の役割分担を守る（平均レベルは切片 $${\alpha}$$ が持つ）。

2. $${\texttt{sigma}}$$ の事前分布は弱すぎず強すぎず（HalfNormal など）にして、データから学習させる。

3. フィット後は 残差の空間相関（モランの $${I}$$ 等）を見て、取り残しがないかを確認。
   うまく入っていれば、残差の空間相関は小さくなります。

🚗🚗🚗

最後に（ひとことで ICAR）

まとめ

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

🔷 尤度
観測データはポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は対数です。
線形予測子は、切片、個体差（ランダム切片、空間相関あり）の和です。

🔷 事前分布
線形予測子の切片 $${\beta}$$ の事前分布は、平均０、標準誤差 100 の正規分布です。

個体差 $${r_j}$$ の事前分布は標準誤差 $${s}$$ の ICAR モデルであり、階層事前分布です。

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

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

次回は 空間相関モデルによる欠損値補完 を学びます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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