変分推論２

A Blog Entry on Bayesian Computation by an Applied Mathematician

$$

$$

要約

EM アルゴリズムとは，目的関数が $$h(x)=\mathrm{E}[H(x,Z)]$$ と表せる場合に対する MM アルゴリズムである．欠測モデルや混合モデルにおいて，最尤推定量を求める用途等の応用がある．このような状況では，（対数）尤度に関して下からの評価を与える代理関数で，特に振る舞いのようもの（$Q$-関数）が見つかるため（ Equation 2 ），代わりにこれを逐次的に最大化することが出来るのである．この手続きは変分 Bayes 法の特別な場合にあたる．

1 最尤推定

クラスタリングを一度さっぱり忘れて，最尤推定を思い出してみる．

定義：最尤推定量¹ (Fisher, 1912)

$\{P_{\theta }{\}}_{\theta \in \mathrm{\Theta }}\subset \mathcal{P}(\mathcal{X})$ を統計モデルで，ある共通の $\sigma$-有限測度 $\mu \in \mathcal{P}(\mathcal{X})$ に関して密度 $\{p_{\theta }{\}}_{\theta \in \mathrm{\Theta }}$ を持つとする．

独立な観測 $X_1,\cdots ,X_n$ の 最尤推定量 とは，モデルの対数尤度 $$\log p_{\theta }$$ を通じて定まる次の目的関数 ${\ell }_n:\mathrm{\Theta }\to (-\mathrm{\infty },0)$ を最大化するような $M$-推定量 をいう： $${\ell }_n(\theta ;X_1,\cdots ,X_n):=\sum _{i=1}^n\log p_{\theta }(X_i),\theta \in \mathrm{\Theta }.$$

1.1 最尤推定と最適化

すなわち，最尤推定量とは最適化問題の解として定式化されるのである．

最大値点であるということは，停留点である必要があるから，微分が零になるという条件を通じて解析的に求まることもある．この $Z$-推定量としての特徴付けは (Carmer, 1946, p. 498) による．

また，計算機的な方法では，Iterative Propertional Fittting や勾配に基づく最適化手法を用いることも考えられる (Robbins and Monro, 1951), (Fletcher, 1987)．

最尤推定量が解析的に求まらない場面には，代表的には欠測モデルなどがある．欠測モデルは，観測される確率変数 $X$ の他に，観測されない確率変数 $Z$ も想定し，その同時分布を考えるモデルである．これにより，$(X,Z)$ 全体には単純な仮定しか置かずとも，$X$ に対して複雑な分布を想定することが可能になるのである．

この場合には，モデルの構造を利用して最尤推定量を求めるための MM アルゴリズム (Sun et al., 2016) の例がある．これが EM アルゴリズム (Dempster et al., 1977) である．

現在でも，その他の MM アルゴリズムが，種々の最適化問題に対する “problem-driven algorithm” であり続けている (T. T. Wu and Lange, 2010)．

1.2 最尤推定と Bayes 推定

最尤推定は，一様事前分布をおいた場合の MAP 推定 とみなせる．この意味で，Bayes 推定の特殊な場合である．

Bayes 推定は MCMC や SMC などのサンプリング法によって統一的に行えるが，殊に MAP 推定に対しては，効率的な最適化法として EM アルゴリズムが使える，ということである．

より一般の Bayes 推定に対応できるような EM アルゴリズムの一般化が，近似アルゴリズムとして存在する．これが次稿で紹介する 変分推論 である．

2 EM アルゴリズム

EM アルゴリズムは，混合モデルに対する最尤推定アルゴリズムである．一般に，目的関数が $$h(x)=\mathrm{E}[H(x,Z)]$$ と表せる場合に対する特殊な MM アルゴリズムである．²

2.1 欠測データと混合モデル

欠測データ (incomplete data) とは，２つの確率変数 $(Z,X)$ について次の図式が成り立つ際の，$Z$ を潜在変数として，$X$ からの観測とみなせるデータをいう (Dempster et al., 1977, p. 1)：

図 1: Missing Data Model / Latent State Model / Completed Model for $X$

これは， $$\begin{array}{}\text{(1)}& p(x|\theta )={\int }_{\mathcal{Z}}p(x,z|\theta )dz\end{array}$$ という形の尤度を持つモデルである．³

これは潜在変数 $Z$ を持つモデルの最も単純な例ともみなせる．特に $Z$ が離散変数である場合，$X$ に対する混合モデルともいう．隠れ Markov モデル はこの発展例である．⁴

