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.

Properties of this class derive their results from information in this dictionary. See docstring of fit_model_spectrum() for details.

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 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 or the algorithm.

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, ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]] | float]

Returns:

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

Notes

The returned dictionary contains the following items if the fit succeeds:

  • acint: estimate of the autocorrelation integral

  • acint_var: variance of the estimate of the autocorrelation integral

  • acint_std: standard error of the estimate 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

  • switch_exponent: exponent used to construct the cutoff

  • criterion: value of the criterion whose minimizer determines the frequency cutoff

  • criterion_expected: expected value of the criterion

  • criterion_var: expected variance of the criterion

  • criterion_zscore: z-score of the criterion

  • ll: log likelihood

  • pars_init: initial guess of the parameters

  • pars: optimized parameters

  • pars_covar: covariance matrix of the parameters

If the model can derive the exponential correlation time, The following properties are also included:

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

  • corrtime_exp_var: variance of the estimated exponential correlation time

  • corrtime_exp_std: standard error of the estimated 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

The ExpPolyModel has the following additional properties:

  • 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

If the fit fails, the following properties are set:

  • criterion: infinity

  • msg: error message

summarize_results(res, uc=None)[source]

Return a string summarizing the Result object.

Parameters: