pyerrors: a python framework for error analysis of Monte Carlo data

We present the pyerrors python package for statistical error analysis of Monte Carlo data. Linear error propagation using automatic differentiation in an object oriented framework is combined with the $\Gamma$-method for a reliable estimation of autocorrelation times. Data from different sources can easily be combined, keeping the information on the origin of error components intact throughout the analysis. pyerrors can be smoothly integrated into the existing scientific python ecosystem which allows for efficient and compact analyses.


Introduction
The Standard Model of particle physics (SM) is very successful in describing the basic constituents of matter and their interactions but falls short to explain a number of observed phenomena such as the matter-antimatter asymmetry of the universe or the existence of dark matter.The search for physics beyond the Standard Model relies on highly precise theoretical calculations of SM observables.Lattice QCD is the only known ab-initio approach to compute these observables in the low energy, hadronic regime of quantum chromodynamics (QCD), describing the strong interactions in the SM.The method is based on a Markov Chain Monte Carlo sampling of the QCD action in Euclidean space, formulated via the path integral formalism.
In recent years, lattice QCD calculations have become a precision tool such that they have a relevant impact on phenomenology and the search for beyond the SM theories, see Reference [1] for a recent review of lattice results.Prime examples are the computation of the strong coupling constant of QCD and of hadronic contributions to the anomalous magnetic moment of the muon [2], where recent studies achieve precision similar to data-driven methods.The main contribution to the systematic uncertainties of these calculations and therefore the main obstacle for obtaining even higher precision is due to the extrapolation to the continuum limit of the theory.
When simulating boxes with ever-increasing resolution, thereby lifting the ultraviolet cutoff, state of the art algorithms experience critical slowing down.This results in a significant increase of the autocorrelations towards the continuum limit [3], which are inherent in Markov Chain Monte Carlo (MCMC) methods.Thus, the analysis of data stemming from MCMC methods requires a careful treatment of these autocorrelations.Otherwise statistical errors might be underestimated, putting the significance of results in question.
The Γ-method as introduced in [4,5] has been shown to decrease systematic errors in the treatment of autocorrelation in comparison to widely used techniques such as binned Jackknife [6,7] or Bootstrap resampling [8], especially in the presence of critical slowing down [3].Whereas the first publicly available implementations [3,5,9] do not allow to combine observables from simulations with different parameters, this is possible in the implementations of [10,11] written in Fortran and Julia. 1n this article, we present the pyerrors package that provides a framework to estimate and propagate uncertainties of Monte Carlo data in python.The guiding design principle of pyerrors is to provide a new data type, Obs, which integrates into the existing scientific python ecosystem as smoothly as possible.Obs can be the argument of most numpy [14] functions, thus allowing for an easy and fast implementation of analyses using the Γ-method.The linear propagation of uncertainties is performed to machine precision via automatic differentiation [10,15].pyerrors has already been used in a set of scientific publications [16][17][18][19][20][21].While this article is focused on data analysis in the context of lattice QCD, the pyerrors package introduced here may also be applied to other kinds of correlated data.
This paper is organized as follows.In Section 2 we briefly summarize how error propagation and error estimation are carried out in the Γ-method approach in combination with automatic differentiation.We proceed with a description of our implementation of these concepts within pyerrors in Section 3. In Section 4, we collect examples to illustrate the basic features of the package.In addition to the core routines that are used for the error computation, a number of other features has been implemented to facilitate frequently encountered tasks such as the computation of effective masses from correlation functions, fits to data with uncertainties and the solution of a Generalized Eigenvalue Problem (GEVP) [22][23][24][25].Some of these routines are presented in Section 4. Users are welcome to contribute to extending the functionality of pyerrors via the repository https://github.com/fjosw/pyerrors.