このように，$X$ の分布を，潜在変数 $Z$ を追加して理解することを，モデルの 完備化 (completion) または 脱周辺化 (demarginalization)，またはデータの拡張 (data augmentation) ともいう．⁵

2.2 EM アルゴリズム

値域 $\mathcal{Z}$ を持つ潜在変数 $Z$ とパラメータ $\theta \in \mathrm{\Theta }$ に関して Equation 1 で表せる尤度関数 $p(x|\theta )$ に関して，Jensen の不等式より，任意の $x,\theta$ で添字づけられた確率密度関数 $q:\mathcal{Z}\to {\mathbb{R}}_{+}$⁶ とパラメータ $\theta \in \mathrm{\Theta }$ について次の評価が成り立つ：

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\log p(x|\theta )& =\log {\int }_{\mathcal{Z}}p(x,z|\theta )dz\\ & \ge {\int }_{\mathcal{Z}}q(z|x,\phi )\log \frac{p(x,z|\theta )}{q(z|x,\phi )}dz\\ & =:F(q,\theta ).\end{array}\end{array}$$

この事実に基づき，$F$ を代理関数として，これを２つの変数 $q,\theta$ について交互に最大化するという手続きを，EM アルゴリズム という．⁷

1. $E$-ステップ：$F$ を $q$ について最大化する． $$\begin{array}{rl}F(q,\theta )& ={\int }_{\mathcal{Z}}q(z|x,\phi )\log \frac{p(z|x,\theta )p(x|\theta )}{q(z|x,\phi )}dz\\ & =\log p(x|\theta )-\mathrm{KL}(q_{\phi },p_{\theta }).\end{array}$$ より，$q(z|x,\phi )=p(z|x,\theta )$ で最大化される．⁸
2. $M$-ステップ：$F$ を $\theta$ について最大化する． $$\begin{array}{rl}F(q,\theta )& ={\int }_{\mathcal{Z}}q(z|x,\phi )\log p(x,z|\theta )dz\\ & -{\int }_{\mathcal{Z}}q(z|x,\phi )\log q(z|x,\phi )dz\\ & =\underset{=:Q(\theta |\phi ,x)}{\underbrace{(q_{\phi }dz|\log p_{\theta })}}+H(q_{\phi })\end{array}$$ より，$Q$ の停留点で最大化される．

総じて，EM アルゴリズムは $p,q$ の KL 乖離度を逐次的に最小化している．

2.3 $E$-ステップの変形

$M$-ステップにおける $F$ の $\theta$ における最大化は $Q$ の $\theta$ による最大化に等価であるから，$E$-ステップは結局，事後分布 $p(z|x,\theta )$ を計算し，これに関する積分である $$Q(\theta |\phi ,x)={\int }_{\mathcal{Z}}p(z|x,\theta )\log p(x,z|\theta )dz$$ を計算する，というステップになる．

モデル $\{p_{\theta }\}$ を複雑にしすぎた場合，この $Q$ の計算は困難で実行不可能になってしまう．解析的に $E$-ステップを実行したい場合，典型的には指数型分布族を仮定する．

そこで，$Q$ を Monte Carlo 推定量で代替して，それを最大化した場合の EM アルゴリズムを MCEM (Monte Carlo EM) という (Wei and Tanner, 1990a), (Wei and Tanner, 1990b)．典型的には Metropolis-Hastings アルゴリズムを用いることになり (Chau et al., 2021)，これがスケーラビリティ問題を産む．⁹

また，この $E$-ステップで必ずしも完全な最大化を達成する必要はない (Neal and Hinton, 1998), (Bishop, 2006, p. 454)．従って，$p(z|x,\theta )$ が複雑すぎる場合，十分近い $q$ を選択してこれに関する積分として $Q$ を近似することが考えられる．特に $p$ を 変分近似 した場合，変分 EM アルゴリズムという (Wainwright and Jordan, 2008, p. 154)．

2.4 $M$-ステップの変形

$Q$ の停留点を探すにあたって，典型的には微分が消える点を探す．

しかしこれが難しい場合，厳密な最大化は行わず，代わりにせめて「現状よりは大きくする」ことを実行するアルゴリズムを用いた場合，これを 一般化 EM アルゴリズム (GEM: Generalized EM) ともいう (Bishop, 2006, p. 454), (Hastie et al., 2009, p. 277)．

