SAEM Burn-In, Stochastic Steps, and Averaging

Why SAEM Settings Need to Be Considered Together

In stochastic approximation expectation-maximization (SAEM), simulated latent variables provide a noisy estimate of quantities needed to update model parameters. The simulation count, step-size schedule, and burn-in therefore affect how much Monte Carlo variation enters the updates and how quickly early behavior influences later estimates. Optional averaging of parameter estimates can make a reported sequence look smoother, but it is a separate operation from the SAEM stochastic update.

These controls have different roles. A larger simulation count can reduce the simulation noise in an iteration, a smaller step size gives that iteration less influence on the running stochastic approximation, and a longer burn-in allows more updates before an analyst summarizes the run. None, by itself, guarantees that the model is identified, the latent-variable simulation is adequate, or the algorithm has reached a stable region. The settings should be assessed together with traces of parameter estimates and the quantities being updated.

This update gives a useful interpretation of the step size: \(\gamma_k\) is the weight assigned to the current Monte Carlo estimate, while \(1-\gamma_k\) is the weight retained from the previous running statistic. A large step size adapts rapidly but lets a noisy estimate have more immediate influence. A small step size damps iteration-to-iteration variation, but also makes the approximation slower to respond when it is still far from a stable region. The parameter update is obtained from the current stochastic approximation according to the model’s maximization step; it is not generally identical to the running statistic itself.

Step-Size Schedules and Their Consequences

A common diminishing schedule, after any chosen initial phase, is

where \(c\) controls the early scale and \(a\) controls the rate of decline. Under standard stochastic-approximation regularity conditions, useful diminishing schedules typically have \(\sum_k \gamma_k=\infty\) and \(\sum_k\gamma_k^2<\infty\). For the displayed power schedule, the first sum diverges when \(a\leq1\), and the second converges when \(2a>1\). These conditions are part of convergence theory, not a guarantee that a finite run is adequate for a particular nonlinear mixed-effects model.

Some implementations instead use an initial phase with \(\gamma_k=1\), followed by a decreasing schedule. With \(\gamma_k=1\), the running statistic is replaced by that iteration’s Monte Carlo estimate, so early estimates can move substantially. Such a phase can help adaptation, but its duration and the subsequent schedule need to be documented. A schedule that stays fixed at a positive step size does not have the same asymptotic averaging behavior as one whose step size diminishes.

Simulation Count: Less Noise at More Computational Cost

The simulation count \(M\) is the number of Monte Carlo draws used to form an iteration’s estimate \(\widehat{s}_k\). If those draws are independent and have variance \(\tau^2\) for a particular statistic, their sample mean has variance \(\tau^2/M\) and standard error \(\tau/\sqrt{M}\). In practice, latent-variable draws may be correlated, especially when generated by a Markov chain. In that case, \(M\) draws do not provide the precision of \(M\) independent draws; autocorrelation reduces the effective information in the simulation.

Increasing \(M\) and decreasing \(\gamma_k\) are not interchangeable. More draws can make the current \(\widehat{s}_k\) less noisy, while a smaller step size reduces the weight of that noisy estimate in the update. More draws can require substantially more computation, and a very small step size can delay movement toward the relevant region. A practical choice depends on whether the observed instability appears to arise from noisy iteration-level estimates, an overly aggressive schedule, poor mixing of the latent-variable simulation, or another model or optimization problem.

Choosing and Evaluating Burn-In

Burn-in is a practical boundary for deciding which iterations to use in a post-run summary or optional average. Early parameter estimates may move sharply as the stochastic approximation responds to its starting values and simulation noise. A burn-in that is too short can include this transient behavior in the summary; one that is unnecessarily long discards iterations and does not repair a problematic run. There is no universal burn-in length that applies across models or datasets.

Inspect traces over the full run, including parameter estimates and important stochastic quantities, rather than choosing the burn-in solely because the last few values look stable. Compare runs with different random seeds when feasible, and check whether estimates after the proposed boundary are broadly consistent. Large persistent trends, abrupt jumps, repeated boundary behavior, or pronounced disagreement across runs are reasons to investigate rather than simply extend the excluded period. A flat trace is also not proof of correctness: a poorly explored or weakly identified model may appear stable.

The term “burn-in” can also refer to early draws discarded from a Markov chain used to simulate latent variables. That is a different decision from excluding early SAEM iterations from an estimate summary. The chain’s mixing and the algorithm’s stochastic-approximation transient should be assessed separately; discarding early SAEM parameter values does not establish that the latent-variable chain mixed well.

Parameter Averaging Is Not the SAEM Update

In SAEM, the stochastic approximation already averages simulated sufficient statistics recursively. Optional parameter averaging instead forms a summary from parameter estimates \(\theta_k\), often over iterations after a selected burn-in \(B\):

This average can reduce visible iteration-to-iteration fluctuation when estimates vary around a stable region. It does not create new information, replace the model’s maximization step, or by itself establish convergence or quantify sampling uncertainty. Because successive iterations are dependent, the variability across the sequence is not automatically a standard error for \(\bar{\theta}_{B:K}\).

Averaging also needs to respect the parameter’s meaning and constraints. For a positive parameter, arithmetic averaging of valid positive estimates remains positive, but an average across a long transient can still be misleading. For a correlation constrained to \([-1,1]\), the arithmetic mean remains in range if every input is in range, but its interpretation and behavior near boundaries still merit inspection. If an algorithm operates on transformed parameters, the average on that transformed scale need not equal the average on the natural scale. State which scale and which iterations were averaged.

The main alternatives serve different purposes. The final iterate is simple and corresponds directly to the last update, but can reflect late stochastic fluctuation. A post-burn-in average can smooth the reported trajectory, but can conceal drift if the run has not stabilized. Increasing simulation count can improve the quality of each Monte Carlo estimate, at additional cost. Changing the step-size schedule changes the balance between adaptation and noise across iterations. These are complementary controls; no single option is a substitute for checking the others.

Diagnostics, Failure Modes, and a Practical Workflow

A sensible settings review begins by recording the schedule, simulation count, burn-in boundary, and any parameter-averaging rule. Then inspect parameter traces and the underlying stochastic quantities. Look for continued directional movement after the proposed burn-in, large jumps that persist despite a decreasing step size, unstable estimates that shrink when simulation effort increases, and disagreement among independent runs. If latent-variable simulation is correlated, assess that simulation’s mixing rather than assuming the nominal draw count reflects independent information.

These diagnostics help distinguish remedies. If per-iteration Monte Carlo noise is prominent, increasing \(M\) may help. If updates are erratic because each current estimate has too much influence, a more gradual step-size schedule may help, though it slows adaptation. If estimates still trend after the proposed burn-in, simply averaging them can hide the trend; investigate initialization, simulation behavior, or model fit before selecting a later boundary. If runs differ substantially, a polished average from one run is not evidence that the discrepancy has been resolved.

Common Pitfall

A common mistake is to say, “The parameter trace is noisy, so I will average it and the estimate will be more precise.” Averaging can reduce visible fluctuations in a sequence, but it does not prove the sequence is centered on the desired estimate, and dependent iterates are not independent replications. The correction is to diagnose why the trace is noisy, assess the simulation and schedule, and report the averaging rule only as a deliberate summary choice. Likewise, discarding more iterations does not correct poor latent-variable mixing or a schedule that has not adapted adequately.