Skip to content
Published on

確率的機械学習とベイズ手法: ベイズ推論からガウス過程、不確実性の定量化まで

Authors

目次

  1. 確率論の基礎とベイズ定理
  2. ベイズ推論: 事前分布と事後分布
  3. MCMC: マルコフ連鎖モンテカルロ
  4. 変分推論とVAE
  5. ガウス過程
  6. ベイズ深層学習と不確実性の定量化
  7. 実践ツール: PyMC、Stan、Pyro
  8. クイズ

確率論の基礎とベイズ定理

ベイズ定理

ベイズ定理は確率的機械学習の根幹です。事象 AA と証拠 BB が与えられたとき:

P(AB)=P(BA)P(A)P(B)P(A|B) = \frac{P(B|A) \cdot P(A)}{P(B)}

機械学習のコンテキストでは、パラメータ θ\theta とデータ DD で表現します:

P(θD)=P(Dθ)P(θ)P(D)P(\theta|D) = \frac{P(D|\theta) \cdot P(\theta)}{P(D)}

  • P(θ)P(\theta): 事前分布(Prior) — データ観測前のパラメータに関する信念
  • P(Dθ)P(D|\theta): 尤度(Likelihood) — パラメータ値が与えられたときのデータ生成確率
  • P(θD)P(\theta|D): 事後分布(Posterior) — データ観測後に更新された信念
  • P(D)P(D): 周辺尤度(Marginal Likelihood) — 正規化定数

主要な確率分布

ガウス分布(Gaussian Distribution)

N(xμ,σ2)=12πσ2exp ⁣((xμ)22σ2)\mathcal{N}(x | \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)

連続データモデリングの基本。中心極限定理により実務で広く使用されます。

ベルヌーイ分布(Bernoulli Distribution)

Bern(xp)=px(1p)1x,x{0,1}\text{Bern}(x | p) = p^x (1-p)^{1-x}, \quad x \in \{0, 1\}

二値分類問題の基本分布。ベイズ二値分類における尤度関数として使用されます。

ディリクレ分布(Dirichlet Distribution)

Dir(pα)=Γ(kαk)kΓ(αk)kpkαk1\text{Dir}(\mathbf{p} | \boldsymbol{\alpha}) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)} \prod_k p_k^{\alpha_k - 1}

多項分布の共役事前分布(Conjugate Prior)。トピックモデリング(LDA)で重要な役割を果たします。

import numpy as np
import scipy.stats as stats

# 主要な確率分布の確認
x = np.linspace(-4, 4, 200)

# ガウス分布
gaussian = stats.norm(loc=0, scale=1)
print(f"ガウス分布の平均: {gaussian.mean():.2f}, 分散: {gaussian.var():.2f}")

# ベルヌーイ分布
p = 0.7
bern_samples = stats.bernoulli(p).rvs(1000)
print(f"ベルヌーイ成功率: {bern_samples.mean():.3f} (期待値: {p})")

# ディリクレ分布のサンプリング
alpha = np.array([2.0, 3.0, 5.0])
dirichlet_samples = stats.dirichlet(alpha).rvs(1000)
print(f"ディリクレサンプル平均: {dirichlet_samples.mean(axis=0)}")
print(f"理論的平均: {alpha / alpha.sum()}")

ベイズ推論: 事前分布と事後分布

MAP vs MLE

最尤推定(MLE, Maximum Likelihood Estimation)

θ^MLE=argmaxθP(Dθ)=argmaxθlogP(Dθ)\hat{\theta}_{MLE} = \arg\max_\theta P(D|\theta) = \arg\max_\theta \log P(D|\theta)

最大事後確率推定(MAP, Maximum A Posteriori)

θ^MAP=argmaxθP(θD)=argmaxθ[logP(Dθ)+logP(θ)]\hat{\theta}_{MAP} = \arg\max_\theta P(\theta|D) = \arg\max_\theta \left[\log P(D|\theta) + \log P(\theta)\right]

MAPはMLEに事前分布の対数を加えたものと同じです。ガウス事前分布 P(θ)=N(0,τ2)P(\theta) = \mathcal{N}(0, \tau^2) を使用すると:

logP(θ)=θ22τ2+const\log P(\theta) = -\frac{\theta^2}{2\tau^2} + \text{const}

これはL2正則化(Ridge Regression)と数学的に同一です。ラプラス事前分布を使用するとL1正則化(Lasso)に対応します。

共役事前分布(Conjugate Prior)

事後分布が事前分布と同じ関数形を持つ場合、共役関係と言います。

尤度共役事前分布事後分布
ベルヌーイ/二項ベータ(Beta)ベータ(Beta)
ガウスガウスガウス
ポアソンガンマ(Gamma)ガンマ(Gamma)
多項ディリクレディリクレ

ベータ-ベルヌーイ共役の例:

事前: P(p)=Beta(α,β)P(p) = \text{Beta}(\alpha, \beta)、尤度: P(Dp)=Bern(p)NP(D|p) = \text{Bern}(p)^N

事後: P(pD)=Beta(α+nH,β+nT)P(p|D) = \text{Beta}(\alpha + n_H, \beta + n_T)

ここで nHn_H は表の回数、nTn_T は裏の回数です。

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

# PyMCによるベイズ線形回帰
np.random.seed(42)
X = np.random.randn(100)
true_alpha, true_beta, true_sigma = 1.5, 2.3, 0.5
y = true_alpha + true_beta * X + np.random.randn(100) * true_sigma

with pm.Model() as linear_model:
    # 事前分布の定義
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # 尤度関数
    mu = alpha + beta * X
    likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y)

    # NUTSサンプラーで事後分布を推定
    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

# 結果のサマリー
summary = az.summary(trace, var_names=['alpha', 'beta', 'sigma'])
print(summary)
# alpha ~ 1.5, beta ~ 2.3, sigma ~ 0.5 付近に収束するはず

MCMC: マルコフ連鎖モンテカルロ

MCMCは解析的に計算困難な事後分布からサンプルを抽出する手法です。

