stacie.estimate module

Algorithm to estimate the autocorrelation integral.

class Result(spectrum, model, cutoff_criterion, props, history)[source]

Bases: object

Container class holding all the results of the autocorrelation integral estimate.

Parameters:
property acint: float

The autocorrelation integral.

property acint_std: float

The uncertainty of the autocorrelation integral.

property corrtime_exp: float

The exponential correlation time.

property corrtime_exp_std: float

The uncertainty of the exponential correlation time.

property corrtime_int: float

The integrated correlation time.

property corrtime_int_std: float

The uncertainty of the integrated correlation time.

cutoff_criterion: CutoffCriterion

The criterion used to select or weight cutoff frequencies.

property fcut: int

The weighted average of the cutoff frequency.

history: list[dict[str]]

History of the cutoff optimization.

Each item is a dictionary returned by fit_model_spectrum(), containing the intermediate results of the fitting process. They are sorted from low to high cutoff frequency.

model: SpectrumModel

The model used to fit the low-frequency part of the spectrum.

property ncut: int

The number of points where the fitting weight is larger than 1/1000.

property neff: int

The weighted average of the effective number of frequencies used in the fit.

props: dict[str]

The properties marginalized over the ensemble of cutoff frequencies.

The following properties documented in fit_model_spectrum() are estimated as weighted averages over the cutoff frequencies:

  • amplitudes_model: model amplitudes at the included frequencies

  • acint: autocorrelation integral

  • acint_std: uncertainty of the autocorrelation integral

  • acint_var: variance of the autocorrelation integral

  • cost_zscore: z-score of the cost function

  • criterion_zscore: z-score of the cutoff criterion

  • fcut: cutoff frequency

  • pars: model parameters

  • pars_covar: covariance matrix of the parameters

Some properties are not averaged over cutoff frequencies:

  • ncut: number of points included in the fit, i.e. with weight larger than WEIGHT_EPS

  • switch_exponent: exponent used to construct the cutoff

  • weights: the weights used to combine the spectrum points in the fit

When using stacie.model.LorentzModel, the following properties are added (derived from the marginalized parameters and their covariance):

  • pars_lorentz: Lorentz parameters (converted from the Padé parameters)

  • pars_lorentz_covar: covariance matrix of the Lorentz parameters

  • corrtime_exp: exponential correlation time

  • corrtime_exp_var: variance of the exponential correlation time

  • corrtime_exp_std: standard error of the exponential correlation time

  • exp_simulation_time: recommended simulation time based on the exponential correlation time

  • exp_block_time: recommended block time based on the exponential correlation time

When using stacie.model.ExpPolyModel, the following additional properties are added (derived from the marginalized parameters and their covariance):

  • log_acint: the logarithm of the autocorrelation integral

  • log_acint_var: variance of the logarithm of the autocorrelation integral

  • log_acint_std: standard error of the logarithm of the autocorrelation integral

spectrum: Spectrum

The input spectrum from which the autocorrelation integral is estimated.

estimate_acint(spectrum, model, *, neff_min=None, neff_max=1000, fcut_min=None, fcut_max=None, fcut_spacing=0.5, switch_exponent=8.0, cutoff_criterion=None, rng=None, nonlinear_budget=100, criterion_high=100, verbose=False, uc=None)[source]

Estimate the integral of the autocorrelation function.

It is recommended to leave the keyword arguments to their default values, except for methodological testing.

This function fits a model to the low-frequency portion of the spectrum and derives an estimate of the autocorrelation (and its uncertainty) from the fit. It repeats this for a range of cutoff frequencies on a logarithmic grid. Finally, an ensemble average over all cutoffs is computed, by using -np.log of the cutoff criterion as weight.

The loop over all cutoff frequencies is performed in scan_frequencies(), while the marginalization over cutoff frequencies is done in marginalize_properties(). The function fit_model_spectrum() performs the actual fitting of the model for a given cutoff frequency.

The cutoff frequency grid is logarithmically spaced, with the ratio between two successive cutoff frequencies given by

\[\frac{f_{k+1}}{f_{k}} = \exp(g_\text{sp} / \beta)\]

where \(g_\text{sp}\) is fcut_spacing and \(\beta\) is switch_exponent.

Parameters:

  • spectrum (Spectrum) – The power spectrum and related metadata, used as inputs for the estimation of the autocorrelation integral. This object can be prepared with the function: stacie.spectrum.compute_spectrum().

  • model (SpectrumModel) – The model used to fit the low-frequency part of the spectrum.

  • neff_min (int | None) – The minimum effective number of frequency data points to include in the fit. (The effective number of points is the sum of weights in the smooth cutoff.) If not provided, this is set to 5 times the number of model parameters as a default.

  • neff_max (int | None) – The maximum number of effective points to include in the fit. This parameter limits the total computational cost. Set to None to disable this stopping criterion.

  • fcut_min (float | None) – The minimum cutoff frequency to use. If given, this parameter can only increase the minimal cutoff derived from neff_min.

  • fcut_max (float | None) – If given, cutoffs beyond this maximum are not considered.

  • fcut_spacing (float) – Dimensionless parameter that controls the spacing between cutoffs in the grid.

  • switch_exponent (float) – Controls the sharpness of the cutoff. Lower values lead to a smoother cutoff, and require fewer cutoff grid points. Higher values sharpen the cutoff, reveal more details, but a finer cutoff grid.

  • cutoff_criterion (CutoffCriterion | None) – The criterion function that is minimized to find the best cutoff frequency and, consequently, the optimal number of points included in the fit. If not given, the default is an instance of stacie.cutoff.CV2LCriterion.

  • rng (Generator | None) – A random number generator for sampling guesses of the nonlinear parameters. If not provided, np.random.default_rng(42) is used. The seed is fixed by default for reproducibility.

  • nonlinear_budget (int) – The number of samples used for the nonlinear parameters, calculated as nonlinear_budget ** num_nonlinear.

  • criterion_high (float) – An high increase in the cutoff criterion value, used to terminate the search for the cutoff frequency.

  • verbose (bool) – Set this to True to print progress information of the frequency cutoff search to the standard output.

  • uc (UnitConfig | None) – Unit configuration object used to format the screen output. If not provided, the default unit configuration is used. See stacie.utils.UnitConfig for details. This only affects the screen output (if any) and not the results!