例えば，大域的最大化の代わりに条件付き最大化を行うこととする方法 ECM (Expectation Conditional Maximization) などがその例である (Meng and Rubin, 1991), (Meng and Rubin, 1993)．(Christian P. Robert and Casella, 2004, p. 200) も参照．

2.5 EM アルゴリズムの有効性

命題：尤度は単調減少する (Dempster et al., 1977)¹⁰

$\{{\hat{\theta }}_{(j)}\}\subset \mathrm{\Theta }$ を EM アルゴリズムの $M$-ステップでの出力列とする．このとき，

$$L({\hat{\theta }}_{(j+1)}|x)\ge L({\hat{\theta }}_{(j)}|x).$$ 等号成立は $$Q({\hat{\theta }}_{(j+1)}|{\hat{\theta }}_{(j)},x)=Q({\hat{\theta }}_{(j)}|{\hat{\theta }}_{(j)},x)$$ の場合のみ．

証明

命題：局所解への収束 (Boyles, 1983)-(C. F. J. Wu, 1983)¹¹

$$Q(\theta |{\theta }_0,x):={\int }_{\mathcal{Z}}p(z|\theta ,x)\log p(\theta |x,z)dz$$ は $\theta ,{\theta }_0\in \mathrm{\Theta }$ について連続であるとする．このとき，EM アルゴリズムの出力 $\{{\hat{\theta }}_{(}j)\}$ は尤度 $p(\theta |x)$ の停留点に単調に収束する．

証明

よって，EM アルゴリズムは局所解には収束する．

しかし，常に尤度が単調増加するという性質上，局所解に囚われてしまった場合，そこから逃れることはないことになる．

大域解に収束することを保証したい場合は，異なる初期値で複数回 EM アルゴリズムを実行するか，擬似除冷 (simulated annealing)¹² などの別の手法を用いることを考える必要がある (Finch et al., 1989)．

3 EM アルゴリズムの実装（Gauss 有限混合モデルの場合）

負担率に確率モデルを置いた場合，ソフト $K$-平均アルゴリズム は EM アルゴリズムになる．

EM アルゴリズムは一般に多峰性に弱いことをここで示す．¹³

3.1 Guass 有限混合モデル

ここでは，以下の，有限な混合モデルで，さらに混合される分布は正規であるものを考える：

定義 (Gaussian finite mixture model)

集合 $[K]$ 上に値を取る隠れ変数 $Z$ の確率質量関数を $(p_k{)}_{k=1}^K$ とする． $$\begin{array}{}\text{(3)}& p(x;({\mu }_k),({\mathrm{\Sigma }}_k),(p_k))=\sum _{k=1}^K{p}_k\varphi (x;{\mu }_k,{\mathrm{\Sigma }}_k)\end{array}$$ として定まるモデル $(p_{({\mu }_k),({\mathrm{\Sigma }}_k),(p_k)})$ を ${\mathbb{R}}^d$ 上の Gauss 有限混合モデル という．

ただし，$\varphi (x$; $\mu$, $\sigma )$ は ${\mathrm{N}}_d(\mu ,{\sigma }^2)$ の密度とした．

Equation 3 は $(X,Z)$ 上の結合分布の族を表しており，そのパラメータは $\theta :=(({\mu }_k),({\mathrm{\Sigma }}_k),(p_k))$ である．さらに，${\theta }_k:=({\mu }_k,{\mathrm{\Sigma }}_k)$ と定める．

3.2 Gauss 有限混合モデルでの EM アルゴリズム

Section 2.2 での議論を，今回の Gauss 有限混合モデルに当てはめてみる．

対数周辺尤度は

$$\begin{array}{rl}\log p(x|\theta )& =\log \left(\sum _{k=1}^K{q}_k\frac{p_k\varphi (x|{\theta }_k)}{q_k}\right)\\ & \ge \sum _{k=1}^K{q}_k\log \left(\frac{p_k\varphi (x|{\theta }_k)}{q_k}\right)\\ & =:F(q_k,\theta ).\end{array}$$

という下界を持つ．

これに基づき，観測 $\{x^{(n)}{\}}_{n=1}^N$ と混合 Gauss モデル 3.1 に対する EM アルゴリズムは次の２段階を繰り返す：