Description of the problem
Data that is obtained from a Markov Chain Monte Carlo simulation exhibits autocorrelation -subsequent measurements are not independent from each other because the underlying distributions evolve by sequential updates.Such autocorrelations in the data reduce the effective number of independent measurements.To obtain an accurate estimate of the statistical uncertainty of an observable, the autocorrelation has to be taken into account appropriately.
In the context of lattice QCD, the correct treatment of autocorrelation becomes more and more important as the field enters the precision era, where sub-percent precision is reached for phenomenologically important quantities computed at fine resolutions below 0.05 fm.Due to critical slowing down [3], a significant increase of autocorrelation times is expected and observed towards the continuum limit for all algorithms that are currently being used for large-scale simulations.Sufficiently long Monte Carlo chains and a precise treatment of autocorrelations are vital to ensure the correctness of error estimates.
Traditionally, the estimation of uncertainties of Monte Carlo observables is performed using resampling techniques such as Jackknife or Bootstrap resampling.In the presence of autocorrelation, primary data is binned to blocks in Monte Carlo time that are large compared to the autocorrelation time of the observable. 2 This may be problematic if long-ranged autocorrelation is only present in complicated derived observables and not visible in primary data that is used to determine the appropriate bin size or if bin sizes are large compared to the total number of measurements.In the latter case, only a small number of bins may be used within the resampling techniques, thus enlarging the error of the error estimation itself.
The aforementioned problems have been addressed in Ref. [5], where the so-called Γ-method has been introduced as a tool to investigate the autocorrelation function Γ(t) of primary and secondary observables and to use it to correctly scale the standard error.When combining data from different sources, e.g., Monte Carlo simulations with different parameters, the treatment within the context of the Γ-method furthermore allows to cleanly separate these sources of uncertainty in arbitrarily complicated derived observables.While this can also be done within resampling approaches, such separation is intrinsic to the Γ-method.Therefore, the management of meta-data and self-contained long-time storage in compliance with the FAIR principles for scientific data management and stewardship (Findability, Accessibility, Interoperability, and Reuse of digital assets) [26] is facilitated.

Error propagation
The general route in Monte Carlo calculations in order to obtain estimators for observables A α is to generate a sufficiently large number N of samples a i α which obey the given probability distribution and take the mean value An observable is furthermore characterized by the fluctuations which contain all relevant information for the computation of the statistical error, as the variance and its more general variant, the autocorrelation function, are invariant under constant shifts of the data.One is often interested in a derived observable which can be an arbitrary function of any number of primary observables The error propagation we use for the consequent error estimation is based on a first-order Taylor series expansion as detailed in Ref. [5] δ where the derivatives are defined as A derived observable is then again characterized by its mean f (ā α ) and its fluctuations δ i f and the distinction between primary and derived observables becomes indistinguishable [27].The numerical estimation of the derivative in eq. ( 5) as proposed in Ref. [5] may be replaced by the derivative obtained from automatic differentiation [10], making the linear error propagation exact to machine precision. 3

Error estimation
In order to estimate the statistical Monte Carlo error of any primary or derived observable one can write the estimator for the autocorrelation function in terms of the fluctuations δ i f as which for N → ∞ becomes the exact autocorrelation function.Note that Γ f (0) = var(δ i f ).The autocorrelation function is expected to exhibit an exponential decay for large t according to where the exponential autocorrelation time τ exp is the largest autocorrelation time in the system.The corresponding integrated autocorrelation time for observable f is estimated via [4] One practical problem is to find an appropriate window W for which all relevant information about the autocorrelation is captured in τ int,f while at the same time the noise contribution for large values of t does not overwhelm the signal [4].In Ref. [5] an automatic windowing procedure was suggested which minimizes the sum of statistical and systematic errors based on the hypothesis that τ exp ∼ Sτ int,f with some factor S.
The failure to detect a coupling to the slowest modes in the system in practical applications, e.g., because the Monte Carlo chain is not long enough, might lead to an underestimation of τ int,f .In Ref. [3] an extension to the estimation of the integrated autocorrelation time was suggested that accounts for the effect of τ exp in the error estimate.The idea is to attach a tail to the autocorrelation function, based on prior knowledge of τ exp .Also in eq. ( 9) one has to find an appropriate summation window W . Instead of using the automatic windowing procedure of Ref. [5], it is custom to use the first value of W for which ρ(t) = Γ f (t)/Γ f (0) is zero within its error.The error of ρ(t) can be approximated by [28] 4 The final error estimate based on the integrated autocorrelation time τ int,f can be computed via which reduces to the standard error in the uncorrelated case for which τ int,f = 0.5.

