Mímisbrunnr知恵の泉

← ベイズ統計 一覧

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：平均場近似と座標上昇法　|　関連：オートエンコーダとVAE（機械学習）

要点（BLUF）

1. なぜ確率的・勾配ベースにするのか

CAVI（平均場近似と座標上昇法）は共役モデルで因子が標準分布になる場合に限り、きれいに回りました。しかし非共役モデルや、変分分布 qq を全共分散ガウスなど自由に取りたい場合は、CAVI の閉じた更新が書けません。そこで ELBO を直接、勾配上昇で最大化します。

∇ϕ ELBO(ϕ)=∇ϕ Eqϕ(θ)[log⁡p(θ,D)−log⁡qϕ(θ)]\nabla_\phi\,\mathrm{ELBO}(\phi)=\nabla_\phi\,\mathbb E_{q_\phi(\theta)}\big[\log p(\theta,D)-\log q_\phi(\theta)\big]

問題は、期待値を取る分布 qϕq_\phi 自体がパラメータ ϕ\phi に依存すること。素朴に微分できません。さらに大規模データでは log⁡p(D∣θ)=∑ilog⁡p(Di∣θ)\log p(D\mid\theta)=\sum_i\log p(D_i\mid\theta) をミニバッチで近似します（SVI の「確率的」はここ）。

2. 再パラメータ化トリック

鍵は、qϕq_\phi からのサンプルを「パラメータに依らない雑音の決定的変換」で書くこと。ガウス qϕ=N(m,s2)q_\phi=\mathcal N(m,s^2) なら

θ=m+s ε,ε∼N(0,1)\theta=m+s\,\varepsilon,\qquad \varepsilon\sim\mathcal N(0,1)

こうすると期待値が ε\varepsilon（ϕ\phi に依らない）についてになり、勾配を期待値の中に入れられます。

∇ϕ Eqϕ[f(θ)]=Eε[∇ϕf(m+sε)]\nabla_\phi\,\mathbb E_{q_\phi}[f(\theta)]=\mathbb E_{\varepsilon}\big[\nabla_\phi f(m+s\varepsilon)\big]

右辺はモンテカルロで素直に推定でき、ff の勾配（自動微分で得られる）を使うので低分散です。対する**スコア関数（REINFORCE）**法 ∇ϕE[f]=E[f ∇ϕlog⁡qϕ]\nabla_\phi\mathbb E[f]=\mathbb E[f\,\nabla_\phi\log q_\phi] は、ff の値をそのまま掛けるため分散が大きい。

3. コード：再パラメータ化で事後を回収する

θ∼N(0,1)\theta\sim\mathcal N(0,1)、yi∼N(θ,σ2)y_i\sim\mathcal N(\theta,\sigma^2) の事後（閉形式 N(1.024,0.3332)\mathcal N(1.024,0.333^2)）を、再パラメータ化トリックの勾配上昇で近似します。

import numpy as np

rng = np.random.default_rng(0)
sigma = 1.0
y = rng.normal(0.8, sigma, size=8); n = len(y); ybar = y.mean()
prec_n = 1.0 + n/sigma**2
mu_n, s_n = (n/sigma**2*ybar)/prec_n, np.sqrt(1/prec_n)    # 閉形式の事後

def grad_logp(th):  return -th + n*(ybar - th)/sigma**2     # ∇_θ[logprior+loglik]

# 再パラメータ化 θ=m+s·ε で ELBO を勾配上昇（log_s で s>0 を保つ）
m, log_s, lr, B = 0.0, 0.0, 0.02, 64
for t in range(4000):
    s = np.exp(log_s)
    eps = rng.standard_normal(B); th = m + s*eps
    g = grad_logp(th)
    g_m = g.mean()                                          # ∂ELBO/∂m
    g_s = (g*eps).mean() + 1/s                              # ∂ELBO/∂s（+エントロピー 1/s）
    m += lr*g_m
    log_s += lr*(g_s*s)                                     # ∂ELBO/∂log_s = ∂ELBO/∂s·s
print(f"再パラメータ化VI: q=N({m:.4f}, {np.exp(log_s):.4f})  / 閉形式 N({mu_n:.4f}, {s_n:.4f})")

出力：

再パラメータ化VI: q=N(1.0169, 0.3334)  / 閉形式 N(1.0240, 0.3333)

出力の意味：勾配上昇だけで、変分パラメータ (m,s)(m,s) が閉形式の事後 N(1.024,0.333)\mathcal N(1.024,0.333) にほぼ到達しました（平均のわずかなずれは確率的勾配の残差）。共役の更新式を一切使わず、∇θlog⁡p(θ,D)\nabla_\theta\log p(\theta,D) さえ計算できれば事後に届く——この汎用性が SVI の価値です。実務では ∇θlog⁡p\nabla_\theta\log p を自動微分で得るので、モデルを書くだけで動きます。

4. 再パラメータ化 vs スコア関数：勾配の分散

なぜ再パラメータ化が好まれるか。同じ ∂ELBO/∂m\partial\mathrm{ELBO}/\partial m を、両方の方法の単サンプル推定で比べ、ばらつきを見ます。

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
sigma = 1.0
y = rng.normal(0.8, sigma, size=8); n = len(y); ybar = y.mean()
grad_logp = lambda th: -th + n*(ybar - th)/sigma**2
logp = lambda th: -0.5*th**2 - 0.5*np.sum((y[:,None]-th[None,:])**2,axis=0)/sigma**2

rng2 = np.random.default_rng(1)
m0, s0 = 0.5, 0.5
eps = rng2.standard_normal(20000)
g_reparam = grad_logp(m0 + s0*eps)                          # 再パラメータ化
th2 = rng2.normal(m0, s0, 20000)
g_score = (logp(th2) - stats.norm(m0,s0).logpdf(th2)) * (th2 - m0)/s0**2   # スコア関数
print(f"∂ELBO/∂m の単サンプル推定の標準偏差:")
print(f"  再パラメータ化 : {g_reparam.std():.3f} (平均{g_reparam.mean():+.3f})")
print(f"  スコア関数     : {g_score.std():.3f} (平均{g_score.mean():+.3f})")
print(f"  → 再パラメータ化の分散は {g_score.std()/g_reparam.std():.1f} 倍小さい")

出力：

∂ELBO/∂m の単サンプル推定の標準偏差:
  再パラメータ化 : 4.473 (平均+4.767)
  スコア関数     : 11.738 (平均+4.679)
  → 再パラメータ化の分散は 2.6 倍小さい

出力の意味：両者の勾配の平均はほぼ同じ（+4.77+4.77 と +4.68+4.68、どちらも不偏）ですが、標準偏差は再パラメータ化 4.474.47 に対しスコア関数 11.7411.74——2.62.6 倍も小さい。低分散な勾配は学習を速く安定させます。高次元になるほどこの差は効き、再パラメータ化が現代の確率的 VI の標準になっている理由です。

5. VAE への橋渡し

再パラメータ化トリックは、**変分オートエンコーダ（VAE）**の学習そのものです（VAEとELBOの共通構造）。VAE では各データ点の潜在変数 zz の変分事後 qϕ(z∣x)=N(μϕ(x),σϕ(x)2)q_\phi(z\mid x)=\mathcal N(\mu_\phi(x),\sigma_\phi(x)^2) をニューラルネットで出し、z=μϕ(x)+σϕ(x)εz=\mu_\phi(x)+\sigma_\phi(x)\varepsilon と再パラメータ化して ELBO を勾配上昇します。エンコーダ・デコーダの重みを自動微分で更新できるのは、この trick のおかげです（実装フレームワークの最新 API は確率的プログラミング概観で、要最新確認）。

⚠️ よくある誤解

関連ノート