1. $E$-ステップ： $$\begin{array}{rl}{r}_k^{(n)}& \leftarrow \mathrm{P}[Z=k|x^{(n)},\theta ]\\ & =\frac{p_k\varphi (x^{(n)}|{\theta }_k)}{\sum _{j=1}^K{p}_j\varphi (x^{(n)}|{\theta }_j)}\end{array}$$ を計算して，$F$ に代入する．
2. $M$-ステップ：$F$ を $\theta$ について最大化する．これは，次の値を計算することに等しい：
   1. $${\mu }_k\leftarrow \frac{\sum _{n=1}^N{r}_k^{(n)}{x}^{(n)}}{\sum _{n=1}^N{r}_k^{(n)}}$$
   2. $${\mathrm{\Sigma }}_k\leftarrow \frac{\sum _{n=1}^N{r}_k^{(n)}(x^{(n)}-{\mu }_k)(x^{(n)}-{\mu }_k{)}^{\mathrm{\top }}}{\sum _{n=1}^N{r}_k^{(n)}}$$
   3. $$p_k\leftarrow \frac{\sum _{n=1}^N{r}_k^{(n)}}{N}$$

$E$-ステップの導出

最初のステップは Bayes の定理から，

$$\begin{array}{rl}\mathrm{P}[Z=k|x,\theta ]& =\frac{p(x|z=k,\theta )p(z=k|\theta )}{p(x|\theta )}\\ & =\frac{p_k\varphi (x|{\theta }_k)}{\sum _{j=1}^K{p}_j\varphi (x|{\theta }_j)}.\end{array}$$ の計算に帰着する．$\{p_k\}$ が一様で，${\sigma }_k=\sigma$ も一様であるとき，これはソフト $K$-平均法における 負担率 $r_{kn}$ に他ならない．このとき， $$\beta =\frac{1}{2{\sigma }^2}.$$

こうして，負担率とは，「データ点がそのクラスターに属するという事後確率」としての意味も持つことが判った．

$M$-ステップの導出¹⁴

2,3,4 はそれぞれ条件 $$\frac{\partial }{\partial {\mu }_k}\sum _{n=1}^N\log p(x_n|\theta )=0$$ $$\frac{\partial }{\partial {\mathrm{\Sigma }}_k}\sum _{n=1}^N\log p(x_n|\theta )=0$$ $$\frac{\partial }{\partial p_k}\sum _{n=1}^N\log p(x_n|\theta )=0$$ から出る．

最大化に Newton-Raphson 法を用いたとも捉えられる．¹⁵

3.3 $K$-平均アルゴリズムとの対応

$E$-ステップが assignment ステップ，$M$-ステップが update ステップに対応する．

ハード $K$-平均法は，歪み尺度 (distortion measure) $$J(r,\mu ):=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}\parallel x_n-{\mu }_k{\parallel }^2$$ を $r,\mu$ のそれぞれについて逐次的に最小化する手法とも見れる．¹⁶

3.4 Gauss 混合モデルの場合¹⁷

$K=2$ での Gauss 混合分布 $$\begin{array}{}\text{(4)}& p\mathrm{N}({\mu }_1,{\sigma }^2)+(1-p)\mathrm{N}({\mu }_2,{\sigma }^2),\end{array}$$ $$p=0.7,\sigma =1,$$ を考える．未知パラメータは $\theta :=({\mu }_1,{\mu }_2)$ である．

実は，混合モデルでは，ここまで単純な例でさえ，尤度は多峰性を持つ．

試しに，$({\mu }_1,{\mu }_2)=(0,3.1)$ として 500 個のデータを生成し，モデル 4 が定める尤度をプロットしてみると，次の通りになる：

<>:42: SyntaxWarning: invalid escape sequence '\m'
<>:43: SyntaxWarning: invalid escape sequence '\m'
<>:44: SyntaxWarning: invalid escape sequence '\m'
<>:42: SyntaxWarning: invalid escape sequence '\m'
<>:43: SyntaxWarning: invalid escape sequence '\m'
<>:44: SyntaxWarning: invalid escape sequence '\m'
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/1951455018.py:42: SyntaxWarning: invalid escape sequence '\m'
  plt.title('Log Likelihood for 500 Observations with $(\mu_1,\mu_2)=(0,3.1)$')
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/1951455018.py:43: SyntaxWarning: invalid escape sequence '\m'
  plt.xlabel('$\mu_1$')
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/1951455018.py:44: SyntaxWarning: invalid escape sequence '\m'
  plt.ylabel('$\mu_2$')

真値 $({\mu }_1,{\mu }_2)=(0,3.1)$ で確かに最大になるが，$({\mu }_1,{\mu }_2)=(2,-0.5)$ 付近で極大値を取っていることがわかる．

3.5 EM アルゴリズムの初期値依存性

EM アルゴリズムはその初期値依存性からランダムな初期値から複数回実行してみる必要がある．モデル 4 の場合，その結果は次のようになる：

class EM_1d:
    """
    Gauss 有限混合モデルに対する EM アルゴリズム

    Parameters:
    - K (int): 混合成分の数．デフォルトは2．
    - max_iter (int): アルゴリズムの最大反復回数．デフォルトは100．
    - tol (float): 収束の閾値．連続する反復での対数尤度の差がこの値以下になった場合，アルゴリズムは収束したと見なされる．デフォルトは1e-4．
    """

    def __init__(self, K=2, init=None, max_iter=100, tol=1e-4):
        self.K = K
        self.max_iter = max_iter
        self.tol = tol

        self.means = None
        self.variances = None
        self.mixing_coefficients = None
        self.log_likelihood_history = []
        self.mean_history = []
        self.initial_value = init

    def expectation(self, X):
        """
        E ステップ

        Parameters:
        - X (ndarray): 観測データ．
        """
        N = X.shape[0]
        r = np.zeros((N, self.K))
        for k in range(self.K):
            pdf = norm.pdf(X, self.means[k], np.sqrt(self.variances[k]))
            r[:, k] = self.mixing_coefficients[k] * pdf
        r /= r.sum(axis=1, keepdims=True)
        return r

    def maximization(self, X, r):
        """
        M ステップ

        Parameters:
        - X (ndarray): 観測データ．
        - r (ndarray): 負担率．
        """
        N = X.shape[0]
        Nk = r.sum(axis=0)
        self.means = (X.T @ r / Nk).T
        self.variances = np.zeros(self.K)
        for k in range(self.K):
            diff = X - self.means[k]
            self.variances[k] = (r[:, k] @ (diff ** 2)) / Nk[k]
        self.mixing_coefficients = Nk / N

    def compute_log_likelihood(self, X):
        """
        対数尤度の計算

        Parameters:
        - X (ndarray): 観測データ．
        """
        log_likelihood = 0
        for x in X:
            log_likelihood += np.log(np.sum([self.mixing_coefficients[k] * norm.pdf(x, self.means[k], np.sqrt(self.variances[k])) for k in range(self.K)]))
        return log_likelihood
    
    def fit(self, X):
        """
        EM アルゴリズムの実行

        Parameters:
        - X (ndarray): 観測データ．
        """
        N = X.shape[0]
        np.random.seed(42)

        if self.initial_value is None:
            random_indeces = np.random.choice(N, self.K, replace=False)
            self.initial_value = X[random_indeces]
        self.means = self.initial_value
        self.initial_value = self.means
        self.variances = np.ones(self.K)
        self.mixing_coefficients = np.ones(self.K) / self.K

        # 反復
        for _ in range(self.max_iter):
            r = self.expectation(X)
            self.maximization(X, r)
            log_likelihood = self.compute_log_likelihood(X)
            self.log_likelihood_history.append(log_likelihood)
            self.mean_history.append(self.means)

            if len(self.log_likelihood_history) >= 2 and np.abs(self.log_likelihood_history[-1] - self.log_likelihood_history[-2]) < self.tol:
                break
        
        return self

<>:32: SyntaxWarning: invalid escape sequence '\m'
<>:33: SyntaxWarning: invalid escape sequence '\m'
<>:34: SyntaxWarning: invalid escape sequence '\m'
<>:32: SyntaxWarning: invalid escape sequence '\m'
<>:33: SyntaxWarning: invalid escape sequence '\m'
<>:34: SyntaxWarning: invalid escape sequence '\m'
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/163584773.py:32: SyntaxWarning: invalid escape sequence '\m'
  axs[1].set_title('Mean Values Progress ($\mu_1,\mu_2$)')
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/163584773.py:33: SyntaxWarning: invalid escape sequence '\m'
  axs[1].set_xlabel('$\mu_1$')
/var/folders/gx/6w78f6997l5___173r25fp3m0000gn/T/ipykernel_74977/163584773.py:34: SyntaxWarning: invalid escape sequence '\m'
  axs[1].set_ylabel('$\mu_2$')

4 Monte Carlo 法による解決

