Pythonで再現する標準ベイズ統計学10章

はじめに

本記事では標準ベイズ統計学の10章で掲載されている図表やモデルをPythonで実装する方法に関して説明します。

10章：非共役事前分布とメトロポリスヘイスティングスアルゴリズム

標準ベイズ統計学の10章では、この課題に対する強力なツールであるメトロポリス・ヘイスティングス（Metropolis-Hastings）アルゴリズムを紹介します。共役事前分布が使えない場合、従来のギブスサンプラーでは事後分布の推定が困難になります。そこで登場するのが、より汎用性の高いメトロポリス・ヘイスティングスアルゴリズムです。このブログでは、このアルゴリズムの基本的な考え方を特にポアソン回帰と経時回帰モデルという二つの実践的な応用例を通じて解説します。

年齢に対する子孫の数

図10.1では、52羽のメスのウタスズメを夏季に調査し、各個体の年齢と新しい子孫の数を記録した結果を、年齢に対する子孫の数の箱ひげ図として表現しています。

import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf

# データを読み込む
sparrows = pd.read_csv('./csv/sparrows.csv')
sparrows["age2"] = sparrows["age"] ** 2

# 必要なデータ列を取得
fledged_corrected = sparrows['fledged']
age_corrected = sparrows['age']

# 年齢を因子として扱い、年齢ごとにデータを再グループ化
sparrows['age_group_corrected'] = pd.Categorical(age_corrected)

# 箱ひげ図を再作成
plt.figure(figsize=(10, 6))
plt.boxplot([fledged_corrected[sparrows['age_group_corrected'] == ag] for ag in sorted(sparrows['age_group_corrected'].unique())], labels=sorted(sparrows['age_group_corrected'].unique()))
plt.xlabel('Age')
plt.ylabel('Fledged')
plt.title('Corrected Boxplot of Fledged by Age')
plt.show()

この図から、2歳のウタスズメが最も高い繁殖成功率の中央値を示し、2歳を超えると子孫の数が減少する傾向が見られます。この傾向は生物学的に説明可能です。1歳の鳥は最初の繁殖期であるため経験が浅く、一方で2歳を超えると健康状態や活動性が全体的に低下するためと考えられます。

格子点による$p(\beta_2|X,y), p(\beta_3|X,y), p(\beta_2, \beta_3|X,y)$の近似

図10.2は、ポアソン回帰モデルのパラメータ$β_2$と$β_3$の事後分布を格子点による近似で表現しています。事前分布$β \sim multivariate normal(0, 100 × I)$を用い、100個の格子点で近似した結果が示されており、$β_2$と$β_3$それぞれの周辺分布と両者の同時分布が描かれています。同時分布のグラフからは、$β_2$と$β_3$の間に強い負の事後相関が確認でき、これは確率が対角線に沿って集中していることを意味します。この結果、格子点による近似では多くの点の密度が0に近くなり、非効率的であることが示唆されています。この非効率性は、より複雑なモデルに対して計算上の困難をもたらす可能性があり、代替手法としてモンテカルロサンプルに基づく近似がより効率的であることが示唆されています。

import numpy as np
import scipy.stats as stats

# 初期パラメータを設定
p = 3
beta0 = np.zeros(p)
S0 = np.diag(np.repeat(100, 3))
gs = 100

# ベータ値のシーケンスを生成
beta1_seq = np.linspace(.27 - 2.5, .27 + 2.5, gs)
beta2_seq = np.linspace(.68 - 2, .68 + 2, gs)
beta3_seq = np.linspace(-.13 - .5, -.13 + .5, gs)

# 後方確率を計算
LPB = np.zeros((gs, gs, gs))
for i in range(gs):
    for j in range(gs):
        for k in range(gs):
            theta = beta1_seq[i] + beta2_seq[j] * age_corrected + beta3_seq[k] * age_corrected**2
            LPB[i, j, k] = stats.norm.logpdf(beta1_seq[i], beta0[0], np.sqrt(S0[0, 0])) + \
                           stats.norm.logpdf(beta2_seq[j], beta0[1], np.sqrt(S0[1, 1])) + \
                           stats.norm.logpdf(beta3_seq[k], beta0[2], np.sqrt(S0[2, 2])) + \
                           np.sum(stats.poisson.logpmf(fledged_corrected, np.exp(theta)))

# 正規化
PB = np.exp(LPB - np.max(LPB))
PB /= np.sum(PB)

#後方確率を計算
PB1 = np.sum(PB, axis=(1, 2))
PB2 = np.sum(PB, axis=(0, 2))
PB3 = np.sum(PB, axis=(0, 1))
PB23 = np.sum(PB, axis=0)

# サンプル数の設定
S = 50000
BETAg = np.zeros((S, 3))

for s in range(S):
    # PB2, PB23, PB の確率を正規化
    p_PB2 = PB2 / PB2.sum()
    i = np.random.choice(gs, p=p_PB2)

    p_PB23 = PB23[i, :] / PB23[i, :].sum()
    j = np.random.choice(gs, p=p_PB23)

    p_PB = PB[:, i, j] / PB[:, i, j].sum()
    k = np.random.choice(gs, p=p_PB)

    BETAg[s, :] = [beta1_seq[k], beta2_seq[i], beta3_seq[j]]


# BETAgの内容を確認
print(BETAg)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

# ベータ値に対する後方確率密度のプロット
plt.figure(figsize=(15, 3.5))

# beta2 のプロット
plt.subplot(1, 3, 1)
plt.plot(beta2_seq, PB2 * len(beta2_seq) / (max(beta2_seq) - min(beta2_seq)), 'b-')
plt.xlabel(r'$\beta_2$')
plt.ylabel(r'$p(\beta_2 | y)$')

# beta3 のプロット
plt.subplot(1, 3, 2)
plt.plot(beta3_seq, PB3 * len(beta3_seq) / (max(beta3_seq) - min(beta3_seq)), 'b-')
plt.xlabel(r'$\beta_3$')
plt.ylabel(r'$p(\beta_3 | y)$')

# 事後予測
Xs = np.column_stack((np.ones(6), np.arange(1, 7), np.arange(1, 7)**2))
eXB_post = np.exp(Xs @ BETAg.T)
qE = np.quantile(eXB_post, [0.025, 0.5, 0.975], axis=1)

# 2次元の最高密度領域(HDR)のプロット
plt.subplot(1, 3, 3)
data = BETAg[:, 1:3]
kde = gaussian_kde(data.T)
x, y = np.mgrid[min(data[:,0]):max(data[:,0]):100j, min(data[:,1]):max(data[:,1]):100j]
positions = np.vstack([x.ravel(), y.ravel()])
z = np.reshape(kde(positions).T, x.shape)
plt.contour(x, y, z)
plt.xlabel(r'$\beta_2$')
plt.ylabel(r'$\beta_3$')

plt.tight_layout()
plt.show()

正規モデルに対するメトロポリスアルゴリズムの結果

メトロポリスアルゴリズムは、事後分布を解析的に得ることが困難な状況で使用される手法です。このアルゴリズムでは、現在の値$θ(s)$から提案値θを生成し、受容率rに基づいてその提案を受け入れるか拒否するかを決定します。受容率rは事後確率の比$p(θ|y)/p(θ(s)|y)$で計算されますが、数値的安定性のために対数スケールで計算されることが多いです。
具体例として、分散既知の正規分布モデルを用いてメトロポリスアルゴリズムを適用しています。$σ² = 1, τ² = 10, μ = 5, n = 5$、および$y = (9.37, 10.18, 9.16, 11.60, 10.33)$というデータに対して、$θ(0) = 0$を初期値とし、提案分布$θ(s+1) \sim normal(θ(s), δ²)$, $δ² = 2$を用いて10,000回のイテレーションを実行しています。
図10.3は、このアルゴリズムの結果を視覚化しています。左図は生成された10,000個のθの値を繰り返し回数の関数として表示しており、初期値から始まってすぐに事後平均10.03付近に収束する様子が観察できます。右図は生成されたθの値のヒストグラムと、理論的な事後分布$normal(10.03, 0.20)$の密度関数を重ねて表示しています。

import scipy.stats as stats

# 初期パラメータ設定
s2 = 1
t2 = 10
mu = 5
np.random.seed(1)
n = 5
y = np.round(np.random.normal(10, 1, n), 2)

mu_n = (np.mean(y) * n / s2 + mu / t2) / (n / s2 + 1 / t2)
t2_n = 1 / (n / s2 + 1 / t2)

# MHアルゴリズム
theta = 0
delta = 2
S = 10000
THETA = []

for s in range(S):
    theta_star = np.random.normal(theta, np.sqrt(delta))
    log_r = (np.sum(stats.norm.logpdf(y, theta_star, np.sqrt(s2))) +
             stats.norm.logpdf(theta_star, mu, np.sqrt(t2))) - \
            (np.sum(stats.norm.logpdf(y, theta, np.sqrt(s2))) +
             stats.norm.logpdf(theta, mu, np.sqrt(t2)))

    if np.log(np.random.uniform()) < log_r:
        theta = theta_star

    THETA.append(theta)

# プロット
plt.figure(figsize=(14, 7))

# THETAの時系列プロット
plt.subplot(1, 2, 1)
plt.plot(range(10, S + 1, 10), THETA[9::10], 'b-')
plt.xlabel('Iteration')
plt.ylabel(r'$\theta$')

# THETAのヒストグラム
plt.subplot(1, 2, 2)
plt.hist(THETA[50:], bins=30, density=True, alpha=0.6)
th = np.linspace(min(THETA), max(THETA), 100)
theta_range = np.linspace(mu_n - 4*np.sqrt(t2_n), mu_n + 4*np.sqrt(t2_n), 100)
posterior_density = stats.norm.pdf(theta_range, mu_n, np.sqrt(t2_n))
plt.plot(theta_range, posterior_density, 'r-', lw=2, label='Theoretical Posterior')
plt.xlabel(r'$\theta$')
plt.ylabel('Density')

plt.tight_layout()
plt.show()

三つの異なる提案分布のもとでのマルコフ連鎖

図10.4は、メトロポリスアルゴリズムにおける提案分布の分散（$δ^2$）の影響を示しています。3つの異なる $δ^2$ の値（1/32, 2, 64）に対するマルコフ連鎖の振る舞いを比較しています。

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from statsmodels.tsa.stattools import acf

# 初期パラメータ設定
s2 = 1
t2 = 10
mu = 5
np.random.seed(1)
y = np.array([9.37, 10.18, 9.16, 11.60, 10.33])

# プロット用の設定
fig, axs = plt.subplots(1, 3, figsize=(15, 5))  # 横に3つ並べる
ACR = []
ACF = []
THETAA = []

# δ^2 の値を1/32, 2, 64に絞る
delta2_values = 2.0**np.array([-5, 1, 6])

# ループを修正して、指定したdelta2のみ計算とプロットを行う
for idx, delta2 in enumerate(delta2_values):
    THETA = []
    acs = 0
    S = 10000
    theta = 0

    for s in range(S):
        theta_star = np.random.normal(theta, np.sqrt(delta2))
        log_r = (np.sum(stats.norm.logpdf(y, theta_star, np.sqrt(s2))) -
                 np.sum(stats.norm.logpdf(y, theta, np.sqrt(s2)))) + \
                (stats.norm.logpdf(theta_star, mu, np.sqrt(t2)) -
                 stats.norm.logpdf(theta, mu, np.sqrt(t2)))

        if np.log(np.random.uniform()) < log_r:
            theta = theta_star
            acs += 1

        THETA.append(theta)

    # THETAの最初の1000点をプロット
    axs[idx].plot(THETA[:1000])
    axs[idx].set_title(f"Delta2 = {delta2}")

    # 受容率と自己相関を計算
    ACR.append(acs / S)
    ACF.append(acf(THETA, nlags=1)[1])

    # THETAAに追加
    THETAA.append(THETA)

plt.tight_layout()
plt.show()

この図は、適切な提案分布の選択がメトロポリスアルゴリズムの効率性に大きく影響することを示しています。$δ^2$ = 2 の場合が最も効率的なサンプリングを実現しており、パラメータ空間を十分に探索しつつ、適度な頻度で新しい値を受け入れています。

ポアソン回帰モデルにおける$\beta_3$のマルコフ連鎖と自己相関分析

スズメの繁殖データを用いてポアソン回帰モデルのベイズ推定を行っています。具体的には、スズメの年齢（$x$）と子孫の数（$Y$）の関係を$log E[Y|x] = β_1 + β_2x + β_3x^2$というモデルで分析しています。メトロポリス・ヘイスティングス（MH）アルゴリズムを使用して、モデルのパラメータ（$β_1$, $β_2$, $β_3$）をベイズ推定しています。結果を2つの主要な図でプロットしています。図10.5では$β_3$の時系列プロットと自己相関関数（ACF）プロットです。図10.6は$β_2$と$β_3$の周辺事後分布の密度プロット、および年齢に対する子孫の数の予測分布です。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import norm, poisson, uniform
from numpy.linalg import inv
from scipy.linalg import cholesky

