Model Spectrum¶

STACIE supports three models for fitting the low-frequency part of the power spectrum. In all models, the value at zero frequency corresponds to the autocorrelation integral.

The ExpPolyModel is the most general; it is an exponential function of a linear combination of simple monomials of the frequency. You can specify the degrees of the monomials, and typically a low degree works fine:
- Degree \(0\) is suitable for a white noise spectrum.
- Degrees \(\{0, 1\}\) can be used to extract useful information from a noisy spectrum.
- Degrees \(\{0, 1, 2\}\) are applicable to spectra with low statistical uncertainty, e.g., averaged over more than 100 inputs.
- An even polynomial with degrees \(\{0, 2\}\) is suitable for spectra that are expected to have a vanishing derivative at zero frequency.
The main advantage of this model is its broad applicability, as it requires little prior knowledge of the functional form of the spectrum.
The PadeModel approximates the spectrum as a rational function, i.e., the ratio of two polynomials. Such rational functions can be parameterized to have well-behaved high-frequency tails, which facilitates the regression.
The LorentzModel is a special case of the Padé model with a Lorentzian peak at the origin plus some white noise. It is equivalent to the Padé model with numerator degrees \(\{0, 2\}\) and denominator degrees \(\{2\}\). This special case is not only implemented for convenience, since it is the most common way of using the Padé model, but also because it allows STACIE to derive the exponential correlation time.

1. ExpPoly Model¶

The ExpPolyModel is defined as:

\[ I^\text{exppoly}(f) = \exp\left(\sum_{s \in S} b_s f^s\right) \]

where \(S\) is the set of polynomial degrees, which must include 0. With this form, \(\exp(b_0)\) corresponds to the integral of the autocorrelation function. When one obtains an estimate \(\hat{b}_0\) and its variance \(\hat{\sigma}^2_{b_0}\), the autocorrelation integral is log-normally distributed with estimated mean and variance:

\[\begin{split} \begin{aligned} \hat{\mathcal{I}} &= \exp\left(\hat{b}_0 + \frac{1}{2}\hat{\sigma}^2_{b_0}\right) \\ \hat{\sigma}^2_{\mathcal{I}} &= \exp\left(2\hat{b}_0 + \hat{\sigma}^2_{b_0}\right) \left(\exp(\hat{\sigma}^2_{b_0}) - 1 \right) \end{aligned} \end{split}\]

To construct this model, you can create an instance of the ExpPolyModel class as follows:

from stacie import ExpPolyModel
model = ExpPolyModel([0, 1, 2])

This model is identified as exppoly(0, 1, 2) in STACIE’s screen output and plots.

2. Pade Model¶

The PadeModel is defined as:

\[ I^\text{pade}(f) = \frac{ \displaystyle \sum_{s \in S_\text{num}} p_s f^s }{ \displaystyle 1 + \sum_{s \in S_\text{den}} q_s f^s } \]

where \(S_\text{num}\) contains the polynomial degrees in the numerator, which must include 0, and \(S_\text{den}\) contains the polynomial degrees in the denominator, which must exclude 0. With this model, \(p_0\) corresponds to the integral of the autocorrelation function, for which we simply have:

\[\begin{split} \begin{aligned} \hat{\mathcal{I}} &= \hat{p}_0 \\ \hat{\sigma}^2_{\mathcal{I}} &= \hat{\sigma}^2_{p_0} \end{aligned} \end{split}\]

To construct this model, you can create an instance of the PadeModel class as follows:

from stacie import PadeModel
model = PadeModel([0, 2], [2])

This model is identified as pade(0, 2; 2) in STACIE’s screen output and plots.

3. Lorentz Model¶

The LorentzModel assumes that for large time lags, \(\Delta_t\), the autocorrelation function (ACF) exhibits an exponential decay superimposed with a white noise background:

\[ C(\Delta_t) = C_0\, \delta(\Delta_t) + C_1 \exp\left(-\frac{|\Delta_t|}{\tau_\text{exp}}\right) \]

where \(C_0\) is the white noise level, \(C_1\) is the amplitude of the exponential decay, and \(\tau_\text{exp}\) is the exponential correlation time. The Lorentz model is essentially the corresponding power spectral density:

\[ I^\text{lorentz}(f) = C_0 + \frac{2 C_1 \tau_\text{exp}}{1 + (2 \pi f \tau_\text{exp})^2} \]

