「自然科学の統計学」を読む　～第9章ベイズ決定 ①事前確率分布・事後確率分布～

はじめに

東京大学教養学部統計学教室編「自然科学の統計学」は，1992年発行のやや古典的な文献であるが，自然科学に関わる統計学的テーマが簡潔にまとめられている。数理統計学の復習も兼ねて，本書を読むこととした。
ただ，基本的なことは他書で学んできたのと，本書自体がかなり細かく説明されているので，本書内の内容や数式を細かく追うというより，実務や統計検定の受験において有用そうなことを選んでまとめてみたい。

本記事は，「第9章ベイズ決定」における，事前確率分布と事後確率分布に関する読書メモである。

第9章ベイズ決定

本章では，ベイズの定理が持つ多くの有用な性質を活用するベイズ統計学を扱う。ベイズ統計では，人間が持つ予想，確信，信念などを主観確率として，判断や決定のために取り入れることができる。

9.2 事前確率分布と事後確率分布

ベイズ統計では，確率分布のパラメータ自体に確率分布を設定する。

十分に情報がないときにおけるパラメータ $\theta$ の確率分布を事前確率分布と呼ぶ。一般に事前確率分布は，

パラメータ $\theta$ の，実験前の予想をなるべく正確に反映するものである
事後確率分布の計算が容易になる

ように設定する。

例 : 新薬開発

開発中の新薬が，偽薬(プラセボ)と比較して有効となる確率を $\theta, \: (0 \lt \theta \lt 1)$ とする。
新薬がプラセボに対して有効であると判定されることを「成功」と表現すると， $n$ 回の試行中， $x$ 回の成功を得る確率は二項分布

$\begin{align} f(x \mid \theta) = {}_n C_x \theta^x (1 - \theta)^{n-x}, \quad x=0,1,..., n \\ \\ \end{align}$

にしたがう。なおこれをパラメータ $\theta$ の関数とみなすと，これは尤度であることが分かる。

$\theta$ の事前確率分布として，ベータ分布

$\begin{align} w(\theta) = Be(\theta \mid \alpha, \beta) = \frac{1}{B(\alpha, \beta)} \theta^{\alpha - 1} (1 - \theta)^{\beta - 1}, \quad (0 \lt \theta \lt 1) \\ \\ \end{align}$

を設定すると，事後確率分布はベイズの定理より，

$\begin{align} w(\theta \mid x) &\propto f(x \mid \theta) w(\theta) \\ \\ &\propto \theta^x (1 - \theta)^{n-x} \times \theta^{\alpha - 1} (1 - \theta)^{\beta - 1} \\ \\ &\propto \theta^{\alpha + x - 1} (1 - \theta)^{\beta + n - x - 1} \\ \\ &\propto Be(\theta \mid \alpha + x - 1, \beta + n - x - 1) \\ \\ \end{align}$

となり，事前確率分布と同様にベータ分布となる。

このように，尤度 $f$ に対して，事前確率分布と事後確率分布が同一種類の分布族になるとき，事前確率分布・事後確率分布の分布族を $f$ の自然な共役分布の族と呼ぶ。

数値計算例

本書P256を例に，事後確率分布の変化を可視化してみる。

データがないとき，事前確率分布には $Be(1, 1)$ を用いる。これは一様分布と同じ形であり，新薬が有効である確率が0から1まで値を取る確率がまんべんなく等しい状態を示している。

この状態から，成功 $S$ と失敗 $F$ に関する以下のデータが得られたとする。

$\begin{align} S, F, S, S, S \\ \\ \end{align}$

このデータを用いると，事後確率分布は，

$\begin{align} Be(2, 1) \rightarrow Be(2, 2) \rightarrow Be(3, 2) \rightarrow Be(4, 2) \rightarrow Be(5, 2) \\ \\ \end{align}$

のように変化する。この変化の様子を可視化したものを以下に示す。

描画用コードはクリックで展開

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

# Define the parameters for the Beta distributions
params = [
    (1, 1), (2, 1), (2, 2),
    (3, 2), (4, 2), (5, 2)
]

# Create a 2x3 grid of subplots
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
axes = axes.flatten() # Flatten the axes array for easy iteration

# Generate x values
x = np.linspace(0, 1, 100)

# Plot each Beta distribution
for i, (a, b) in enumerate(params):
    y = beta.pdf(x, a, b)
    axes[i].plot(x, y, label=f'Be({a}, {b})')
    axes[i].set_title(f'Be({a}, {b})')
    axes[i].set_xlabel('x')
    axes[i].set_ylabel('Probability Density')
    axes[i].legend()
    axes[i].grid(True)

    # Fill the area for Be(5, 2) where x >= 0.5
    if (a, b) == (5, 2):
        x_fill = np.linspace(0.5, 1, 100)
        y_fill = beta.pdf(x_fill, a, b)
        axes[i].fill_between(x_fill, y_fill, color='orange', alpha=0.5)


plt.tight_layout()
plt.show()

最終的に得られる事後確率分布は，あくまでパラメータ $\theta$ の確率分布であるので， $\theta$ に関する特定の値を評価したい場合，この確率分布の期待値や上側確率を計算する必要がある。
最終的に得られた事後確率分布 $Be(5, 2)$ について， $\theta \geq 0.5$ となる確率を求めると， $P(\theta \geq 0.5) = 57/64$ となり，かなり確からしいと言える。

まとめと感想

今回は「第9章ベイズ決定」における，事前確率分布と事後確率分布についてまとめた。

ベイズ統計においては，データを用いたパラメータ推定を，事後確率分布を用いて行なう。そのため事前確率分布の設計の指針は，事前の予想を正確に反映するものであること，また事後確率分布の計算が容易になること，の2点が挙げられていた。ただこれは基本的な方針であり，実業務における利用シーンでは尤度も事前確率分布もかなり複雑になり，事後確率分布に含まれる積分計算が実行できなくなるため，各種の数値計算手法が提案されている。

一方で，基本的な分布を用いて，事後確率が更新されていく様子を可視化して理解することは重要である。今回示した例では，初めは事前確率分布はフラットな形をしていたが，データが増えるにつれて確率密度関数の形状が偏り，特定の値の確率値が大きくなっていった。このように，基本的な例においてデータが増えることにより事後確率分布が変化していく様子を理解しておくと，事後確率分布の計算が複雑になってもイメージが付きやすいのではないかとかんがえられる。

本記事を最後まで読んでくださり，どうもありがとうございました。