SAEM Versus Monte Carlo EM

Two Ways to Approximate the EM E-Step

In a nonlinear mixed-effects model, the observed data are available but the individual latent effects are not. The EM framework handles this by taking an expectation over latent effects conditional on the observed data and current parameter values. In many nonlinear models, that expectation cannot be evaluated exactly. Monte Carlo EM (MCEM) and stochastic approximation EM (SAEM) both use simulation, but they use it in different ways: MCEM repeatedly estimates an expectation for the current iteration, whereas SAEM recursively updates an estimate of the sufficient statistics needed by the M-step.

Let $y$ denote all observed data, $b$ the vector of latent random effects, and $\theta$ the population parameters. Write $\ell_c(\theta;y,b)$ for the complete-data log-likelihood, as if both $y$ and $b$ were observed. The conditional distribution $p(b\mid y,\theta_k)$ is the latent-effect distribution given the observed data at the current parameter value $\theta_k$.

Definition: The EM auxiliary function at iteration $k$ is $$ Q(\theta\mid\theta_k) = E_{\theta_k}\!\left[\ell_c(\theta;y,b)\mid y\right]. $$ MCEM approximates this conditional expectation anew at each iteration using a Monte Carlo sample. SAEM instead uses a stochastic-approximation recursion to update an estimate of the sufficient statistics that determine this expectation.

A common complete-data representation is

\ell_c(\theta;y,b) = c(y,b)-\psi(\theta)+\phi(\theta)^\mathsf{T}T(y,b),

where $T(y,b)$ is a vector of complete-data sufficient statistics, and $c$, $\psi$, and $\phi$ are functions that do not depend on the latent effects in the same way. In this case, the part of $Q$ that depends on $\theta$ is determined by $E_{\theta_k}[T(y,b)\mid y]$. This is why SAEM can update a statistic vector rather than repeatedly estimating the entire expected log-likelihood as a function of every candidate $\theta$.

What Each Algorithm Simulates

At an MCEM iteration, the algorithm obtains $M_k$ draws $b^{(k,1)},\ldots,b^{(k,M_k)}$ from the conditional distribution at $\theta_k$, or uses a simulation method whose draws target that distribution. It forms

\widehat Q_{M_k}(\theta\mid\theta_k) = \frac{1}{M_k} \sum_{m=1}^{M_k} \ell_c\!\left(\theta;y,b^{(k,m)}\right).

The approximate M-step then maximizes this estimated function, or a numerical approximation to it. The simulation and averaging are repeated at the next iteration using the updated parameter value. Independent conditional draws are one possible approach; correlated draws from a suitable Markov chain can also be used, provided the chain targets the required conditional distribution and Monte Carlo uncertainty is assessed appropriately.

SAEM uses a different recursion. Suppose the complete-data likelihood has the sufficient-statistic form above, and let $S_k$ estimate the conditional expectation of $T(y,b)$. After generating a latent-effect draw $b_k$ conditional on $y$ at the current parameter value, the basic update is

S_k=S_{k-1}+\gamma_k\{T(y,b_k)-S_{k-1}\}, \qquad 0<\gamma_k\leq 1.

The parameter update uses $S_k$ in place of the exact conditional expectation of $T$. The conditional distribution changes as the parameter estimates change, so the target expectation is not fixed throughout the run. The simulation mechanism must produce draws from the appropriate conditional distribution as the parameters evolve; in practice this may use a Markov chain whose transition kernel has that conditional distribution as its invariant distribution.

Key distinction: MCEM spends simulation effort estimating the current E-step with a fresh batch of draws. SAEM uses each conditional draw to make a gain-weighted correction to a running sufficient-statistic estimate. SAEM therefore does not require a large independent batch at every iteration, but its result depends on both the simulation chain and the gain sequence.

Worked Examples

Worked Example: One MCEM Average and One SAEM Update

Consider a hypothetical Gaussian random-intercept model with one observation $y=1$, latent effect $b\sim N(0,\omega^2)$, and $y\mid b\sim N(b,\sigma^2)$. For this illustration, $\omega^2=1$ and $\sigma^2=1$ are the current values. A sufficient statistic for the part of the complete-data log-likelihood involving $\omega^2$ is $T(b)=b^2$.