Metropolis-Hastingsアルゴリズム

  1. 現在の状態 θ(t)\theta^{(t)} から提案分布 q(θθ(t))q(\theta'|\theta^{(t)}) で候補 θ\theta' を生成
  2. 受容確率を計算:

α=min ⁣(1,P(θD)q(θ(t)θ)P(θ(t)D)q(θθ(t)))\alpha = \min\!\left(1, \frac{P(\theta'|D) \cdot q(\theta^{(t)}|\theta')}{P(\theta^{(t)}|D) \cdot q(\theta'|\theta^{(t)})}\right)

  1. uU(0,1)u \sim U(0,1) をサンプリングし、u<αu < \alpha なら θ(t+1)=θ\theta^{(t+1)} = \theta'、そうでなければ θ(t+1)=θ(t)\theta^{(t+1)} = \theta^{(t)}
import numpy as np

def metropolis_hastings(log_posterior, initial, n_samples=10000, proposal_std=0.5):
    """Metropolis-Hastings MCMC実装"""
    samples = np.zeros((n_samples, len(initial)))
    current = np.array(initial, dtype=float)
    accepted = 0

    for i in range(n_samples):
        # 提案分布から候補をサンプリング
        proposal = current + np.random.randn(len(current)) * proposal_std

        # 受容確率の計算(対数空間で)
        log_ratio = log_posterior(proposal) - log_posterior(current)
        accept_prob = min(1.0, np.exp(log_ratio))

        # 受容/棄却の決定
        if np.random.rand() < accept_prob:
            current = proposal
            accepted += 1

        samples[i] = current

    print(f"受容率: {accepted / n_samples:.3f}")
    return samples

# ガウス混合分布からのサンプリング
def log_posterior_mixture(theta):
    log_p1 = -0.5 * ((theta[0] - 2)**2 + (theta[1] - 2)**2)
    log_p2 = -0.5 * ((theta[0] + 2)**2 + (theta[1] + 2)**2)
    return np.log(0.5 * np.exp(log_p1) + 0.5 * np.exp(log_p2) + 1e-10)

samples = metropolis_hastings(log_posterior_mixture, [0.0, 0.0], n_samples=50000)

# バーンイン除去後に分析
burn_in = 5000
posterior_samples = samples[burn_in:]
print(f"事後平均: {posterior_samples.mean(axis=0)}")
print(f"事後標準偏差: {posterior_samples.std(axis=0)}")

ギブスサンプリング

各パラメータを他のパラメータを条件として順次サンプリングします:

θi(t+1)P(θiθi(t),D)\theta_i^{(t+1)} \sim P(\theta_i | \theta_{-i}^{(t)}, D)

条件付き分布が簡単にサンプリングできる場合に特に効率的です。

HMC (Hamiltonian Monte Carlo)

物理学のハミルトン力学を活用して、高次元空間を効率的に探索します。

位置 θ\theta に運動量 rr を導入して、ハミルトニアン:

H(θ,r)=logP(θD)+12rTM1rH(\theta, r) = -\log P(\theta|D) + \frac{1}{2} r^T M^{-1} r

リープフロッグ積分器で軌跡をシミュレーションして、無相関なサンプルを効率的に生成します。PyMCのNUTS(No-U-Turn Sampler)はHMCの自動化バージョンです。


変分推論とVAE

ELBO (Evidence Lower BOund)

事後分布 P(θD)P(\theta|D) が直接計算困難な場合、近似分布 q(θ)q(\theta) を最適化します。

KLダイバージェンスを最小化:

KL[q(θ)P(θD)]=Eq[logq(θ)]Eq[logP(θD)]\text{KL}[q(\theta) \| P(\theta|D)] = \mathbb{E}_{q}[\log q(\theta)] - \mathbb{E}_{q}[\log P(\theta|D)]

整理すると:

logP(D)=ELBO(q)+KL[q(θ)P(θD)]\log P(D) = \text{ELBO}(q) + \text{KL}[q(\theta) \| P(\theta|D)]

ELBO(q)=Eq[logP(Dθ)]KL[q(θ)P(θ)]\text{ELBO}(q) = \mathbb{E}_{q}[\log P(D|\theta)] - \text{KL}[q(\theta) \| P(\theta)]

ELBOを最大化することは、周辺尤度の下界を最大化することと等価です。

平均場近似(Mean-Field Approximation)

パラメータが独立と仮定:

q(θ)=iqi(θi)q(\theta) = \prod_i q_i(\theta_i)

qiq_i に対する最適解:

logqi(θi)=Eqi[logP(D,θ)]+const\log q_i^*(\theta_i) = \mathbb{E}_{q_{-i}}[\log P(D, \theta)] + \text{const}

VAE (Variational Autoencoder)

VAEは変分推論を深層学習に適用した生成モデルです。

VAEのELBO:

L=Eqϕ(zx)[logpθ(xz)]KL[qϕ(zx)p(z)]\mathcal{L} = \mathbb{E}_{q_\phi(z|x)}[\log p_\theta(x|z)] - \text{KL}[q_\phi(z|x) \| p(z)]

  • 第1項: 再構成損失 (Reconstruction Loss)
  • 第2項: 正則化項 (KLダイバージェンス)

再パラメータ化トリック(Reparameterization Trick):

zN(μ,σ2)z \sim \mathcal{N}(\mu, \sigma^2) から直接サンプリングすると逆伝播が不可能です。代わりに:

z=μ+σϵ,ϵN(0,I)z = \mu + \sigma \odot \epsilon, \quad \epsilon \sim \mathcal{N}(0, I)

import torch
import torch.nn as nn

class VAE(nn.Module):
    def __init__(self, input_dim=784, latent_dim=20, hidden_dim=400):
        super().__init__()
        # エンコーダ
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU()
        )
        self.fc_mu = nn.Linear(hidden_dim, latent_dim)
        self.fc_logvar = nn.Linear(hidden_dim, latent_dim)

        # デコーダ
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
            nn.Sigmoid()
        )

    def encode(self, x):
        h = self.encoder(x)
        return self.fc_mu(h), self.fc_logvar(h)

    def reparameterize(self, mu, logvar):
        """再パラメータ化トリック: z = mu + eps * std"""
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        recon = self.decoder(z)
        return recon, mu, logvar

def vae_loss(recon_x, x, mu, logvar):
    # 再構成損失 (BCE)
    recon_loss = nn.functional.binary_cross_entropy(recon_x, x, reduction='sum')
    # KLダイバージェンス: -0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    kl_div = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return recon_loss + kl_div

ガウス過程

定義とカーネル関数

ガウス過程(GP)は関数の分布です。任意の有限集合の関数値が多変量ガウス分布に従う場合、GPと呼びます:

f(x)GP(m(x),k(x,x))f(x) \sim \mathcal{GP}(m(x), k(x, x'))

  • m(x)m(x): 平均関数 (通常0に設定)
  • k(x,x)k(x, x'): カーネル関数 (共分散関数)

主要なカーネル関数:

RBF (Squared Exponential) カーネル: kSE(x,x)=σf2exp ⁣(xx222)k_{SE}(x, x') = \sigma_f^2 \exp\!\left(-\frac{\|x - x'\|^2}{2\ell^2}\right)

Maternカーネル (nu=5/2): kM52(x,x)=σf2(1+5xx+5xx232)exp ⁣(5xx)k_{M52}(x, x') = \sigma_f^2\left(1 + \frac{\sqrt{5}\|x-x'\|}{\ell} + \frac{5\|x-x'\|^2}{3\ell^2}\right)\exp\!\left(-\frac{\sqrt{5}\|x-x'\|}{\ell}\right)

周期カーネル: kper(x,x)=σf2exp ⁣(2sin2(πxx/p)2)k_{per}(x, x') = \sigma_f^2 \exp\!\left(-\frac{2\sin^2(\pi|x-x'|/p)}{\ell^2}\right)

GPR (ガウス過程回帰)

観測データ (X,y)(\mathbf{X}, \mathbf{y}) が与えられたとき、新しい入力 X\mathbf{X}_* に対する予測:

fX,y,XN(fˉ,cov(f))\mathbf{f}_* | \mathbf{X}, \mathbf{y}, \mathbf{X}_* \sim \mathcal{N}(\bar{\mathbf{f}}_*, \text{cov}(\mathbf{f}_*))

fˉ=K(X,X)[K(X,X)+σn2I]1y\bar{\mathbf{f}}_* = K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 I]^{-1}\mathbf{y}