Comparison to other error estimation techniques
In this section we compare the error estimation strategy which we outlined above to two of the most used techniques in the lattice field theory community -Jackknife [6,7] and Bootstrap [8] resampling.All three approaches have in common that they are associative, so analyses can be split up in multiple smaller steps in order to obtain more complicated derived observables.In contrast to Jackknife or Bootstrap resampling the Γ-method requires the application of non-linear functions in this process only once.In turn, one requires the derivative of the result with respect to every input observable.For a large number of samples and complicated functions the computation of the derivatives can be cheaper than evaluating the function for every sample as one would do it in a traditional resampling approach [10]. 5 benefit of the Γ-method approach is that the systematic error in the estimation of the autocorrelation that comes with the truncation of the sum in eq. ( 8) decreases asymptotically as exp(−W/τ exp ).The systematic error from a binning procedure on the other hand only decreases asymptotically as τ exp /B where B is the bin size as demonstrated in Ref. [5].For sufficiently many samples N the ratio of the systematic to the statistical error becomes negligible for the Γ-method whereas the ratio stays constant in case of a binning analysis.
Another more practical advantage of the Γ-method error analysis strategy is that all the required information for a dedicated autocorrelation analysis is always available for any derived observable by construction.In the case of binning analyses combined with Jackknife or Bootstrap resampling one has to make some choice for the bin sizes of each ensemble, possibly determined from a series of standard observables.This can be dangerous in cases where very long autocorrelation times arise in specific derived observables but were not visible in the primary observables that contribute to it (for examples of such cases see Ref. [10]).A dedicated study of the evolution of covariances with the bin size including an extrapolation to infinite bin size, as outlined in Appendix D of Ref. [29], may resolve this problem in binning analyses.
One advantage of Bootstrap resampling over both Jackknife resampling and the linear error propagation we advocate is that it can capture non-linear distributions and thus non-symmetric confidence intervals.However, since a fundamental assumption of Monte Carlo simulations is the central limit theorem, a linear approximation should reproduce the correct distribution in the infinite sample limit.

Implementation
At the core of the pyerrors implementation stands the Obs class which provides the user with a new python data type for Monte Carlo observables.An Obs object is initialized with a set of Monte Carlo samples and a corresponding string which serves as unique identifier for the Monte Carlo ensemble. 6In addition the user can explicitly provide the set of configurations on which the samples were obtained also allowing for irregular or gapped Monte Carlo histories.We refer the reader to Appendix A for an explanation of our implementation of irregular Monte Carlo chains.
Elementary mathematical operations are overloaded for the Obs class.The user can, e.g., multiply an Obs object with an integer, a floating point number or another Obs object and the underlying error propagation is taken care of intrinsically following eq.( 4), where the required derivatives are obtained via automatic differentiation [10].Overloading also works for most numpy functions such as log or exp allowing the user to design an analysis without having to think about cross-or autocorrelations.Obs objects can also be arguments of iterative algorithms like non-linear least squares minimization or root finding.In order to obtain the required derivatives for arbitrary functions we use the autograd package [15].For complex valued observables we also provide a CObs class in analogy to the python complex type.
The raw results of Markov Chain Monte Carlo calculations in lattice QCD are often not single observables but Euclidean time dependent correlation functions of two (or more) operators.To conveniently handle correlation functions in analyses we provide the Corr class which can be considered an array class for Obs objects.In addition, the class has a set of regularly used methods for the manipulation of correlator objects.The user can, e.g., obtain the effective mass or discrete derivatives of a correlator by calling a class method which returns a new Corr object.The Corr class can also handle matrix valued correlation functions and may be used to solve a GEVP [22][23][24][25].Having access to data types which map onto the relevant concepts simplifies analyses, makes them more readable and less prone to error.The Corr class for example eliminates the need for any explicit loops over the temporal argument of the correlation function.

Error estimation
After arriving at a desired arbitrary function of primary observables (which can be obtained in multiple steps as explained in the previous sections) the error of this quantity can be estimated by calling the method gamma method of the Obs object.We implemented both, the automatic windowing procedure of Ref. [5] and the extended variant proposed in Ref. [3].The computationally most expensive part in estimating the Monte Carlo error along the lines described in Section 2.3 is the computation of the autocorrelation function according to eq. ( 6). 7 The computation can be sped up by making use of the Wiener-Khintchine theorem [30,31] and the real fast Fourier transform algorithm (rfft) [32] with computational complexity O(N log N ).To reproduce the exact expression of eq. ( 6) a linear convolution can be used which requires extending the vector δ f by t max zeros to ensure that there are no problems with the periodicity assumed by the rfft routine, In practice we compute the autocorrelation function up to t max = N/2.
The error estimation for a single Monte Carlo ensemble can be easily extended to analyses to which multiple, independent Monte Carlo ensembles contribute.Within a pyerrors Obs we keep track of the fluctuation of each ensemble independently and uniquely identify it with a string.The error contributions from independent ensembles can then be estimated individually via eq.( 11) and added in quadrature.This procedure allows to cleanly separate error contributions from independent ensembles and prevents contamination with spurious correlation which might arise in combined Bootstrap analyses.
Many analyses make use of information that is provided as external input.Examples for such external information in the context of lattice field theory are experimentally determined quantities such as meson masses or decay constants which are used to calibrate the theory.The treatment of (correlated) external observables via Gaussian error propagation is straightforward in the context of the Γ-method with automatic differentiation.The knowledge of the covariance matrix C for a set of external observables p i together with the Jacobian matrix for the computation of a secondary quantity O = f (p i ) which is defined via may be used to determine the contribution of the p i to the total error of O via This computation is exact in the sense of Gaussian error propagation.The information on the derivative of any secondary observable with respect to p i as well as on the covariance between the p i may be propagated through the entire analysis.The clean separation of Monte Carlo data from different independent ensembles as well as correlated constants as external input allows one to monitor the individual contributions to the squared error.

Examples
In this section we provide a few minimal examples which demonstrate the capabilities of the pyerrors framework.For a full account of the functionality see the online documentation [33].

Mathematical operations
An Obs object which carries the information of a single Monte Carlo observable can be constructed by providing two lists, one containing the samples, the other the names of the individual Monte Carlo chains on which the observable has support on In [1]: import pyerrors as pe my_obs = pe .Obs ([ sample_data ] , [ " A1 " ]) Optionally the identifiers of the configurations can be explicitly specified and tracked via the argument idl.For data on every second configuration between 1 and 2000 one could specify In [2]: If the argument idl is not given it is assumed that the first sample corresponds to configuration 1 and that there are no gaps in between configurations.An Obs object can then be modified with arbitrary elementary operations as well as numpy functions which are overloaded for the Obs class In [3]: import numpy as np derived_obs = np .log ( my_obs ** 2) After arriving at the derived observable of interest the error can be estimated via the gamma method.
In [4]: derived_obs .gamma_method () print ( derived_obs ) 0.2611(43) The default behavior of the gamma method is to use the automatic windowing procedure of Ref. [5] with the parameter S = 2.0.The user can call the method with a different value for S as argument or with an estimate of the exponential autocorrelation time via the argument tau exp in order to refine the estimation of the autocorrelation.If S = 0 is specified, the autocorrelation analysis is deactivated and the standard error is calculated instead.

Correlation functions and the effective mass
Lattice data which is defined on a set of timeslices can be represented using the Corr class In [5]: corr = pe .Corr ( list_of_obs ) Corr objects support arithmetic operations with numbers, Obs and other Corrs via operator overloading.They also implement the gamma method and other commonly used operations.As an example, we can obtain the log effective mass of the correlation function C(x 0 ) by calling the method m eff which returns another Corr object In [6]: am_eff = corr .m_eff ( variant = ' log ') am_eff .gamma_method () The new Corr object am eff has the same length as corr.The last entry, which is not defined in this case, is filled with None.We also implemented other variants of the effective mass, e.g., the cosh effective mass for periodic boundary conditions.After having identified an appropriate region in which the effective mass appears to form a plateau the corresponding value can be extracted by calling the method plateau with the desired range as argument.The method returns an Obs object and we retain access to all tools for error analysis In [7]: am = am_eff .plateau ( plateau_range = [5 , 18]) am .gamma_method () print ( am ) 0.0539(50)