Equation 1 の逆向きの関係 $$\begin{array}{rl}p(\theta |x)& ={\int }_{\mathcal{Z}}p(z,\theta |x)dz\\ & ={\int }_{\mathcal{Z}}p(\theta |z,x)p(z|x)dz\end{array}$$ も成り立つという 階層構造 (hierarchical structure) を持つモデルにおいて，Bayes 推論が Gibbs サンプラーによって実行できる (C. P. Robert, 1996)．¹⁸

このような欠測モデルの文脈で Gibbs サンプラーを用いる手法は，データ拡張 の名前でも知られる (Tanner and Wong, 1987)．

加えて，初期値依存性や局所解へのトラップが懸念されるという EM アルゴリズムの問題点を，MCMC はいずれも持ち合わせていない．

さらに，混合数 $K$ に関する検定も構成できる (Mengersen and Robert, 1996) など，Gibbs サンプラーひとつで確率モデルに関する種々の情報を取り出せる．

最尤推定の代わりに Bayes 推定を行なっているため，データ数が少なくとも，過学習の問題が起こりにくいという利点もある．

Bayes 階層モデルは複雑なモデルに対する表現力が高く，地球科学をはじめとして多くの応用分野で使われている (Hrafnkelsson, 2023)．

4.1 Gibbs サンプリング

高次元な確率変数 $(U_1,\cdots ,U_K)$ のシミュレーションを行いたい場合，直接行うのではなく，条件付き分布 $p(u_k|u_{-k})$ からのサンプリングを繰り返すことでこれを行うことが出来る．¹⁹

1. 任意の初期値 $U_1^{(0)},\cdots ,U_K^{(0)}$ を与える．
2. 各 $k\in [K]$ について， $$U_k^{(t)}\sim p(u_k|U_{-k}^{(t-1)})$$ をサンプリングする．
3. 十分時間が経過した際，アルゴリズムの出力 $(U_1^{(t)},\cdots ,U_K^{(t)})$ は $(U_1,\cdots ,U_K)$ と同分布になる．

実際，$\{(U_1^{(t)},\cdots ,U_K^{(t)}){\}}_{t\in \mathbb{N}}$ はエルゴード的な Markov 連鎖を定め，定常分布 $p(U_1,\cdots ,U_K)$ を持つ．

命題 (Jean Diebolt and Robert, 1994)²⁰

$p(u_1|u_2)$ が正，または $p(u_2|u_1)$ が正ならば，Markov 連鎖 $\{U_1^{(t)}\},\{U_2^{(t)}\}$ はいずれもエルゴード的で，不変分布 $p(u_1|u_2),p(u_2|u_1)$ を持つ．

4.2 確率的 EM アルゴリズム

Gibbs サンプリングアルゴリズムは，EM アルゴリズム2.2の変形とみなせる：

1. $E$-ステップ：EM アルゴリズムでは $$Q(\theta |\vartheta ,x):=(p_{\vartheta }dz|\log p_{\theta })$$ を評価するところであったが，Gibbs サンプリングでは，$p(z|x,\vartheta )$ のサンプリングを行う．
2. $M$-ステップ：EM アルゴリズムでは $$\underset{\vartheta \in \mathrm{\Theta }}{\mathrm{argmax}}Q(\vartheta |\theta ,x)=(p_{\theta }dz|\log p_{\vartheta })$$ を求めるところであったが，Gibbs サンプリングでは，$p(\theta |z,x)$ のサンプリングを行う．

これは $E$-ステップでの $Q$ 関数の評価が困難であるとき，$p(z|x,\theta )$ からのサンプリングでこれを回避できるという美点もある．

4.2.1 EM アルゴリズムへの部分的な適用：$E$-ステップ

またこの美点のみを用いて，$p(z|x,\theta )$ からサンプリングをして $Q$ の Monte Carlo 推定量 $$Q(\theta )=\frac{1}{M}\sum _{m=1}^M\log p(x,z^{(m)}|\theta )$$ を計算し，$M$-ステップとしてこれを最大化して $\{{\hat{\theta }}_j\}$ を得るという 確率的 EM アルゴリズム (Stochastic EM) も考えられる (Celeux and Diebolt, 1985)．²¹

この場合，$\{{\hat{\theta }}_j\}$ は多くの場合エルゴード的な Markov 連鎖を定めるが，これがどこに収束するかの特定が難しい (J. Diebolt and Ip, 1996)．

4.2.2 EM アルゴリズムへの部分的な適用：$M$-ステップ