cov(f)=K(X,X)K(X,X)[K(X,X)+σn2I]1K(X,X)\text{cov}(\mathbf{f}_*) = K(\mathbf{X}_*, \mathbf{X}_*) - K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 I]^{-1}K(\mathbf{X}, \mathbf{X}_*)

import torch
import gpytorch

class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# データの準備
torch.manual_seed(42)
train_x = torch.linspace(0, 1, 20)
train_y = torch.sin(train_x * 2 * 3.14159) + torch.randn(20) * 0.1

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

# ハイパーパラメータの最適化 (MLL最大化)
model.train()
likelihood.train()
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(200):
    optimizer.zero_grad()
    output = model(train_x)
    loss = -mll(output, train_y)
    loss.backward()
    optimizer.step()

# 予測
model.eval()
likelihood.eval()
test_x = torch.linspace(0, 1, 100)
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    predictions = likelihood(model(test_x))
    mean = predictions.mean
    lower, upper = predictions.confidence_region()

print(f"学習済みlengthscale: {model.covar_module.base_kernel.lengthscale.item():.4f}")
print(f"学習済みnoise: {likelihood.noise.item():.4f}")

ベイズ深層学習と不確実性の定量化

不確実性の種類

  • 認識論的不確実性(Epistemic Uncertainty): モデル自体の不確実性。データが増えると減少。
  • 偶然論的不確実性(Aleatoric Uncertainty): データに固有のノイズ。データをいくら増やしても減少しない。

MCドロップアウト

テスト時にもドロップアウトを有効にしてベイズ推論を近似します。Gal & Ghahramani (2016)はドロップアウトがガウス過程の近似と数学的に同等であることを証明しました。

P(yx,X,Y)1Tt=1TP(yx,ω^t)P(y^*|x^*, X, Y) \approx \frac{1}{T} \sum_{t=1}^{T} P(y^*|x^*, \hat{\omega}_t)

import torch
import torch.nn as nn
import numpy as np

class BayesianMLP(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim, dropout_p=0.1):
        super().__init__()
        self.dropout_p = dropout_p
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, output_dim)
        )

    def forward(self, x):
        return self.net(x)

def mc_dropout_predict(model, x, n_samples=100):
    """MCドロップアウトで予測不確実性を推定"""
    model.train()  # ドロップアウトを有効に維持
    predictions = []
    with torch.no_grad():
        for _ in range(n_samples):
            pred = model(x)
            predictions.append(pred)

    predictions = torch.stack(predictions)  # (n_samples, batch, output)
    mean = predictions.mean(dim=0)
    # 不確実性 = 予測値の分散
    uncertainty = predictions.var(dim=0)
    return mean, uncertainty

# 使用例
model = BayesianMLP(input_dim=10, hidden_dim=64, output_dim=1, dropout_p=0.1)
x_test = torch.randn(32, 10)
mean_pred, uncertainty = mc_dropout_predict(model, x_test, n_samples=200)
print(f"予測平均のshape: {mean_pred.shape}, 不確実性: {uncertainty.shape}")
print(f"平均不確実性: {uncertainty.mean().item():.4f}")

