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 analysis of data that is based on Markov Chain Monte Carlo (MCMC) methods requires a careful treatment of the inherent autocorrelation. This is especially important in the presence of so-called slow modes that may couple to observables of interest, causing long-ranged autocorrelation. If such correlations are not properly treated, statistical errors might be underestimated, putting the significance of results in question. The Γ-method as introduced in [1,2] has been shown to decrease systematic errors in the treatment of autocorrelation in comparison to widely used techniques such as binned Jackknife [3,4] or Bootstrap resampling [5].
Despite its theoretical advantages, especially in the presence of critical slowing down [6], the Γ-method seems to be used less frequently than other resampling techniques within the lattice field theory community. Whereas the first publicly available implementations [2,6,7] do not allow to combine observables from simulations with different parameters, this is possible in the implementations of [8,9] written in Fortran and Julia. 1 In 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 [12] 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 [8,13]. 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) [14][15][16][17]. 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 [6], 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. [2], 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) [18] is facilitated.

Error propagation
The general route in Monte Carlo calculations in order to obtain estimators for observables α is to generate a sufficiently large number N of samples a i α which obey the given probability distribution and take the mean valuē 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. [2] δ 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 [19]. The numerical estimation of the derivative in eq. (5) as proposed in Ref. [2] may be replaced by the derivative obtained from automatic differentiation [8], 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 [1] τ 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 [1]. In Ref. [2] 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. [6] 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. [2], it is custom to use the first value of W for which is zero within its error. The error of ρ(t) can be approximated by [20] 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 [3,4] and Bootstrap [5] 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 is typically cheaper than evaluating the function for every sample as one would do it in a traditional resampling approach [8]. 5 A clear 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. [2]. 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 case of binning analyses combined with customary Jackknife or Bootstrap resampling one has to make some choice for the bin sizes entering the analysis at a very early stage and one might be prone to sticking to a bin size that gave reasonable estimate for a series of standard observables. This is particularly 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. [8]).
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 one underlying assumption of Monte Carlo simulations is the central limit theorem, a linear approximation seems to be sufficient in most cases.

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. 6 In 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 [8]. 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 [13]. 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 [14][15][16][17]. 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. [2] and the extended variant proposed in Ref. [6]. 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 [21,22] and the real fast Fourier transform algorithm (rfft) [23] with computational complexity O(N log N ). To reproduce the exact expression of eq. (6) a linear convolution can be used which requires to extend the vector δ f by t max zeros where t max is the maximal Monte Carlo time separation t one is interested in, 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 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 [24].

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 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.

0.2611(43)
The default behavior of the gamma method is to use the automatic windowing procedure of Ref. [2] 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 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

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 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. 8 In the following example we convert the previously extracted mass to physical units using a lattice spacing which is known as external input. 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 In [13]: m . plot_piechart () the corresponding output of cell [13] is shown in Fig. 1.

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  Figure 1: Output of cell [13], visualizing the relative contributions to the squared total error. 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 The output of cell [16] is shown in Fig. 2 Figure 2: Output of cell [16], visualizing the extraction of the energy of the first excited state via a GEVP.
In the example above the eigenvalues are sorted for every timeslice independently. We also implement additional methods from Ref. [25] 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 [18]. 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 contain information about the author, the creation date and location and the version of pyerrors. The file content may be further described by a python string, list or dictionary 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. the available configurations or if data on single configurations is missing. To generalize the notation in the main text, let us consider a set of N subsequent configurations E = {1, 2, . . . , N − 1, N } and the subset of configurations E α ⊂ E, which has been used for the measurements of the primary observable a i α . For the correct bookkeeping of configurations, we introduce the quantity which may be used to compute the number of configurations in E α by We now define the mean value, previously introduced in eq. (1), bȳ generalizes eq. (4). To account for the reduced number of configurations in the estimation of the autocorrelation function, we need to redefine Γ f (t) from eq. (6) with an adjusted norm N α (t) which replaces the factor N − t according to Equations (10) and (11) are then defined with the normalization factor N α which gives the number of configurations that enter the analysis instead of N .
In the limit E α = E all formulae in the main text are recovered.
In terms of the implementation in pyerrors, eqs. (12) and (13) have to be adapted via , (A.9) to allow for the inclusion of irregular Monte Carlo chains. The rfft is employed to compute N α (t) efficiently.