Return type:

Result

Returns:

result – The inputs, intermediate results and outputs of the algorithm.

finalize_properties(props, model)[source]

Add remaining properties in-place.

Parameters:

  • props (dict[str]) – The properties dictionary to finalize. This is either the output of fit_model_spectrum() or the marginalized properties obtained by marginalize_properties(). This dictionary is modified in-place to add model-specific properties and to compute standard errors from variances.

  • model (SpectrumModel) – The model used to fit the spectrum.

fit_model_spectrum(spectrum, model, fcut, switch_exponent, cutoff_criterion, rng, nonlinear_budget)[source]

Optimize the parameter of a model for a given spectrum and cutoff frequency.

Parameters:

  • spectrum (Spectrum) – The spectrum object containing the input data.

  • model (SpectrumModel) – The model to be fitted to the spectrum.

  • fcut (float) – The cutoff frequency (in frequency units) used to construct the weights.

  • switch_exponent (float) – Controls the sharpness of the cutoff. Lower values lead to a smoother cutoff. Higher values sharpen the cutoff.

  • cutoff_criterion (CutoffCriterion) – The criterion function that is minimized to find the optimal cutoff (and thus determine the number of points to include in the fit).

  • rng (Generator) – A random number generator for sampling guesses of the nonlinear parameters.

  • nonlinear_budget (int) – The number of samples to use for the nonlinear parameters is nonlinear_budget ** num_nonlinear

Return type:

dict[str]

Returns:

props – A dictionary containing various intermediate results of the cost function calculation. See Notes for details.

Notes

The returned dictionary contains at least the following items, irrespective of whether the fit succeeds or fails:

  • fcut: cutoff frequency used

  • ncut: number of points included in the fit, i.e. with weight larger than WEIGHT_EPS

  • switch_exponent: exponent used to construct the cutoff

  • neff: effective number of points used in the fit (sum of weights)

  • pars_init: initial guess of the parameters

  • criterion: value of the cutoff criterion, or infinity if the fit fails.

  • msg: error message, if the fit fails

If the fit succeeds, the following additional statistical estimates are also set:

  • acint: autocorrelation integral

  • acint_var: variance of the autocorrelation integral

  • acint_std: standard error of the autocorrelation integral

  • cost_value: cost function value

  • cost_grad: cost gradient vector (if deriv >= 1)

  • cost_hess: cost Hessian matrix (if deriv == 2)

  • cost_hess_scales: Hessian rescaling vector, see robust_posinv.

  • cost_hess_rescaled_evals: Rescaled Hessian eigenvalues

  • cost_hess_rescaled_evecs: Rescaled Hessian eigenvectors

  • cost_expected: expected value of the cost function

  • cost_var: expected variance of the cost function

  • cost_zscore: z-score of the cost function

  • criterion_expected: expected value of the cutoff criterion

  • criterion_var: expected variance of the cutoff criterion

  • criterion_zscore: z-score of the cutoff criterion

  • ll: log likelihood

  • pars: model parameters

  • pars_covar: covariance matrix of the model parameters

When using stacie.model.LorentzModel, the following estimates are added:

  • pars_lorentz: Lorentz parameters (converted from the Padé parameters)

  • pars_lorentz_covar: covariance matrix of the Lorentz parameters

  • corrtime_exp: exponential correlation time, the slowest time scale in the sequences

  • corrtime_exp_var: variance of the exponential correlation time

  • corrtime_exp_std: standard error of the exponential correlation time

  • exp_simulation_time: recommended simulation time based on the exponential correlation time

  • exp_block_time: recommended block time based on the exponential correlation time

When using stacie.model.ExpPolyModel, the following estimates are added:

  • log_acint: the logarithm of the autocorrelation integral

  • log_acint_var: variance of the logarithm of the autocorrelation integral

  • log_acint_std: standard error of the logarithm of the autocorrelation integral

marginalize_properties(spectrum, model, history, *, switch_exponent=8.0, cutoff_criterion=None)[source]

Marginalize the properties over the ensemble of cutoff frequencies.

:param See estimate_acint() for parameter descriptions other than history.: :param The history parameter is the list of dictionaries returned by scan_frequencies().:

Return type:

Result

Returns:

result – The inputs, intermediate results and outputs of the algorithm. This object is returned by the function estimate_acint().

Parameters:
scan_frequencies(spectrum, model, *, neff_min=None, neff_max=1000, fcut_min=None, fcut_max=None, fcut_spacing=0.5, switch_exponent=8.0, cutoff_criterion=None, rng=None, nonlinear_budget=100, criterion_high=100, verbose=False, uc=None)[source]

Scan over cutoff frequencies and fit a model for each cutoff.

:param See estimate_acint() for parameter descriptions.:

Return type:

list[dict[str]]

Returns:

history – A list of dictionaries, one for each cutoff frequency, each containing various intermediate results of the fitting.

Parameters:
summarize_results(res, uc=None)[source]

Return a string summarizing the Result object.

Parameters: