Exercises

If $\varphi$ is concave and $X$ is a random variable with $\mathbb{E}|X| < \infty$, then $\varphi(\mathbb{E}[X]) \geq \mathbb{E}[\varphi(X)]$.

The logarithm is concave on $(0,\infty)$. Hence $\log(\mathbb{E}[X]) \geq \mathbb{E}[\log X]$ for positive $X$. Equality holds iff $X$ is almost surely constant. $\blacksquare$

$\mathcal{F}(q,\boldsymbol{\theta}) = \int q(\mathbf{z})\log\frac{p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})}{q(\mathbf{z})}\,d\mathbf{z} = \int q(\mathbf{z})\log\frac{p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})\,f(\mathbf{y};\boldsymbol{\theta})}{q(\mathbf{z})}\,d\mathbf{z}$.

$= \log f(\mathbf{y};\boldsymbol{\theta}) - \int q(\mathbf{z})\log\frac{q(\mathbf{z})}{p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})}\,d\mathbf{z} = \log f(\mathbf{y};\boldsymbol{\theta}) - D(q\Vert p(\cdot\mid\mathbf{y};\boldsymbol{\theta}))$. Rearranging gives the stated identity. $\blacksquare$

$\sum_k \gamma_{ik}^{(t)} = \frac{\sum_k \pi_k^{(t)}\,\mathcal{N}(\mathbf{y}_i;\boldsymbol{\mu}_k^{(t)},\boldsymbol{\Sigma}_{k}^{(t)})} {\sum_j \pi_j^{(t)}\,\mathcal{N}(\mathbf{y}_i;\boldsymbol{\mu}_j^{(t)},\boldsymbol{\Sigma}_{j}^{(t)})} = 1$. $\blacksquare$

$\mathcal{L}(\boldsymbol{\pi},\lambda) = \sum_k N_k \log\pi_k - \lambda\!\left(\sum_k \pi_k - 1\right)$.

$\partial \mathcal{L}/\partial \pi_k = N_k/\pi_k - \lambda = 0$ $\Rightarrow$ $\pi_k = N_k/\lambda$.

$\sum_k \pi_k = 1 \Rightarrow \lambda = \sum_k N_k = n$. Hence $\pi_k^{(t+1)} = N_k/n$. Non-negativity holds automatically since $N_k \geq 0$. $\blacksquare$

$\gamma_{i1}^{(t)} = \frac{\exp(-(y_i-\mu_1^{(t)})^2/2)}{\exp(-(y_i-\mu_1^{(t)})^2/2)+\exp(-(y_i-\mu_2^{(t)})^2/2)}$ $= \sigma\!\big((\mu_1^{(t)}-\mu_2^{(t)})\,y_i - \tfrac12((\mu_1^{(t)})^2-(\mu_2^{(t)})^2)\big)$ (logistic sigmoid); $\gamma_{i2}^{(t)} = 1 - \gamma_{i1}^{(t)}$.

$\mu_k^{(t+1)} = \sum_i \gamma_{ik}^{(t)} y_i / \sum_i \gamma_{ik}^{(t)}$ for $k=1,2$.

Swapping $(\mu_1,\mu_2)$ swaps $(\gamma_{i1},\gamma_{i2})$ for every $i$, hence swaps the M-step outputs. The symmetric fixed point $\mu_1 = \mu_2$ is stationary — breaking symmetry in initialization is essential. $\blacksquare$

Assume dominated convergence applies (bound on $|\nabla_{\boldsymbol{\theta}} p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})|$ uniform in $\boldsymbol{\theta}$ locally). Then $\nabla f(\mathbf{y};\boldsymbol{\theta}) = \int \nabla p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})\,d\mathbf{z}$.

$\nabla \log f(\mathbf{y};\boldsymbol{\theta}) = \frac{\nabla f(\mathbf{y};\boldsymbol{\theta})}{f(\mathbf{y};\boldsymbol{\theta})} = \int \frac{p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})}{f(\mathbf{y};\boldsymbol{\theta})} \nabla \log p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})\,d\mathbf{z}$.

The fraction equals $p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})$, so the right-hand side is the posterior expectation of the complete score, as claimed. $\blacksquare$