# データ
fledged = sparrows['fledged']
age = sparrows['age']
age2 = sparrows['age2']

# MCMCアルゴリズムのセットアップ
n = len(fledged)
p = 3  # パラメータの数（intercept, age, age^2）
X = np.column_stack((np.ones(n), age, age2))
pmn_beta = np.zeros(p)
psd_beta = np.repeat(10, p)
var_prop = np.var(np.log(fledged + 1/2)) * inv(X.T @ X)
# beta = np.zeros(p)
beta = [0.2766, 0.6817, -0.1345]
S = 10000
BETA = np.zeros((S, p))
ac = 0
np.random.seed(100)

# 多変量正規分布からのサンプリング関数
def rmvnorm(n, mu, Sigma):
    return np.random.multivariate_normal(mean=mu, cov=Sigma, size=n)[0]

# MCMCアルゴリズムの実行
for s in range(S):
    # 新しいbetaを提案
    beta_p = rmvnorm(n=n, mu=beta, Sigma=var_prop)

    # 対数尤度比を計算
    lhr = (np.sum(poisson.logpmf(fledged, np.exp(X @ beta_p))) -
           np.sum(poisson.logpmf(fledged, np.exp(X @ beta))) +
           np.sum(norm.logpdf(beta_p, pmn_beta, psd_beta)) -
           np.sum(norm.logpdf(beta, pmn_beta, psd_beta)))

    # 受容/拒否の判定
    if np.log(uniform.rvs()) < lhr:
        beta = beta_p
        ac += 1

    BETA[s, :] = beta

print(ac / S)

# プロット設定
plt.figure(figsize=(15, 3.5))

# 時系列プロット
plt.subplot(1, 3, 1)
j = 2  # beta[3]に対応するインデックス
thin = np.arange(0, S, int(S/1000))  # 間引き
plt.plot(thin, BETA[thin, j], 'b-')
plt.axhline(y=np.mean(BETA[:, j]), color='red', linestyle='-')
plt.xlabel('Iteration')
plt.ylabel(r'$\beta_3$')

# 自己相関関数（ACF）プロット
plt.subplot(1, 3, 2)
plot_acf(BETA[:, j], ax=plt.gca(), color='gray')
plt.xlabel('Lag')
plt.ylabel('ACF')

plt.subplot(1, 3, 3)
plot_acf(BETA[thin, j], ax=plt.gca(), color='gray')
plt.xlabel('Lag/10')
plt.ylabel('ACF')

plt.tight_layout()
plt.show()

ポアソン回帰モデルのパラメータ推定とウタスズメの年齢別繁殖成功率の予測結果

# プロット設定
plt.figure(figsize=(15, 6.5))

# beta2の密度プロット
plt.subplot(1, 3, 1)
density_beta2 = gaussian_kde(BETA[:, 1])
x_beta2 = np.linspace(min(BETA[:, 1]), max(BETA[:, 1]), 1000)
plt.plot(x_beta2, density_beta2(x_beta2), 'b-', lw=2)
plt.xlabel(r'$\beta_2$')
plt.ylabel(r'$p(\beta_2 | y)$')

# beta3の密度プロット
plt.subplot(1, 3, 2)
density_beta3 = gaussian_kde(BETA[:, 2])
x_beta3 = np.linspace(min(BETA[:, 2]), max(BETA[:, 2]), 1000)
plt.plot(x_beta3, density_beta3(x_beta3), 'b-', lw=2)
plt.xlabel(r'$\beta_3$')
plt.ylabel(r'$p(\beta_3 | y)$')

# 予測値のプロット
plt.subplot(1, 3, 3)
Xs = np.column_stack((np.ones(6), np.arange(1, 7), np.arange(1, 7)**2))
eXB_post = np.exp(Xs @ BETA.T)
qE = np.quantile(eXB_post, [0.025, 0.5, 0.975], axis=1)
ages = np.arange(1, 7)
plt.plot(ages, qE[0, :], 'k-', lw=1)
plt.plot(ages, qE[1, :], 'k-', lw=2)
plt.plot(ages, qE[2, :], 'k-', lw=1)
plt.xlabel('Age')
plt.ylabel('Number of Offspring')

