ホーム » コードあり » 2019 Q5(3)

2019 Q5(3)

投稿者: Qchan in コードあり, ベイズ統計, ラプラス分布, 復習3周目, 統計検定1級 2019年統計数理投稿日: 2024年3月11日最終更新日: 2025年2月4日

前問の事後確率密度関数の最大値を取る期待値μを求めました。

追記

■2025/2/4 $\max\left(, \xi\right)$ や $\min\left(, \xi\right)$ の導き方

$h(\mu) = -\frac{1}{2} n (\mu - \bar{y})^2 - \lambda |\mu - \xi|$

を変形し

$h(\mu) = -\frac{1}{2} n \{(\mu - \xi) - (\bar{y} - \xi)\}^2 - \lambda |\mu - \xi|$

となる。

ここで、 $(\mu - \xi)$ と $(\bar{y} - \xi)$ の符号を考えると、同符号のとき、第一項は異符号のときと比べて大きくなる。なぜなら

AとBが同符号のときの $-(A-B)^2$ のほうが

AとBが異符号のときの $-(A-B)^2$ よりも大きいからだ。

一方で第二項 $- \lambda |\mu - \xi|$ は変化しない。よって、

$\hat{\mu} - \xi \quad$ と、 $\bar{y} - \xi$ は同符号となるから

$\bar{y} > \xi$ のときは $\hat{\mu} > \xi$ と言えるので

$\hat{\mu}$ は、 $\xi$ を下回らないので

$\hat{\mu} = \max\left(\bar{y} - \frac{\lambda}{n}, \xi\right)$ のように書ける。

同様に

$\bar{y} < \xi$ のときは $\hat{\mu} < \xi$ と言えるので

$\hat{\mu}$ は、 $\xi$ を上回らないので

$\hat{\mu} = \min\left(\bar{y} + \frac{\lambda}{n}, \xi\right)$ のように書ける。

まとめると

$\hat{\mu} = \begin{cases} \max \left( \bar{y} - \frac{\lambda}{n}, \xi \right) & (\bar{y} > \xi) \\ \xi & (\bar{y} = \xi) \\ \min \left( \bar{y} + \frac{\lambda}{n}, \xi \right) & (\bar{y} < \xi) \end{cases}$

となる。

コード

対数事後分布 h(μ)と、その第1項（観測データに基づく項）および第2項（事前分布に基づく項）のグラフを描画します。パラメータ ξ を 1 または 5、スケールパラメータ λ を 25 または 50 として、それぞれの組み合わせで描画します。

# 2019 Q5(3)  2024.9.29

import numpy as np
import matplotlib.pyplot as plt

# 対数事後分布の第1項: 観測データに基づく項
def h_term1(mu, y):
    n = len(y)
    y_mean = np.mean(y)
    return -0.5 * n * (mu - y_mean)**2

# 対数事後分布の第2項: 事前分布に基づく項
def h_term2(mu, lambd, xi):
    return -lambd * np.abs(mu - xi)

# 対数事後分布 h(μ) の定義
def h(mu, y, lambd, xi):
    return h_term1(mu, y) + h_term2(mu, lambd, xi)

# μハット（事後分布のモード）を計算
def mu_hat(y_mean, lambd, xi, n):
    if y_mean > xi:
        return max(y_mean - lambd / n, xi)
    elif y_mean < xi:
        return min(y_mean + lambd / n, xi)
    else:
        return xi

# パラメータ設定のリスト (xi, lambd) の組み合わせ
params = [(1, 25), (1, 50), (5, 25), (5, 50)]
n_samples = 20  # サンプル数
mu_values = np.linspace(0, 6, 500)  # μの範囲

# グラフを4つ表示するための設定
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for i, (xi, lambd) in enumerate(params):
    # 観測データを生成して、平均を3に調整
    np.random.seed(43)
    y = np.random.normal(loc=3, scale=1, size=n_samples)
    y_mean = np.mean(y)
    
    # μハットを計算
    mu_hat_value = mu_hat(y_mean, lambd, xi, n_samples)

    # 各項の計算
    h_values = h(mu_values, y, lambd, xi)
    h_term1_values = h_term1(mu_values, y)
    h_term2_values = h_term2(mu_values, lambd, xi)

    # サブプロットに描画
    ax = axes[i//2, i%2]
    ax.plot(mu_values, h_values, label='h(μ)', color='orange')
    ax.plot(mu_values, h_term1_values, label='第1項: -n/2 (μ - ȳ)²', color='blue')
    ax.plot(mu_values, h_term2_values, label='第2項: -λ|μ - ξ|', color='green')
    
    # μハットを縦線でプロット
    ax.axvline(mu_hat_value, color='red', linestyle='--', label=f'μハット = {mu_hat_value:.2f}')

    # グラフのタイトル
    ax.set_title(f'ξ = {xi}, λ = {lambd}, ȳ = {y_mean:.2f}, n = {n_samples}')
    ax.set_xlabel('μ')
    ax.set_ylabel('値')
    ax.legend()
    ax.grid(True)

# グラフを描画
plt.tight_layout()
plt.show()

対数事後分布 h(μ)が最大値を取るμハットが

であることが確認できました。

カテゴリー

投稿一覧