- Authors

- Name
- Youngju Kim
- @fjvbn20031
目次
- 確率論の基礎とベイズ定理
- ベイズ推論: 事前分布と事後分布
- MCMC: マルコフ連鎖モンテカルロ
- 変分推論とVAE
- ガウス過程
- ベイズ深層学習と不確実性の定量化
- 実践ツール: PyMC、Stan、Pyro
- クイズ
確率論の基礎とベイズ定理
ベイズ定理
ベイズ定理は確率的機械学習の根幹です。事象 と証拠 が与えられたとき:
機械学習のコンテキストでは、パラメータ とデータ で表現します:
- : 事前分布(Prior) — データ観測前のパラメータに関する信念
- : 尤度(Likelihood) — パラメータ値が与えられたときのデータ生成確率
- : 事後分布(Posterior) — データ観測後に更新された信念
- : 周辺尤度(Marginal Likelihood) — 正規化定数
主要な確率分布
ガウス分布(Gaussian Distribution)
連続データモデリングの基本。中心極限定理により実務で広く使用されます。
ベルヌーイ分布(Bernoulli Distribution)
二値分類問題の基本分布。ベイズ二値分類における尤度関数として使用されます。
ディリクレ分布(Dirichlet Distribution)
多項分布の共役事前分布(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)
最大事後確率推定(MAP, Maximum A Posteriori)
MAPはMLEに事前分布の対数を加えたものと同じです。ガウス事前分布 を使用すると:
これはL2正則化(Ridge Regression)と数学的に同一です。ラプラス事前分布を使用するとL1正則化(Lasso)に対応します。
共役事前分布(Conjugate Prior)
事後分布が事前分布と同じ関数形を持つ場合、共役関係と言います。
| 尤度 | 共役事前分布 | 事後分布 |
|---|---|---|
| ベルヌーイ/二項 | ベータ(Beta) | ベータ(Beta) |
| ガウス | ガウス | ガウス |
| ポアソン | ガンマ(Gamma) | ガンマ(Gamma) |
| 多項 | ディリクレ | ディリクレ |
ベータ-ベルヌーイ共役の例:
事前: 、尤度:
事後:
ここで は表の回数、 は裏の回数です。
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アルゴリズム
- 現在の状態 から提案分布 で候補 を生成
- 受容確率を計算:
- をサンプリングし、 なら 、そうでなければ
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)}")
ギブスサンプリング
各パラメータを他のパラメータを条件として順次サンプリングします:
条件付き分布が簡単にサンプリングできる場合に特に効率的です。
HMC (Hamiltonian Monte Carlo)
物理学のハミルトン力学を活用して、高次元空間を効率的に探索します。
位置 に運動量 を導入して、ハミルトニアン:
リープフロッグ積分器で軌跡をシミュレーションして、無相関なサンプルを効率的に生成します。PyMCのNUTS(No-U-Turn Sampler)はHMCの自動化バージョンです。
変分推論とVAE
ELBO (Evidence Lower BOund)
事後分布 が直接計算困難な場合、近似分布 を最適化します。
KLダイバージェンスを最小化:
整理すると:
ELBOを最大化することは、周辺尤度の下界を最大化することと等価です。
平均場近似(Mean-Field Approximation)
パラメータが独立と仮定:
各 に対する最適解:
VAE (Variational Autoencoder)
VAEは変分推論を深層学習に適用した生成モデルです。
VAEのELBO:
- 第1項: 再構成損失 (Reconstruction Loss)
- 第2項: 正則化項 (KLダイバージェンス)
再パラメータ化トリック(Reparameterization Trick):
から直接サンプリングすると逆伝播が不可能です。代わりに:
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と呼びます:
- : 平均関数 (通常0に設定)
- : カーネル関数 (共分散関数)
主要なカーネル関数:
RBF (Squared Exponential) カーネル:
Maternカーネル (nu=5/2):
周期カーネル:
GPR (ガウス過程回帰)
観測データ が与えられたとき、新しい入力 に対する予測:
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)はドロップアウトがガウス過程の近似と数学的に同等であることを証明しました。
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)
重みに分布を付与して完全なベイズ処理を行います:
変分推論で事後分布を近似:
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']))
確率的プログラミングライブラリの比較
| ライブラリ | バックエンド | 推論手法 | 最適用途 |
|---|---|---|---|
| PyMC | JAX/NumPy | NUTS, VI, SMC | 柔軟なベイズモデル |
| Stan | C++ | HMC/NUTS | 研究グレードの推論 |
| Pyro | PyTorch | SVI, MCMC | 深層確率モデル |
| NumPyro | JAX | NUTS, SVI | 高速GPU推論 |
| TensorFlow Probability | TF | HMC, VI | 本番TFパイプライン |
クイズ
Q1. MAP推定とMLEの違い、およびMAPが正則化(regularization)と同等である理由は?
答え: MAP = MLE + log事前分布
解説: MLEは を最適化し、MAPは を最適化します。対数を取ると: MAP = 対数尤度 + 対数事前分布です。ガウス事前分布 を使用すると、対数事前分布は となり、L2正則化(Ridge)と同一のペナルティを与えます。ラプラス事前分布を使用すると となり、L1正則化(Lasso)と同一です。つまり、正則化はパラメータに関するベイズ的な事前信念を頻度論的に表現したものです。
Q2. MCMCでバーンイン期間(burn-in period)が必要な理由は?
答え: 初期値への依存性を除去するために必要です。
解説: MCMCチェーンの初期状態は任意に設定されており、事後分布から遠い空間にある可能性があります。マルコフ連鎖が定常分布(stationary distribution)である事後分布に収束(混合)するまでに一定の時間が必要です。バーンイン期間はこの収束前の初期サンプルを破棄するプロセスです。バーンインなしで初期サンプルを含めると、推論が初期値に偏る(バイアス)ことになります。NUTSなどの現代的なサンプラーでは、ウォームアップ段階がこの役割を担います。
Q3. VAEの再パラメータ化トリック(reparameterization trick)が必要な理由は?
答え: 逆伝播可能な勾配を作るためです。
解説: VAEの潜在変数 を直接サンプリングすると、サンプリング操作は微分不可能な確率ノードになります。逆伝播時に を計算できません。再パラメータ化トリックは 、 と変換して確率性を入力 に分離します。これにより から 、 までの決定論的な経路が生まれ、逆伝播が可能になります。
Q4. ガウス過程でカーネル関数の選択が予測に与える影響は?
答え: カーネルは関数空間における事前信念をエンコードします。
解説: カーネル関数 は入力間の類似度を定義し、学習する関数の構造的な仮定を反映します。RBFカーネルは無限回微分可能な滑らかな関数を仮定し、長さスケール が大きいほど遠い点も相関が高くなります。Maternカーネルは有限の微分可能性を仮定し、より粗い関数を表現します。周期カーネルは周期的なパターンを仮定します。カーネルのハイパーパラメータは周辺尤度(Marginal Likelihood)の最大化によって自動学習されます。不適切なカーネルの選択は過小適合や過大適合につながります。
Q5. MCドロップアウトがベイズ推論の近似となる数学的根拠は?
答え: ドロップアウトは深層ガウス過程に対する変分推論と数学的に同等です。
解説: Gal & Ghahramani (2016)は、任意の深さと非線形性を持つニューラルネットワークにドロップアウトを適用することが、深層ガウス過程に対する変分推論と数学的に同等であることを証明しました。具体的には、各レイヤーの重みにベルヌーイ分布を事前分布として設定して変分推論を実行すると、ドロップアウト学習の目的関数と一致します。したがって、テスト時にドロップアウトを有効にしたまま 回順伝播すると、事後分布のモンテカルロ積分になります。ドロップアウト率はモデルの不確実性を調整するハイパーパラメータとしてチューニングが必要です。