plt.tight_layout()
plt.show()

気温とCO2濃度のデータ

図10.7は、気温と二酸化炭素（CO2）濃度の関係を示しています。左の図は、気温とCO2濃度の時系列データを標準化して表示しています。右の図は、CO2濃度と現在の気温の差を散布図で表示しています。

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# データ
icecore = pd.read_csv('./csv/icecore.csv')

# データの標準化
icecore['standardized_CO2'] = (icecore['co2'] - icecore['co2'].mean()) / icecore['co2'].std()
icecore['standardized_temp'] = (icecore['tmp'] - icecore['tmp'].mean()) / icecore['tmp'].std()

# 設定
plt.figure(figsize=(15, 3))

# プロット1: 標準化された年別測定
plt.subplot(1, 3, 1)
plt.plot(icecore['year'], icecore['standardized_CO2'], label='CO2', color='gray')
plt.plot(icecore['year'], icecore['standardized_temp'], label='temp', color='black')
plt.xlabel('year')
plt.ylabel('standardized measurement')
plt.ylim(-2.5, 3)
plt.legend()

# プロット2: CO2 vs 温度
plt.subplot(1, 3, 2)
plt.scatter(icecore['co2'], icecore['tmp'])
plt.xlabel('CO2 (ppmv)')
plt.ylabel('temperature difference (deg C)')

plt.tight_layout()

最小二乗推定の残差分析

図10.8は、CO2濃度と気温の線形回帰モデルに関する残差診断を示しています。左の図は残差のヒストグラムを表示しており、この結果から正規性からの顕著な逸脱は見られないことがわかります。これは回帰モデルの誤差項が正規分布に従うという仮定を支持しています。右の図は残差の自己相関を示しています。

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
import statsmodels.api as sm
import pandas as pd
import numpy as np

# 線形回帰モデルのフィッティング
X = sm.add_constant(icecore['co2']) 
y = icecore['tmp']  
model = sm.OLS(y, X)
results = model.fit()

# 図とサブプロットの設定
plt.figure(figsize=(14, 7))  
plt.subplot(1, 2, 1)  # 1行2列のグリッドの最初のサブプロット
plt.hist(results.resid, color='blue', edgecolor='black', bins=20)
plt.xlabel('residual')
plt.ylabel('frequency')

plt.subplot(1, 2, 2)  # 1x2グリッドの2番目のサブプロット
plot_acf(results.resid, lags=40, alpha=0.05, use_vlines=True, ax=plt.gca())
plt.xlabel('lag')
plt.ylabel('Autocorrelation')

# レイアウトを調整
plt.tight_layout()  
plt.show() 

マルコフ連鎖によって生成された$\rho$の最初の1000個の値

図10.9の左図はメトロポリスヘイスティングスアルゴリズムによって生成された最初の1000個の値$\{\rho^{(1)} \dots \rho^{(1000)} \}$を図示しています。また、右の図は自己相関を示しています。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.graphics.tsaplots import plot_acf
from scipy.stats import norm, gamma
from scipy.linalg import solve

icecore = pd.read_csv('./csv/icecore.csv')

# MCMC関数
def rmvnorm(n, mu, Sigma):
    return np.random.multivariate_normal(mu, Sigma, n)

def mcmc(S, odens):
    n = len(icecore)
    y = icecore['tmp'].values
    X = np.column_stack((np.ones(n), icecore['co2'].values))
    DY = np.abs(np.subtract.outer(range(n), range(n)))

    lmfit = np.linalg.lstsq(X, y, rcond=None)[0]
    beta = lmfit
    s2 = np.var(y - X @ beta)
    phi = np.corrcoef(y[:-1], y[1:])[0, 1]
    nu0, s20 = 1, 1
    T0 = np.diag([1/1000] * 2)

    OUT = []
    ac = 0

    for s in range(S):
        Cor = phi ** DY
        iCor = np.linalg.inv(Cor)
        V_beta = np.linalg.inv(X.T @ iCor @ X / s2 + T0)
        E_beta = V_beta @ (X.T @ iCor @ y / s2)
        beta = rmvnorm(1, E_beta.flatten(), V_beta)[0]

        s2 = 1 / np.random.gamma((nu0 + n) / 2, 2 / (nu0 * s20 + (y - X @ beta).T @ iCor @ (y - X @ beta)))

        phi_p = abs(np.random.uniform(phi - 0.1, phi + 0.1))
        phi_p = min(phi_p, 2 - phi_p)
        lr = -0.5 * (np.linalg.slogdet(phi_p ** DY)[1] - np.linalg.slogdet(phi ** DY)[1] +
                     np.sum(np.diag((y - X @ beta).reshape(-1, 1) @ (y - X @ beta).reshape(1, -1) @
                                    (np.linalg.inv(phi_p ** DY) - np.linalg.inv(phi ** DY)))) / s2)

        if np.log(np.random.uniform()) < lr:
            phi = phi_p
            ac += 1

        if s % odens == 0:
            OUT.append(np.concatenate([beta, [s2, phi]]))

    return np.array(OUT), ac / S

