ベイズ推定

ベイズ推定（Bayesian inference） は、未知パラメータを確率変数として扱い（伝統的統計学ではパラメータは定数）、データを観測する前の知識（事前分布）を、観測データによって更新する推定方法

確率分布再訪¶

確率分布¶

確率分布（確率密度関数） $p(x)$ は、次の形で表すことができる：

p(x) = \frac{1}{Z}\, f(x)

(1)

ここで

$f(x)$ ：カーネル（kernel）
$Z$ ：正規化定数（normalizing constant）

である。

例：正規分布

カーネルを

f(x) = \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)

(2)

正規化定数を

Z = \sqrt{2\pi\sigma^2}

(3)

とおくと、正規分布は次のように表すことができる

p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) =\frac{1}{Z}\, f(x)

(4)

カーネル¶

確率分布の カーネル（kernel） とは、正規化定数を除いた、確率分布の「形」だけを表す関数である。

記号的には

p(x) \propto f(x)

(5)

と書かれることが多く、この $\propto$ は「 $p(x)$ は $f(x)$ に比例する（定数倍を除いて同じ）」という意味を持つ。

正規化定数¶

カーネル $f(x)$ がどれほど「確率っぽい形」をしていても、確率の定義上、積分（和）が1でなければ確率分布ではない。そこで、 正規化定数（normalizing constant） $Z$ は、確率分布の条件 $\int p(x)\,dx = 1$ を満たすために導入される。

$Z$ は $Z = \int f(x)\,dx$ （離散の場合は $Z=\sum_x f(x)$ ）と定義される。

特に高次元では

Z = \int_{\mathbb{R}^d} f(x)\,dx

(6)

の計算が極めて困難、あるいは解析的に不可能になる。

そのため実務では

$Z$ を明示的に計算しない
$Z$ を含まない比だけを使う

という戦略が取られる。

直感的な理解¶

カーネルは「どこが高く、どこが低いか」を決める
正規化定数は「全体をどれだけ縮めるか」を決める

したがって、相対的な大小関係だけが重要な場面では、カーネルだけ分かっていれば十分である。

例えばMAP推定では

\arg\max_\theta p(\theta \mid y) = \arg\max_\theta f(\theta)

(7)

となるため、正規化定数は完全に不要である。

ベイズ推定での例¶

事後分布は

p(\theta \mid y) = \frac{p(y \mid \theta)p(\theta)}{p(y)}

(8)

と書けるが、通常は

p(\theta \mid y) \propto p(y \mid \theta)p(\theta)

(9)

と表現される。

ここで

カーネル： $p(y \mid \theta)p(\theta)$
正規化定数： $p(y)=\int p(y \mid \theta)p(\theta)\,d\theta$

であり、 $p(y)$ は 周辺尤度（evidence） と呼ばれる。

例：コイン投げ¶

表が出たら1、裏が出たら0として以下のようなデータがとれたとする。

[1,0,1,1,0,1,0,1,1,0]

事前分布にベータ分布、尤度をベルヌーイ分布とする

import numpy as np
import pymc as pm
import arviz as az

# 観測データ
y_obs = np.array([1,0,1,1,0,1,0,1,1,0], dtype=int)  # 0/1
N = y_obs.size

with pm.Model() as model:
    p = pm.Beta("p", alpha=1, beta=1)
    y = pm.Bernoulli("y", p=p, observed=y_obs)
    idata = pm.sample(2000, tune=1000, chains=4, progressbar=False, random_seed=42)

az.summary(idata, var_names=["p"])

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

Loading...

点推定¶

EAP (expected a posteriori, 事後期待値)¶

事後分布の期待値を推定値とする方法

\hat{\theta}_{\text {eap }} =E[\theta \mid \boldsymbol{x}] =\int \theta f(\theta \mid \boldsymbol{x}) d \theta=\int \theta \frac{f(\boldsymbol{x} \mid \theta) f(\theta)}{f(\boldsymbol{x})} d \theta

(10)

eap = idata.posterior.mean(dim=["chain", "draw"])
float(eap["p"])

0.5793009566171891

MAP (maximum a posteriori, 事後確率最大値）¶

事後分布の最大値（mode）を推定値とする

\hat{\theta}_{\text {map }}=\max _\theta f(\theta \mid \boldsymbol{x})

(11)

with model:
    map_est = pm.find_MAP()
float(map_est["p"])

Loading...

0.5999999988336109

MED (posterior median, 事後中央値）¶

事後分布の中央値、すなわち、累積分布関数が0.5になる点を推定値とするもの

F\left(\hat{\theta}_{m e d} \mid x\right)=\int^{\hat{\theta}_{\text {med }}} f(\theta \mid \boldsymbol{x}) d \theta=1 / 2

(12)

事後分散¶

EAP推定量の散らばりの目安には、事後分布の分散である 事後分散（posterior variance） と、その平方根をとった 事後標準偏差（posterior standard deviation） が用いられる。

V[\theta] =V[\theta \mid \boldsymbol{x}] =E\left[\left(\theta-\hat{\theta}_{e a p}\right)^2 \mid \boldsymbol{x}\right] =\int\left(\theta-\hat{\theta}_{e a p}\right)^2 f(\theta \mid \boldsymbol{x}) d \theta

(13)

std = idata.posterior.std(dim=["chain", "draw"])
float(std["p"])

0.13651722214261366

予測分布¶

パラメータの分布ではなく、データの予測値を出すもの。

条件付き予測分布¶

モデルの分布を $f(x\mid \theta)$ とした場合、 条件付き予測分布（conditional predictive distribution） は $f(x^* \mid \hat\theta)$ である。 $\hat\theta$ はEAP推定量など事後分布からの推定量を用いる。

事後予測分布¶

事後予測分布（posterior predictive distribution）

f(x^* \mid \boldsymbol{x})=\int_{-\infty}^{+\infty} f(x^* \mid \theta) f(\theta \mid \boldsymbol{x}) d \theta

(14)

事後分布 $f(\theta \mid \boldsymbol{x})$ におけるモデル分布 $f(x^* \mid \theta)$ の期待値

変分ベイズ¶

MCMCより高速な近似計算で事後分布を求める方法。

正則な一般化線形モデルなど比較的単純なモデルであれば適するようだが、複雑なモデルになるほど近似誤差（MCMCとの差）が大きくなる

【徹底解説】変分ベイズをはじめからていねいに | Academaid

確率分布 再訪¶