This model can be expressed as a special case of the Padé model, with numerator degrees \(\{0, 2\}\) and denominator degrees \(\{2\}\), and is fitted accordingly. The Padé model corresponds to a Lorentzian peak only if \(q_2 > 0\) and \(p_0\, q_2 > p_2\). In this case, \(\tau_\text{exp}\) is related to the width of the peak (\(2 \pi \tau_\text{exp}\)) in the power spectrum. When these conditions are met after the regression, the parameters of the Padé model are converted as follows:

\[\begin{split} \begin{aligned} \hat{C}_0 &= \frac{\hat{p}_2}{\hat{q}_2} \\ \hat{C}_1 &= \frac{\pi}{\sqrt{\hat{q}_2}}\left(\hat{p}_0 - \frac{\hat{p}_2}{\hat{q}_2}\right) \\ \hat{\tau}_\text{exp} &= \frac{\sqrt{\hat{q}_2}}{2 \pi} \end{aligned} \end{split}\]

The covariance of these parameters is derived from the covariance of the Padé parameters using first-order error propagation. We first introduce the following Jacobian matrix:

\[\begin{split} \begin{aligned} \mathbf{J} &= \begin{pmatrix} \frac{\partial C_0}{\partial p_0} & \frac{\partial C_0}{\partial p_2} & \frac{\partial C_0}{\partial q_2} \\ \frac{\partial C_1}{\partial p_0} & \frac{\partial C_1}{\partial p_2} & \frac{\partial C_1}{\partial q_2} \\ \frac{\partial \tau_\text{exp}}{\partial p_0} & \frac{\partial \tau_\text{exp}}{\partial p_2} & \frac{\partial \tau_\text{exp}}{\partial q_2} \end{pmatrix} = \begin{pmatrix} 0 & \frac{1}{q_2} & -\frac{p_2}{q_2^2} \\ \frac{\pi}{\sqrt{q_2}} & -\frac{\pi}{q_2^{3/2}} & \frac{\pi}{2 q_2^{3/2}}\left(\frac{3p_2}{q_2} - p_0\right) \\ 0 & 0 & \frac{1}{4 \pi \sqrt{q_2}} \end{pmatrix} \end{aligned} \end{split}\]

and use it to compute the covariance of the Lorentz parameters:

\[ \hat{\mathbf{C}}_\text{lorentz} = \mathbf{J} \, \hat{\mathbf{C}}_\text{pade(0, 2; 2)} \, \mathbf{J}^T \]

Note that this model is also applicable to data whose short-time correlations are not exponential, as long as the tail of the ACF decays exponentially. Such deviating short-time correlations will only affect the white noise level \(\hat{C}_0\) and features in the PSD at higher frequencies, which are ignored by STACIE.

The implementation of the Lorentz model has the following advantages over the equivalent Padé model:

The exponential correlation time and its uncertainty are computed.
If no exponential correlation time can be computed, i.e. when \(\hat{q}_2 \le 0\) and \(\hat{p}_0\, \hat{q}_2 \le \hat{p}_2\), the fit is not retained for the final average over all cutoff grid points. (These are typically poor fits.)
After the optimization of the model parameters at a given frequency cutoff, the following two heuristics are applied to exclude or downweight fits with a high relative error of the exponential correlation time:
- If the relative error of the exponential correlation time is 100 times larger than that of the autocorrelation integral, the fit is discarded.
- In all other cases, a penalty is added to the cutoff criterion, which is proportional to the ratio of the relative error of the exponential correlation time and the autocorrelation integral.
This heuristic is based on the observation that large relative errors of the exponential correlation time often indicate poor fits for which the MAP estimate and the Laplace approximation of the posterior distribution tend to be unreliable. By taking a ratio of relative errors, the penalty is dimensionless and insensitive to the overall uncertainty of the spectrum.

The result of these two heuristics is that the final estimate of the exponential correlation time, after averaging over all frequency cutoffs, is more robust and reliable.

Note that the hyperparameters of the penalty can be altered, but it is recommended to keep their default values.

To construct this model, you can create an instance of the LorentzModel class as follows:

from stacie import LorentzModel
model = LorentzModel()

This model is identified as lorentz() in STACIE’s screen output and plots.