Let $s(\boldsymbol{\theta},\mathbf{z}) = \nabla \log p(\mathbf{y},\mathbf{z};\boldsymbol{\theta})$. Then $\nabla \ell(\boldsymbol{\theta}) = \int p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})s(\boldsymbol{\theta},\mathbf{z})\,d\mathbf{z}$. Differentiating once more using the product rule: $\nabla^2 \ell = \int \nabla p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})\,s^\top\,d\mathbf{z} + \int p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})\,\nabla s\,d\mathbf{z}$.

$\nabla p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta}) = p(\mathbf{z}\mid\mathbf{y};\boldsymbol{\theta})\, (s(\boldsymbol{\theta},\mathbf{z}) - \nabla \ell(\boldsymbol{\theta}))$.

Plugging in and simplifying yields Louis's formula. The first term is the expected complete-data Hessian; the second is (minus) the posterior covariance of the complete-data score. $\blacksquare$

For each censored $i$ (i.e. $y_i^+ = 0$), the latent variable $y_i$ satisfies $y_i \leq 0$ and $y_i \sim \mathcal{N}(\theta,1)$.

For censored $i$: $\hat{y}_i^{(t)} = \mathbb{E}[y_i \mid y_i \leq 0, \theta^{(t)}] = \theta^{(t)} - \frac{\phi(-\theta^{(t)})}{\Phi(-\theta^{(t)})}$ (the mean of a Gaussian truncated to $(-\infty,0]$), where $\phi,\Phi$ are the standard normal pdf/cdf. For uncensored $i$: $\hat{y}_i = y_i^+$.

$\theta^{(t+1)} = \frac{1}{n}\sum_i \hat{y}_i^{(t)}$ — simple sample mean over the imputed values.

False. EM converges to a stationary point (local max, saddle point, or sometimes a boundary). The log-likelihood for mixtures and most latent-variable models is non-concave, so local maxima are generic. Multi-start is the standard remedy.

Fix $\pi_1 = 1/2$ (say), $\boldsymbol{\mu}_1 = \mathbf{y}_1$, $\boldsymbol{\Sigma}_{1} = \varepsilon \mathbf{I}$, and keep the other component with any finite parameters covering all data.

$\mathcal{N}(\mathbf{y}_1;\mathbf{y}_1,\varepsilon \mathbf{I}) = (2\pi\varepsilon)^{-d/2}\to\infty$ as $\varepsilon\to 0$. Hence $\log f(\mathbf{y}_1;\boldsymbol{\theta}) \geq \log(\pi_1 (2\pi\varepsilon)^{-d/2})\to\infty$, and the total log-likelihood (a sum over samples) also diverges. Remedy: regularize the covariances. $\blacksquare$

Holding $\{\boldsymbol{\mu}_k\}$ fixed, reassigning each $i$ to $\arg\min_k \|\mathbf{y}_i - \boldsymbol{\mu}_k\|^2$ minimizes its contribution to $J$, so $J$ does not increase.

Holding $\{c_i\}$ fixed, $J$ becomes a separable sum of quadratics in $\boldsymbol{\mu}_k$; each is minimized by the mean of its cluster. So the update does not increase $J$ either. Combining, $J$ is monotonically non-increasing. Since $J \geq 0$ and there are finitely many possible assignments, Lloyd's algorithm converges in finitely many iterations. $\blacksquare$

$\mathbb{E}_{\mathbf{x}\mid\mathbf{y};\boldsymbol{\alpha}^{(t)}}[\tfrac12 \log\alpha_j - \tfrac12 \alpha_j x_j^2] = \tfrac12 \log\alpha_j - \tfrac12 \alpha_j (\mu_j^{(t),2} + \boldsymbol{\Sigma}_{jj}^{(t)})$.

Differentiating with respect to $\alpha_j$: $1/(2\alpha_j) - \tfrac12(\mu_j^{(t),2} + \boldsymbol{\Sigma}_{jj}^{(t)}) = 0$ $\Rightarrow \alpha_j^{(t+1)} = 1/(\mu_j^{(t),2} + \boldsymbol{\Sigma}_{jj}^{(t)})$. $\blacksquare$