ベイズニューラルネットワーク(BNN)

重みに分布を付与して完全なベイズ処理を行います:

WN(0,σw2I)W \sim \mathcal{N}(0, \sigma_w^2 I)

変分推論で事後分布を近似:

qϕ(W)=N(μϕ,diag(σϕ2))q_\phi(W) = \mathcal{N}(\mu_\phi, \text{diag}(\sigma_\phi^2))

import torch
import pyro
import pyro.distributions as dist
from pyro.nn import PyroModule, PyroSample

class BayesianLinear(PyroModule):
    def __init__(self, in_features, out_features):
        super().__init__()
        self.in_features = in_features
        self.out_features = out_features

        # 重みとバイアスに事前分布を定義
        self.weight = PyroSample(
            dist.Normal(0., 1.).expand([out_features, in_features]).to_event(2)
        )
        self.bias = PyroSample(
            dist.Normal(0., 10.).expand([out_features]).to_event(1)
        )

    def forward(self, x):
        return x @ self.weight.T + self.bias

class BayesianNN(PyroModule):
    def __init__(self, in_dim, hidden_dim, out_dim):
        super().__init__()
        self.layer1 = BayesianLinear(in_dim, hidden_dim)
        self.layer2 = BayesianLinear(hidden_dim, out_dim)

    def forward(self, x, y=None):
        x = torch.relu(self.layer1(x))
        mu = self.layer2(x).squeeze(-1)
        sigma = pyro.sample("sigma", dist.Uniform(0., 1.))
        with pyro.plate("data", x.shape[0]):
            obs = pyro.sample("obs", dist.Normal(mu, sigma), obs=y)
        return mu

実践ツール: PyMC、Stan、Pyro

NumPyroによるベイズロジスティック回帰

import jax
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS

def logistic_regression_model(X, y=None):
    """NumPyroベイズロジスティック回帰モデル"""
    n_features = X.shape[1]

    # 事前分布
    alpha = numpyro.sample("alpha", dist.Normal(0., 10.))
    beta = numpyro.sample("beta", dist.Normal(jnp.zeros(n_features), jnp.ones(n_features)))

    # ロジットの計算
    logits = alpha + X @ beta

    # 尤度
    with numpyro.plate("obs", X.shape[0]):
        numpyro.sample("y", dist.Bernoulli(logits=logits), obs=y)

# NUTSサンプラーで推論を実行
def run_mcmc(X_train, y_train, num_samples=2000, num_warmup=1000):
    kernel = NUTS(logistic_regression_model)
    mcmc = MCMC(kernel, num_warmup=num_warmup, num_samples=num_samples)
    rng_key = jax.random.PRNGKey(42)
    mcmc.run(rng_key, X_train, y_train)
    return mcmc.get_samples()

# 使用例 (データ生成)
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler

X, y = make_classification(n_samples=200, n_features=4, random_state=42)
X = StandardScaler().fit_transform(X)
X_jax = jnp.array(X)
y_jax = jnp.array(y)

samples = run_mcmc(X_jax, y_jax)
print(f"推定されたalpha: {samples['alpha'].mean():.3f} +/- {samples['alpha'].std():.3f}")
print(f"推定されたbeta: {samples['beta'].mean(axis=0)}")

PyMCによる階層的ベイズモデル

import pymc as pm
import numpy as np

# 階層的ベイズモデル (Hierarchical Bayesian Model)
# 複数グループのデータを同時にモデリング
n_groups = 5
n_per_group = 20
np.random.seed(42)

group_means = np.random.randn(n_groups) * 2
y_obs = np.concatenate([
    np.random.randn(n_per_group) + group_means[g]
    for g in range(n_groups)
])
group_idx = np.repeat(np.arange(n_groups), n_per_group)

with pm.Model() as hierarchical_model:
    # ハイパー事前分布 (Hyperprior)
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_global = pm.HalfNormal('sigma_global', sigma=5)

    # グループ別平均 (部分プーリング, Partial Pooling)
    mu_group = pm.Normal('mu_group', mu=mu_global, sigma=sigma_global, shape=n_groups)

    # 観測尤度
    sigma_obs = pm.HalfNormal('sigma_obs', sigma=2)
    y = pm.Normal('y', mu=mu_group[group_idx], sigma=sigma_obs, observed=y_obs)

    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

import arviz as az
print(az.summary(trace, var_names=['mu_group', 'mu_global', 'sigma_global']))

