  Time-dependent ensemble averages, i.e., trajectory-based averages of some
observable, are of importance in many fields of science. A crucial objective
when interpreting such data is to fit these averages (for instance, squared
displacements) with a function and extract parameters (such as diffusion
constants). A commonly overlooked challenge in such function fitting procedures
is that fluctuations around mean values, by construction, exhibit temporal
correlations. We show that the only available general purpose function fitting
methods, correlated chi-square method and the weighted least squares method
(which neglects correlation), fail at either robust parameter estimation or
accurate error estimation. We remedy this by deriving a new closed-form error
estimation formula for weighted least square fitting. The new formula uses the
full covariance matrix, i.e., rigorously includes temporal correlations, but is
free of the robustness issues, inherent to the correlated chi-square method. We
demonstrate its accuracy in four examples of importance in many fields:
Brownian motion, damped harmonic oscillation, fractional Brownian motion and
continuous time random walks. We also successfully apply our method, weighted
least squares including correlation in error estimation (WLS-ICE), to particle
tracking data. The WLS-ICE method is applicable to arbitrary fit functions, and
we provide a publically available WLS-ICE software.