The conditional distribution is normal with variance $(1/\omega^2+1/\sigma^2)^{-1}=0.5$ and mean $0.5y=0.5$. Thus its exact conditional second moment is $E[b^2\mid y]=\operatorname{Var}(b\mid y)+E[b\mid y]^2=0.5+0.5^2=0.75$.

Suppose four hypothetical MCEM draws from this conditional distribution are $0,1,1,2$. The estimated sufficient-statistic expectation is $$ \widehat E[b^2\mid y] = \frac{0^2+1^2+1^2+2^2}{4} = \frac{6}{4} = 1.50. $$ This is a Monte Carlo estimate, not the exact value $0.75$; with only four draws, sampling error can be substantial.

Now consider a SAEM update with previous statistic estimate $S_{k-1}=0.80$, one hypothetical conditional draw $b_k=1.2$, and gain $\gamma_k=0.25$. Then $$ S_k = 0.80+0.25(1.2^2-0.80) = 0.80+0.25(0.64) = 0.96. $$ The two results differ because MCEM averages its current batch, while SAEM adjusts a value carrying information from previous iterations. Neither number is guaranteed to equal the exact conditional expectation on a particular iteration.

Worked Example: Monte Carlo Error in an MCEM E-Step

Suppose, in a hypothetical MCEM E-step, four independent conditional draws produce statistic values $T=b^2$ of $0.25, 1, 2.25,$ and $4$. The Monte Carlo estimate is $$ \bar T=\frac{0.25+1+2.25+4}{4}=1.875. $$

The sample standard deviation of these four values is $$ s_T= \sqrt{\frac{(0.25-1.875)^2+(1-1.875)^2+(2.25-1.875)^2+(4-1.875)^2}{3}} = \sqrt{2.6875} \approx 1.639. $$ The estimated Monte Carlo standard error of the average is therefore $$ SE(\bar T)=\frac{s_T}{\sqrt{4}}=\frac{1.639}{2}\approx 0.8195. $$

This standard error is large relative to the estimate $1.875$, so this small batch gives a noisy E-step estimate. If 16 independent draws had a similar standard deviation, the corresponding approximate standard error would be $1.639/\sqrt{16}=0.4098$, half as large. This is a planning comparison, not a guarantee that a particular batch of 16 draws will have exactly that standard error.

Worked Example: A Decreasing SAEM Gain

For a hypothetical scalar sufficient statistic, use the gain schedule $\gamma_k=k^{-0.7}$. At $k=10$, $\gamma_{10}=10^{-0.7}\approx0.1995$. If $S_9=0.900$ and the newly simulated statistic is $T(y,b_{10})=1.50$, then $$ S_{10}=0.900+0.1995(1.50-0.900) =0.900+0.1197 \approx1.020. $$

At $k=20$, $\gamma_{20}=20^{-0.7}\approx0.1228$. If the next simulated statistic is $0.600$, then $$ S_{20}=1.020+0.1228(0.600-1.020) =1.020-0.0516 \approx0.9684. $$ The later draw still changes the running estimate, but it receives less weight than a draw at iteration 10 would have received. For this schedule, $\sum_k\gamma_k$ diverges because $0.7\leq1$, while $\sum_k\gamma_k^2$ converges because $1.4>1$. These are standard stochastic-approximation conditions that allow continued learning while limiting accumulated noise under suitable regularity assumptions.

Why the Convergence Behavior Differs

For independent draws from the correct conditional distribution, an MCEM average of a statistic $T$ has Monte Carlo variance $\operatorname{Var}(T\mid y,\theta_k)/M_k$. Increasing $M_k$ reduces this variance. If a fixed, small number of draws is used at every iteration, however, the estimated E-step retains nonzero simulation noise. The parameter sequence may fluctuate around a solution rather than settle exactly. Increasing the number of draws as the fit proceeds, or otherwise ensuring Monte Carlo error becomes sufficiently small, can support convergence to a stationary point under suitable regularity conditions. There is no universal schedule that guarantees this for every model and sampling method.

SAEM manages simulation noise through the gain sequence. A commonly used condition is $\sum_k\gamma_k=\infty$ and $\sum_k\gamma_k^2<\infty$, together with regularity, stability, and adequate conditional simulation. The first condition prevents the algorithm from ceasing to respond to new information; the second limits the accumulated effect of random updates. Under appropriate conditions, SAEM can converge to a stationary point of the observed-data likelihood. This is a qualified result, not a promise that every run reaches the global maximum or that every finite update improves the observed likelihood.

