A spatiotemporal framework for estimating trial-to-trial variability in evoked response (ER) data is presented. Spatial and temporal bases capture the aspects of the response that are consistent across trials, while the basis expansion coefficients represent the variable components of the response. We focus on the simplest case of constant spatiotemporal response shape and varying amplitude across trials. Two different constraints on the amplitude evolution are employed to effectively integrate the individual responses and improve robustness at low SNR. The linear dynamical system response constraint estimates the current trial amplitude as an unknown constant scaling of the estimate in the previous trial plus zero-mean Gaussian noise with unknown variance. The independent response constraint estimates response amplitudes across trials as independent Gaussian random variables having unknown mean and variance. We develop a generalized expectation-maximization algorithm to obtain the maximum-likelihood (ML) estimates of the signal waveform, noise covariance matrix, and unknown constraint parameters. ML source localization is achieved by scanning the likelihood over different sets of spatial bases. We demonstrate the variability estimation and source localization effectiveness of the proposed algorithms using both real and simulated ER data.