# MCMC実行
np.random.seed(1)
OUT_1000, _ = mcmc(1000, 1)

# 図10.9
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].plot(OUT_1000[:, 3])
axes[0].set_xlabel('スキャン')
axes[0].set_ylabel('ρ')
plot_acf(OUT_1000[:, 3], lags=40, ax=axes[1])
axes[1].set_xlabel('lag')
plt.tight_layout()
plt.show()

25,000スキャンのマルコフ連鎖から25回のスキャンごとに抽出した$\rho$の値

図10.10は、マルコフ連鎖モンテカルロ法（MCMC）におけるシンニング（thinning）の効果を示しています。25,000回のスキャンのうち、25スキャンごとに結果を保存しています。これにより、元の100,000個のパラメータ値から1,000個のサンプルに減らしています。25,000回のスキャンのうち、25スキャンごとに結果を保存しています。これにより、元の100,000個のパラメータ値から1,000個のサンプルに減らしています。左側のプロットは、シンニングを適用した後のパラメータの値の推移を示しています。これにより、マルコフ連鎖の挙動を視覚化しています。右側のプロットは、シンニングを適用した後のサンプルの自己相関を示しています。シンニングを行うことで、サンプル間の自己相関が大幅に減少していることが確認できます。

OUT_25000, _ = mcmc(25000, 25)

# 図10.10
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].plot(OUT_25000[:, 3])
axes[0].set_xlabel('スキャン(25回ごと)')
axes[0].set_ylabel('ρ')
plot_acf(OUT_25000[:, 3], lags=40, ax=axes[1])
axes[1].set_xlabel('ラグ(25回ごと)')
plt.tight_layout()
plt.show()

$\beta_2$（傾きパラメータ）の事後分布と回帰直線の事後平均

図10.11の左図では$\beta_2$（傾きパラメータ）の事後密度に対するモンテカルロ近似を示しており、右図ではOLSとGLSの比較を行なっています。

# 図10.11
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# KDEを使用したプロット
kde = stats.gaussian_kde(OUT_25000[:, 1])
x_range = np.linspace(OUT_25000[:, 1].min(), OUT_25000[:, 1].max(), 200)
axes[0].plot(x_range, kde(x_range))
axes[0].fill_between(x_range, kde(x_range), alpha=0.5)
axes[0].set_xlabel('$β_2$')
axes[0].set_ylabel('事後周辺密度')

# 散布図
axes[1].scatter(icecore['co2'], icecore['tmp'])
axes[1].set_xlabel('$CO_2$')
axes[1].set_ylabel('気温')
b0, b1 = OUT_25000[:, :2].mean(axis=0)
axes[1].plot(icecore['co2'], b0 + b1 * icecore['co2'], 'k', linewidth=2, label='GLS推定')
lmfit = np.linalg.lstsq(np.column_stack((np.ones_like(icecore['co2']), icecore['co2'])), icecore['tmp'], rcond=None)[0]
axes[1].plot(icecore['co2'], lmfit[0] + lmfit[1] * icecore['co2'], 'gray', linewidth=2, label='OLS推定')
axes[1].legend()

plt.tight_layout()
plt.show()

最後に

本ブログでは、標準ベイズ統計学の第10章で扱われる非共役事前分布とメトロポリス・ヘイスティングスアルゴリズムについて、Pythonを用いた実装と可視化を行いました。次は11章の「線形・一般化線形混合効果モデル」について実装を行いたいと思います。