If $N_k^{(t)} = \sum_i \gamma_{ik}^{(t)} = 0$, the M-step divisions $\boldsymbol{\mu}_k^{(t+1)}, \boldsymbol{\Sigma}_{k}^{(t+1)}$ become $0/0$ and the component is dead — it cannot recover because zero responsibility gives zero prior in the next E-step. In practice, maintain $\gamma_{ik}^{(t)} \geq \epsilon$ by regularizing $\pi_k$, or drop the collapsed component and restart with $K-1$.

$\ell(\boldsymbol{\theta}^{(t)})$ is non-decreasing (Theorem 8.2) and bounded on $\Lambda_0$, hence converges to some $\ell^{\star}$.

Let $\boldsymbol{\theta}^{\star}$ be a limit point, $\boldsymbol{\theta}^{(t_k)}\to\boldsymbol{\theta}^{\star}$. By continuity of the EM map and uniqueness of the accumulation's log-likelihood $\ell^{\star}$, $\boldsymbol{\theta}^{\star}$ is a fixed point: $\boldsymbol{\theta}^{\star} \in \arg\max Q(\cdot\mid\boldsymbol{\theta}^{\star})$.

Interior fixed points satisfy $\nabla_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}^{\star}\mid\boldsymbol{\theta}^{\star}) = 0$. Fisher's identity gives $\nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}^{\star}) = 0$. $\blacksquare$

$\gamma_i^{(t)} = \Pr(Z_i = 1\mid y_i;\pi^{(t)})$ as above.

$\pi^{(t+1)} = \frac{1}{n}\sum_{i=1}^n \gamma_i^{(t)}$ — the average soft count of samples assigned to component 1.

For censored $j$ of observation $i$, let $\hat{y}_{ij}^{(t)} = \mathbb{E}[y_{ij}\mid |y_{ij}|>c, \theta_j^{(t)}] = \theta_j^{(t)} + \frac{\phi(c-\theta_j^{(t)}) - \phi(-c-\theta_j^{(t)})}{1 - \Phi(c-\theta_j^{(t)}) + \Phi(-c-\theta_j^{(t)})}$.

$\boldsymbol{\theta}^{(t+1)} = \frac{1}{n}\sum_i \hat{\mathbf{y}}_i^{(t)}$ (component-wise, using imputed or observed values).

$\max_{q_i}$ of $\mathcal{F}$ has unconstrained solution $q_i^{\star}(k) \propto \pi_k\,\mathcal{N}(\mathbf{y}_i;\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_{k}) = \gamma_{ik}$.

$\max_{\boldsymbol{\theta}}$ of $\mathcal{F}$ with $q = \prod_i q_i$ reduces to maximizing $Q(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)})$, which recovers the Theorem 8.4 updates. $\blacksquare$

Partially correct. The E-step sets $q^{(t+1)} = p(\cdot\mid\mathbf{y};\boldsymbol{\theta}^{(t)})$, so the KL drops to zero with respect to the current $\boldsymbol{\theta}^{(t)}$. The M-step then moves $\boldsymbol{\theta}$, which may increase the KL to the new posterior. The quantity that monotonically increases is $\ell(\boldsymbol{\theta})$, not the negative KL.

Write the EM map as $\boldsymbol{\theta}^{(t+1)} = M(\boldsymbol{\theta}^{(t)})$. At the fixed point $\boldsymbol{\theta}^{\star}$, the Jacobian is $DM(\boldsymbol{\theta}^{\star}) = \mathbf{I} - \mathbf{I}_c^{-1}(\boldsymbol{\theta}^{\star})\mathbf{I}_o(\boldsymbol{\theta}^{\star})$ (Dempster-Laird-Rubin 1977, §3.2).

The convergence rate is the spectral radius of $DM(\boldsymbol{\theta}^{\star})$. When little information is missing ($\mathbf{I}_o \approx \mathbf{I}_c$) the rate is small (fast convergence). When much is missing ($\mathbf{I}_o \ll \mathbf{I}_c$) the rate approaches 1 (slow convergence) — which motivates acceleration methods.

Log-likelihood on the training data increases monotonically with $K$ (more flexibility $\Rightarrow$ more overfitting). In the limit $K=n$, each component is a Dirac at a data point and the likelihood diverges.

Use a criterion that penalizes complexity: BIC $= -2\ell + p\log n$ (where $p$ is the number of free parameters), AIC $= -2\ell + 2p$, or cross-validated likelihood. For GMMs, BIC is standard.