確率的プログラミングライブラリの比較

ライブラリバックエンド推論手法最適用途
PyMCJAX/NumPyNUTS, VI, SMC柔軟なベイズモデル
StanC++HMC/NUTS研究グレードの推論
PyroPyTorchSVI, MCMC深層確率モデル
NumPyroJAXNUTS, SVI高速GPU推論
TensorFlow ProbabilityTFHMC, VI本番TFパイプライン

クイズ

Q1. MAP推定とMLEの違い、およびMAPが正則化(regularization)と同等である理由は?

答え: MAP = MLE + log事前分布

解説: MLEは argmaxP(Dθ)\arg\max P(D|\theta) を最適化し、MAPは argmaxP(Dθ)P(θ)\arg\max P(D|\theta)P(\theta) を最適化します。対数を取ると: MAP = 対数尤度 + 対数事前分布です。ガウス事前分布 P(θ)exp(λθ2)P(\theta) \propto \exp(-\lambda\|\theta\|^2) を使用すると、対数事前分布は λθ2-\lambda\|\theta\|^2 となり、L2正則化(Ridge)と同一のペナルティを与えます。ラプラス事前分布を使用すると λθ1-\lambda\|\theta\|_1 となり、L1正則化(Lasso)と同一です。つまり、正則化はパラメータに関するベイズ的な事前信念を頻度論的に表現したものです。

Q2. MCMCでバーンイン期間(burn-in period)が必要な理由は?

答え: 初期値への依存性を除去するために必要です。

解説: MCMCチェーンの初期状態は任意に設定されており、事後分布から遠い空間にある可能性があります。マルコフ連鎖が定常分布(stationary distribution)である事後分布に収束(混合)するまでに一定の時間が必要です。バーンイン期間はこの収束前の初期サンプルを破棄するプロセスです。バーンインなしで初期サンプルを含めると、推論が初期値に偏る(バイアス)ことになります。NUTSなどの現代的なサンプラーでは、ウォームアップ段階がこの役割を担います。

Q3. VAEの再パラメータ化トリック(reparameterization trick)が必要な理由は?

答え: 逆伝播可能な勾配を作るためです。

解説: VAEの潜在変数 zN(μϕ(x),σϕ2(x))z \sim \mathcal{N}(\mu_\phi(x), \sigma_\phi^2(x)) を直接サンプリングすると、サンプリング操作は微分不可能な確率ノードになります。逆伝播時に z/ϕ\partial z / \partial \phi を計算できません。再パラメータ化トリックは z=μϕ(x)+σϕ(x)ϵz = \mu_\phi(x) + \sigma_\phi(x) \odot \epsilonϵN(0,I)\epsilon \sim \mathcal{N}(0, I) と変換して確率性を入力 ϵ\epsilon に分離します。これにより zz から μϕ\mu_\phiσϕ\sigma_\phi までの決定論的な経路が生まれ、逆伝播が可能になります。

Q4. ガウス過程でカーネル関数の選択が予測に与える影響は?

答え: カーネルは関数空間における事前信念をエンコードします。

解説: カーネル関数 k(x,x)k(x, x') は入力間の類似度を定義し、学習する関数の構造的な仮定を反映します。RBFカーネルは無限回微分可能な滑らかな関数を仮定し、長さスケール \ell が大きいほど遠い点も相関が高くなります。Maternカーネルは有限の微分可能性を仮定し、より粗い関数を表現します。周期カーネルは周期的なパターンを仮定します。カーネルのハイパーパラメータは周辺尤度(Marginal Likelihood)の最大化によって自動学習されます。不適切なカーネルの選択は過小適合や過大適合につながります。

Q5. MCドロップアウトがベイズ推論の近似となる数学的根拠は?

答え: ドロップアウトは深層ガウス過程に対する変分推論と数学的に同等です。

解説: Gal & Ghahramani (2016)は、任意の深さと非線形性を持つニューラルネットワークにドロップアウトを適用することが、深層ガウス過程に対する変分推論と数学的に同等であることを証明しました。具体的には、各レイヤーの重みにベルヌーイ分布を事前分布として設定して変分推論を実行すると、ドロップアウト学習の目的関数と一致します。したがって、テスト時にドロップアウトを有効にしたまま TT 回順伝播すると、事後分布のモンテカルロ積分になります。ドロップアウト率はモデルの不確実性を調整するハイパーパラメータとしてチューニングが必要です。