Gibbs サンプリングの考え方を $M$-ステップにのみ導入し，$M$-ステップを完全に最大化するのではなく「条件付き最大化」に置き換えても，EM アルゴリズム本来の収束性は保たれる．$\theta =({\theta }_1,{\theta }_2)$ と分解できる際に，いずれか片方ずつのみを最大化する，などである．これを ECM (Expectation Conditional Maximization) アルゴリズムという (Meng and Rubin, 1991), (Meng and Rubin, 1993)．

$M$-ステップのみを確率的にすることで，EM アルゴリズムの局所解へのトラップを改善することができる．そのような手法の例に，SAME (State Augmentation for Marginal Estimation) (Doucet et al., 2002) などがある．

4.3 諸言

欠測モデル2.1のように，一般に グラフィカルモデル として知られる，局所的な関係のみから指定されるモデルや潜在変数を持つモデルでは，Gibbs サンプリングにより効率的に結合分布からサンプリングができる．

MCMC はグラフィカルモデルを用いた Bayes 推論の，強力な武器である．²²

References

Bishop, C. M. (2006). Pattern recognition and machine learning. Springer New York.

Boyles, R. A. (1983). On the convergence of the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 45(1), 47–50.

Carmer, H. (1946). Mathematical methods of statistics. Princeton University Press.

Celeux, G., and Diebolt, J. (1985). The SEM algorithm: A probabilistic teacher algorithm derived from the EM algorithm for the mixture problem. Computational Statistics Quarterly, 2, 73–82.

Chau, N. H., Moulines, É., Rásonyi, M., Sabanis, S., and Zhang, Y. (2021). On stochastic gradient langevin dynamics with dependent data streams: The fully nonconvex case. SIAM Journal on Mathematics of Data Science, 3(3), 959–986.

Chopin, N., and Papaspiliopoulos, O. (2020). An introduction to sequential monte carlo. Springer Cham.

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1), 1–22.

Diebolt, J., and Ip, E. (1996). Markov chain monte carlo in practice. In W. R. Gilks, S. Richardson, and D. Spiegelhalter, editors, pages 259–274. Chapman; Hall.

Diebolt, Jean, and Robert, C. P. (1994). Estimation of finite mixture distributions through bayesian sampling. Journal of the Royal Statistical Society. Series B (Methodological), 56(2), 363–375.

Doucet, A., Godsill, S. J., and Robert, C. P. (2002). Marginal maximum a posteriori estimation using markov chain monte carlo. Statistics and Computing, 12, 77–84.

Finch, S. J., Mendell, N. R., and Thode Jr., H. C. (1989). Probabilistic measures of adequacy of a numerical search for a global maximum. Journal of the American Statistical Association, 84(408), 1020–1023.

Fisher, R. A. (1912). On an absolute criterion for fitting frequency curves. Messenger of Mathematics, 41, 155–160.

Fletcher, R. (1987). Practical methods of optimization. John Wiley & Sons.

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The elements of statistical learning. Springer New York.

Hrafnkelsson, B. (Ed.). (2023). Statistical modeling using bayesian latent gaussian models: With applications in geophysics and encironmental sciences. Springer Cham.

Johnston, T., Makras, N., and Sabanis, S. (2024). Taming the interacting particle langevin algorithm – the superlinear case.

MacKay, D. J. C. (2003). Information theory, inference and learning algorithms. Cambridge University Press.

Meng, X.-L., and Rubin, D. B. (1991). Using EM to obtain asymptotic variance-covariance matrices: The SEM algorithm. Journal of the American Statistical Association, 86(416), 899–909.

Meng, X.-L., and Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80(2), 267–278.

Mengersen, K. L., and Robert, C. P. (1996). Testing for mixtures: A bayesian entropic approach. In Bayesian statistics 5: Proceedings of the fifth valencia international meetings, pages 255–276.

Neal, R. M., and Hinton, G. E. (1998). Learning in graphical models. In M. I. Jordan, editor, pages 355–368. Springer Dordrecht.

Robbins, H., and Monro, S. (1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22(3), 400–407.

Robert, C. P. (1996). Markov chain monte carlo in practice. In W. R. Gilks and D. Spiegelhalter, editors, pages 441–464. Chapman & Hall, London.

Robert, Christian P. (2007). The bayesian choice: From decision-theoretic foundations to computational implementation. Springer New York.

Robert, Christian P., and Casella, G. (2004). Monte carlo statistical methods. Springer New York.

Sun, Y., Babu, P., and Palomar, D. P. (2016). Majorization-minimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3), 794–816.

Tanner, M. A., and Wong, W. H. (1987). The calculation of posterior distributions by data augmentation. Journal of the American Statistical Association, 82(398).

Wainwright, M. J., and Jordan, M. I. (2008). Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1-2), 1–305.

Wei, G. C. G., and Tanner, M. A. (1990a). A monte carlo implementation of the EM algorithm and the poor man’s data augmentation algorithm. Journal of the American Statistical Association, 85(411), 699–704.

Wei, G. C. G., and Tanner, M. A. (1990b). Posterior computations for censored regression data. Journal of the American Statistical Association, 85(411), 829–839.

Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. The Annals of Statistics, 11(1), 95–103.

Wu, T. T., and Lange, K. (2010). The MM alternative to EM. Statistical Science, 25(4), 492–505.

甘利俊一. (1989). 神経回路網モデルとコネクショニズム,Vol. 22. 東京大学出版会.

鎌谷研吾. (2020). モンテカルロ統計計算. 講談社.

Footnotes

1. (Carmer, 1946, p. 498) によると，(Fisher, 1912) が初出であるが，以前に Gauss がその特別な形を用いていた．また，(Carmer, 1946, p. 499) での定義はこことは違い，尤度関数の停留点（＝尤度方程式の解）と定義している．

2. (Christian P. Robert and Casella, 2004, p. 174) 式(5.8)．

3. (Christian P. Robert and Casella, 2004, p. 174) 式(5.7)．$p$ を完備化された尤度 (completed likelihood) ともいう．$X$ を incomplete data, $(X,Z)$ を complete data ともいう (Bishop, 2006, pp. 433, p.440)．

4. 特に，隠れ Markov モデルの文脈では，EM アルゴリズムは Baum-Welch アルゴリズム とも呼ばれる (Chopin and Papaspiliopoulos, 2020, p. 70)．

5. それぞれ，(Christian P. Robert, 2007, p. 330)，(Christian P. Robert and Casella, 2004, p. 176)，(Hastie et al., 2009, p. 276)，

6. 正確には確率核 $Q:\mathcal{X}\times \mathrm{\Theta }\to \mathcal{Z}$．

7. この $F$ は多く $Q$ とも表され，$Q$-関数ともいう．$p(x|\theta )$ やその対数は 証拠 (evidence) ともいうので，$F$ は 証拠下界 (ELBO: Evidence Lower BOund) ともいう．

8. 式変形は (Bishop, 2006, p. 450) も参照．この $p(z|x,\theta )$ は観測 $x$ の下での，潜在変数 $z$ の条件付き分布である．しかし，このように双方を最大化ステップと見る変分法的な見方が出来るのである (Wainwright and Jordan, 2008, pp. 153–154), (Neal and Hinton, 1998), (Hastie et al., 2009, p. 277)．よって，この $E$-ステップも，GEM のように，必ずしも完全な最大化を達成する必要はないことがわかる (Neal and Hinton, 1998), (Bishop, 2006, p. 454)．例えば変分近似を行った場合，変分 EM アルゴリズムができあがる (Wainwright and Jordan, 2008, p. 154)．

9. (Johnston et al., 2024) にも言及あり．

10. (Christian P. Robert and Casella, 2004, p. 177) 定理5.15，(Christian P. Robert, 2007, p. 334) 演習6.52．

11. (Christian P. Robert and Casella, 2004, p. 178) 定理5.16．

12. この用語は (甘利俊一, 1989, p. 141) の 模擬除冷 の表現に触発された．

13. (Christian P. Robert and Casella, 2004) も参照．

14. (Bishop, 2006, pp. 436–439) も参照．

15. (MacKay, 2003, p. 303)．

16. (Bishop, 2006, p. 424)．

17. (Christian P. Robert and Casella, 2004, pp. 181–182) 例5.19 も参照．

18. (Christian P. Robert, 2007, p. 307) も参照．

19. $u_{-k}:=u_{1:(k-1),(k+1):K}$ とした．

20. (Christian P. Robert, 2007, p. 309) 補題6.3.6，(鎌谷研吾, 2020, p. 139) 定理5.7．

21. (Christian P. Robert and Casella, 2004, p. 200) 5.5.1 節も参照．

22. (Christian P. Robert, 2007, p. 318) も参照．