The gain does not compensate for a simulation chain that fails to explore the conditional distribution. If successive draws are highly correlated, the information in them is less than in an equal number of independent draws. Poor mixing can impair MCEM estimates and SAEM updates alike. MCEM can often respond by taking more effective draws per E-step, while SAEM relies on a chain that explores adequately across iterations as well as on a suitable gain schedule.

Practical Comparison, Diagnostics, and Limitations

Feature	MCEM	SAEM
What is updated from simulation?	A Monte Carlo estimate of the current E-step objective or its sufficient-statistic expectations	A running estimate of sufficient statistics
Simulation pattern	A batch of draws is used at each parameter iteration	One or a modest number of conditional draws can update the recursion at each iteration
Main control on Monte Carlo noise	Number and effective independence of draws per E-step	Gain sequence, conditional-chain behavior, and accumulated stochastic updates
Typical finite-run behavior	Noisy E-step estimates can make parameter updates fluctuate	Individual updates are stochastic; decreasing gains reduce the influence of later noise

Diagnostics should match the source of randomness. For MCEM, examine whether the E-step estimates are stable enough for the parameter update: compare estimates across batches or assess their Monte Carlo standard errors, accounting for autocorrelation when draws come from a chain. If estimates change materially when the simulation effort increases, the E-step is not yet precise enough for confident interpretation. For SAEM, inspect whether the conditional simulation explores plausible latent-effect values and whether the sufficient-statistic updates settle as the gain decreases. Persistent drift can indicate parameter instability, poor mixing, or an unsuitable gain schedule rather than a need to run more iterations blindly.

Neither algorithm guarantees a global optimum in a nonlinear mixed-effects model. Results can depend on initialization, the shape of the likelihood, and how well the conditional latent-effect distribution is explored. Both methods also require computational effort to simulate latent effects; SAEM's smaller per-iteration simulation requirement does not automatically mean lower total cost, especially if its chain mixes slowly or many iterations are needed. Conversely, MCEM's larger batch can be wasteful early in a run when parameter estimates are far from stable.

Common mistake: “SAEM uses only one draw, so its E-step must be inaccurate; MCEM uses a batch, so it must be more reliable.” The first claim overlooks SAEM's recursive averaging across iterations, and the second overlooks MCEM's finite-sample error. The relevant questions are whether the conditional simulations explore adequately, whether MCEM's batch is large enough for its E-step, and whether SAEM's gain schedule is appropriate for its evolving target.

The choice is therefore not simply “many draws” versus “few draws.” MCEM is attractive when a sufficiently accurate batch expectation can be computed at each iteration and the batch size can be managed effectively. SAEM is attractive when the sufficient-statistic representation is available and recursive updates are computationally useful. Their shared dependence on conditional simulation means that neither should be trusted solely because its parameter estimates appear numerically stable.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

Question: What does MCEM average at iteration $k$?
Solution: It averages complete-data log-likelihood contributions, or the sufficient statistics that determine them, across draws targeting $p(b\mid y,\theta_k)$. The result estimates the current conditional expectation.
Question: If independent MCEM draws have statistic variance $9$, what is the Monte Carlo standard error of their average for $M=25$?
Solution: The standard error is $\sqrt{9/25}=3/5=0.600$, in the units of the statistic.
Question: With $S_{k-1}=2$, $T(y,b_k)=5$, and $\gamma_k=0.2$, calculate the SAEM update.
Solution: $S_k=2+0.2(5-2)=2.6$.
Question: Why can a fixed, small MCEM batch produce persistent parameter fluctuation?
Solution: Its conditional-expectation estimate retains Monte Carlo error at each iteration. That continuing noise can carry into the parameter updates rather than disappearing as the run proceeds.
Question: Does a decreasing SAEM gain ensure convergence if the conditional simulation chain does not explore the target distribution?
Solution: No. Gain conditions alone are insufficient; the latent-effect simulation must adequately target and explore the relevant conditional distribution, and other regularity and stability conditions are needed.

What You'll Learn