Least square fitting of data
Suppose we want to fit the correlation function with an exponential decay f (t) = Z 2 e −amt .We define a fit function with two parameters as follows We can then perform an uncorrelated least squares fit by calling the corresponding class method In [9]: fit_result = corr .fit ( fitf , fitrange = [5 , 26]) From this fit result, we can get the values of Z and am and estimate their errors In [10]: fit_result .gamma_method () Z , am_fit = fit_result print (Z , am_fit ) 0.0243(10) 0.0536(46) Once again, the results are Obs objects which can be reused in subsequent analysis steps.pyerrors also supports more complicated fitting procedures like fully correlated fits or fits involving errors in the independent variables via the total least squares method.

Error contributions from different sources
For quantities with an error not originating from Monte Carlo simulations we can use the Covobs class which is initialized with a value and a squared error. 8In the following example we convert the previously extracted mass to physical units using a lattice spacing which is known as external input.

The generalized eigenvalue problem
Corr objects can also be used to represent matrices of observables defined at different time slices of the lattice.This is very useful, if one wants to solve a GEVP.We consider a correlator matrix which consists of symmetric N × N matrices at every timeslice t, where O i , O j are different operators with the same quantum numbers, e.g., operators with different levels of Gaussian smearing applied.
In order to extract the individual energy levels we can solve the GEVP for a reference time t 0 on every timeslice In [14]: evs = C .GEVP ( t0 =2) We can then obtain the corresponding eigenvalues on each timeslice via the relation C i (t) = v i (t, t 0 ) T C(t)v i (t, t 0 ) by calling the projected method and calculate the corresponding effective masses in one line In [15]: am_eff_0 = C .projected ( evs [0]).m_eff () am_eff_0 .tag = " Ground state " am_eff_1 = C .projected ( evs [1]).m_eff () am_eff_1 .tag = " First excited state " After extracting a plateau value we can visualize the effective masses together with the plateau The output of cell [16] is shown in Fig. 2. In the example above the eigenvalues are sorted for every timeslice independently.We also implement additional methods from Ref. [34] to sort the eigenvectors, if they are independently computed on different timeslices.

Data export and storage
To allow for reusability of stored data, metadata and the data itself have to be processed and stored together in accordance with the FAIR principles [26].One vital piece of information, the origin of the data, is propagated throughout an entire analysis with pyerrors.The export and import of native pyerrors structures such as Corr and Obs or lists and numpy.ndarrays of Obs is implemented via compressed JSON files.In addition to the data, these All relevant data and metadata can then be restored at any later stage within the Obs objects In [18]: Z , am_fit = pe .input .json .load_json ( ' fit_data ') The JSON serialization also facilitates the storage of pyerrors structures in relational databases.

Conclusions
pyerrors combines well-known techniques for the estimation and propagation of uncertainties in the context of Monte Carlo data.The framework allows to implement complex analyses in the existing scientific python ecosystem.This enables the user to profit from the known benefits of the Γ-method without having to care about the propagation of uncertainties, which is performed in the background.While we have explained the features of pyerrors in the context of lattice field theories, the implementation may be helpful for any kind of Monte Carlo data that exhibits autocorrelation.

Figure 1 :
Figure 1: Output of cell [13], visualizing the relative contributions to the squared total error.

Figure 2 :
Figure2: Output of cell[16], visualizing the extraction of the energy of the first excited state via a GEVP.
After estimating the error of the mass via the gamma method we can display more details about the result by calling the details method When more than one Monte Carlo ensemble or external quantity contribute to the final error we can visualize the relative contributions to the squared error in the following way A2Lattice spacing from[Ref.Y]