{'fraction_non_alphanumeric': 0.033645693886657745, 'fraction_numerical': 0.012940651494868362, 'mean_word_length': 4.711693342848404, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The long road to calibrated prediction uncertainty in computational chemistry

Uncertainty quantification (UQ) in computational chemistry (CC) is still in its infancy. Very few CC methods are designed to provide a confidence level on their predictions, and most users still rely improperly on the mean absolute error as an accuracy metric. The development of reliable uncertainty quantification methods is essential, notably for computational chemistry to be used confidently in industrial processes. A review of the CC-UQ literature shows that there is no common standard procedure to report nor validate prediction uncertainty. I consider here analysis tools using concepts (calibration and sharpness) developed in meteorology and machine learning for the validation of probabilistic forecasters. These tools are adapted to CC-UQ and applied to datasets of prediction uncertainties provided by composite methods, Bayesian Ensembles methods, machine learning and a posteriori statistical methods.


I. INTRODUCTION
As stated in recent perspective articles 1-3 , uncertainty quantification (UQ) in computational chemistry (CC) is still in its early stages of development. For instance, in electronic structure theory, at the exception of the BEEF-type methods 4 , none of the methods implemented in popular CC codes provides an uncertainty or a confidence index on the calculated properties. Confidence is generally based on benchmark studies which provide performance indices, such as the ubiquitous mean absolute error (MAE). However, it is well established that, lacking a probabilistic interpretation, the MAE should not be used as an uncertainty proxy [5][6][7] . As will be seen below, prediction uncertainty (PU) might be much more complex to estimate than usual performance indices.
This difficulty to quantify our confidence in model predictions has far reaching consequences.
It is an obstacle for CC to stand on par with, or to replace, physical measurements 8 or to be used in decision making 9 . It has also a strong impact in multi-scale simulation, where propagation of uncertainty through the scales is necessary to assess the reliability of predictions [10][11][12] , or in iterative learning to minimize the cost of high-level calculations [13][14][15] . In benchmarking, a very sensible way to select among a set of levels of theory would be to pick one with a fit-to-purpose PU. It is also worth to note that uncertainty provides a metric for the comparison of measurement values, a prerequisite to the production of reproducible results 16 . In all such applications, a PU estimate has to be fair : underestimation is potentially dangerous (overconfidence), and overestimation is wasteful. Reaching a good balance is the main challenge of UQ.
For the purpose of the present study, I sorted UQ approaches in two classes, according to their embedding level within the CC method: • Embedded UQ methods produce PU concurrently with property predictions. This is an heterogeneous class which presently encompasses the above-mentioned BEEF or Bayesian Ensemble approach 4 and several bottom-up correction methods with detailed uncertainty budgets (the Type B methods of Ruscic 5 ), such as the Feller-Peterson-Dixon method 17,18 or the ATOMIC protocol [19][20][21] . Some machine learning methods, such as Bayesian neural networks which are designed to provide PUs along with prediction, also belong here 22,23 .
These methods are akin to the "molecule-specific" UQ concept 3 .
Very often, the validation of PU values provided by these CC-UQ methods is based on the qualitative appreciation of the agreement of prediction uncertainty with the amplitude of errors 28,38 .
For the Bayesian ensemble methods, calibration has been assessed by a visual examination of the normality of z-scores histograms 39 , or by visual evaluation of the width of error bars 4,40 . When quantitative tools have been used, notably for the validation of a-posteriori methods, a popular test statistic is the prediction interval coverage probability (PICP) 23,41 . PICP tests have been done either on the full dataset 21,29,42 , or using a splitting scheme between calibration and validation sets 6,42 . Globally, there does not seem to be (yet) a consensus in the community on the adequate validation vocabulary, concepts and tools, which would be necessary to compare the merits of different CC-UQ strategies.
In the past decades, probabilistic forecasting for meteorology has been the object of fundamental developments of validation concepts and methods 43,44 . Two major concepts resulting from these studies are calibration (reliability), and sharpness (resolution). A probabilistic prediction method is said to be calibrated if the confidence of predictions matches the probability of being correct for all confidence levels 45 . A calibrated method is sharp if it produces the tightest possible confidence intervals 46 . Sharpness is conditional on calibration (a method cannot be sharp if it is not calibrated). With the emergence of PU estimation in Machine Learning, several calibration and sharpness metrics, and graphical checks are now also used in this field 22, 23,46,47 .
Most CC methods are based on deterministic algorithms, but our lack of knowledge about their prediction errors links them to probabilistic forecasters. Part of these validation tools have recently been introduced to the CC field by Tran et al. 22 to compare the calibration and sharpness of several ML algorithms trained to predict adsorption energies. This is, to my knowledge the only study of this kind in the application field of interest. My intent in this study is to adapt and apply calibration/sharpness validation methods to a wider CC domain, and to evaluate their pertinence in various scenarios, in the expectation that such methods could be more generally used in the community.
In the next section (Sect. II), I review the specifics of CC-UQ (error sources, UQ methods) and previous validation practices. I present then (Sect.III) the panel of validation tools for probabilis-tic forecasters (calibration and sharpness metrics and graphical checks), and I propose adapted methods for the typical CC-UQ outputs. In Section IV, these tools are applied to a variety of datasets. The main features of this study are reported and discussed in the Conclusions (Sect. V).

II. THE CC-UQ CONTEXT
A major obstacle to UQ in computational chemistry is the predominant role of systematic errors.
As uncertainty is a non-negative parameter expected to quantify unpredictable errors 48 , estimation (and correction) of predictable or systematic errors is a challenging, but necessary preliminary step to CC-UQ. The international reference guide for metrology (a.k.a. 'the GUM' 48 ) states that "It is assumed that the result of a measurement has been corrected for all recognized significant systematic effects and that every effort has been made to identify such effects." (GUM, Sect. 3.2.4). As a consequence, the partition between unpredictable and predictable errors depends on the efforts that a modeler is ready or able to invest in the analysis.
I consider next the main sources of errors in CC (Sect. II A) and the existing UQ methods, based on different approaches to systematic error corrections (Sect. II B).

A. Error sources in Computational Chemistry
The main error sources in CC have been reviewed recently 1,2 , and can be tagged as numerical, parametric, and model errors. I come back briefly on these categories in order to discuss the expected errors distributions, which is an essential ingredient to define UQ validation methods.

Numerical Errors
These include errors due to finite arithmetics implementation of CC codes and, for stochastic methods, to the random errors resulting from finite sampling. For properly implemented and converged methods, numerical errors can often be assumed to be well controlled and negligible against other error sources 8 (except for instances of numerical chaos 49,50 ). However, for most algorithms we are still missing rigorous error bounds 51,52 , and estimation of the amplitude of numerical errors due to finite arithmetics requires computational approaches. In the Monte Carlo framework, one observes the effects of tiny perturbations of model inputs 8,49 and/or arithmetic operations [53][54][55] , without need to alter the codes. The use of interval arithmetics has also been proposed 56 , but it might require in-depth recoding, barely an option for legacy CC codes. This is however not a fatality for newly developed codes, as shown recently by Herbst et al. 52 with the development of a DFT code in the Julia language accepting arbitrary floating point types, including intervals. For some parameters, the amplitude of numerical errors can also be estimated by systematic convergence studies 57 .
One can consider numerical errors as random variables, but their distribution is not necessarily normal. Non-normal distributions occur, for instance, for rounding errors in floating point arithmetics 58 , or stochastic errors for QoIs in molecular dynamics 9 .

Parametric errors
These occur in methods involving the statistical estimation of parameters with respect to reference data (e.g., semi-empirical force-fields, extrapolation schemes in composite methods, semiempirical methods, statistical corrections...). An essential property of parametric errors is that their amplitude should decrease when the size of the reference dataset increases. The amplitude of parametric errors is typically estimated by Monte Carlo uncertainty propagation, from a probability density function of the parameters obtained, for instance, by Bayesian inference 59,60 .
The distribution of parametric errors depends on the distribution of the uncertain parameters and on the functional form of the model with respect to these parameters. Here again, normality cannot be assumed as a default distribution feature.

Model errors
For computational chemistry, model errors include level-of-theory errors (e.g., the choice of a density functional approximation, DFA) and representation errors (e.g., the choice of a basis set or discretization grid). They are often the dominant contributions to a CC uncertainty budget.
Compared to numerical and parametric errors, which are mostly aleatoric, model errors are systematic. The estimation and correction of model errors is usually handled either by bottom-up approaches of physics-based corrections 5 , or by corrective statistical models (so-called discrepancy functions 32 ) involving the comparison with a dataset of reference results 61 . It is important to note that, in contrast to parametric errors, model errors do not decrease when the size of the reference dataset increases.
Model errors can present any type of distribution 7 . However, a recent study on the shape of error distributions for different QoIs and DFAs showed that the distributions are most often unimodal, with various levels of asymmetry, and that they are in most cases heavier-tailed than a normal distribution 62 .
B. Main error correction and PU estimation methods I review here the main methods appearing in the CC-UQ literature to correct systematic errors and estimate prediction uncertainty, focusing on the level of information they are able to provide, which goes form a standard uncertainty to a full distributions of predictions.

Bottom-up correction methods
Proceeding from a chosen level-of-theory/representation and applying successive corrections based on physical criteria is a systematic way to establish an uncertainty budget: each correction, being imperfect, coming with some uncertainty. This kind of approach is referred to as "Type B methods" by Ruscic 5 , in reference to the method defined in the GUM 48 for the estimation of uncertainty in absence of a statistical sample (statistical analysis of data samples defines Type A methods 5,63 ). It appears however that some form of Type A analysis lies at the heart of some bottom-up correction methods. For instance, the Feller-Peterson-Dixon method considers two main uncertainty components: the complete basis-set (CBS) and the zero-point energy (ZPE) corrections. For the CBS correction, uncertainty is estimated as the half-spread of the corrections provided by five different extrapolation models 17 . For ZPE, it is taken as half the difference in the last step of a series of approximations 64 . I note for further reference that in the FPD approach uncertainties are combined linearly instead of quadratically, a worst case scenario ignoring possible error compensations 65 . This is expected to provide a "crude but conservative uncertainty", from which the estimation of prediction intervals is not straightforward and complicates the interpretation of calibration tests.
In contrast, the ATOMIC composite method aims "to provide realistic corrections and uncertainty estimates corresponding to intervals of 95% confidence..." 19 . Here, a mixture of linear and quadratic uncertainty combination is used in deriving the final uncertainty. Bakowies 20 summarizes the interest of the bottom-up correction approach: "The proposed model is a welcome alternative to statistical assessment, first because it does not depend on comparison with experiment, second because it recognizes the expected scaling of error with system size, and third because it provides a detailed account of the importance of various contributions to overall error and uncertainty."

A-posteriori PU estimation
A-posteriori approaches elaborate a statistical model independent from the CC method to estimate its PU, generally after correction of the original model predictions for systematic errors 6,24,66 .
As for bottom-up corrections, different PU estimation models are necessary for different predicted properties by a given CC method.
The a-posteriori approach requires three elements: 1. A reference dataset on which to calibrate the statistical model, the ideal being a large set (at least several hundreds) of high-accuracy values.
2. A statistical model to correct the (systematic) trends in the model errors. This optional step might involve from a simple shift to machine learning models.
3. A statistical method to estimate the prediction uncertainty from the residual errors, the complexity of which depends on their distribution.
If a trend correction is applied (step 2), it is important to note that the a-posteriori approach does not provide a PU for the original CC method, but for its corrected version. This has been clearly illustrated for the correction of vibrational frequencies by scaling factors 24 , or for the linear correction of DFT-predicted Mössbauer isomer shifts 29 .
One can point out two main obstacles to the success of the a-posteriori approach, affecting the first and third steps: • The reference dataset should contain as many as possible high-accuracy data, which, depending on the studied property, might be difficult to achieve, notably for experimental data. For the application of simple correction models, the errors dataset should also be homogeneous, in the sense that it should not contain multiple contradictory trends preventing the success of step 2. This might be less stringent for machine learning correction models, which should be able to correct for multiple trends, provided a relevant set of input features. Also, care should be taken that the prediction uncertainty is corrected from reference data uncertainty 6,63 , a problem absent from bottom-up correction methods.
• The ability to establish reliable PUs relies heavily on the distribution of the corrected errors. Skewed or heteroscedastic error distributions might be specifically challenging for the definition of an uncertainty. The difficulty to establish and communicate asymmetric and/or heteroscedastic uncertainty might even be considered as a criterion to reject some CC methods for PU estimation.
A-posteriori methods have been mostly used to estimate standard uncertainties or prediction intervals. Depending on the approach, a full distribution might be accessible (e.g., for regression models, Gaussian processes...).

Bayesian ensemble methods
In the DFT framework, the Bayesian Error Estimation density Functional (BEEF) family of methods has been designed to provide uncertainty on its predictions 39,67-70 .
In a first step, the parameters of a DFA are calibrated by Bayesian inference against a reference dataset. However, the resulting parametric uncertainties are typically insufficient to cover the amplitude of prediction errors (remember that the amplitude of parametric errors decreases with the size of the calibration dataset).
To restore the validity of the statistical model, the variance-covariance matrix of the reference data is scaled by a 'temperature' factor such that the mean prediction variance matches the variance of prediction errors 40 . The resulting probability density function (pdf) of the DFA parameters is captured as an ensemble of parameters values which can then be used to estimate PU on any relevant property by Monte Carlo sampling. Note that, as Bayesian ensembles are calibrated to cover the amplitude of prediction errors, their prediction uncertainty is affected by reference data uncertainty. When comparing predictions to experimental values, care should be taken not to count experimental uncertainty twice.
In the available literature the results of BEEF models are summarized by a mean value and a standard deviation and the ensembles are not used for validation.
UQ is central to many automatized applications in screening or design of efficient compounds, and several ML algorithms are now available to estimate prediction uncertainty, either through distribution-or ensemble-based methods. A recent comparison of the calibration and sharpness of these methods by Tran et al. 22 reveals a wide spectrum of reliability of ML-UQ methods.

III. METHODS FOR THE VALIDATION OF PREDICTION UNCERTAINTY
The general framework for the validation of probabilistic predictions is first presented to define the calibration and sharpness concepts and the associated statistical test and metrics. The cases with restricted information, in the shape of an uncertainty (standard or expanded) instead of a full distribution, are then considered, for which validation methods are derived from the general framework. The concepts are illustrated on synthetic data.

A. General case: probabilistic predictions
A probabilistic prediction provides a distribution over the values that can be taken by the quantity of interest (QoI), V . The predicted cumulative density function (CDF) for V is noted F (V ) and the reciprocal quantile function F −1 (p) ≡ F −1 p , where p is a probability. For the sake of validation, predictions are made for a series of M test systems for which one has reference values . For each reference system i, one has thus a predicted CDF F i (V ) and the corresponding quantile function F −1 p,i . Few UQ methods provide prediction distributions, and approximations of CDF functions are often available through ensembles of values, representative of the predictive distribution. Such ensembles are generated by ML-UQ methods, Bayesian Ensemble BEEF-type methods or stochastic algorithms.

Definitions: calibration and sharpness
A method is considered to be calibrated (or reliable) if the confidence of predictions matches the probability of being correct for all confidence levels 22,45 . Formally, the empirical and predicted CDFs should be identical 46 where 1(x) is the indicator function for proposition x. This equation can be generalized to pre- where is the 100p % prediction interval for the predicted value V i .
However, this calibration condition is not sufficient to ensure that the prediction uncertainties or confidence intervals are useful. Eqs. 1-2 provide a global calibration assessment over the test set, but does not guarantee that calibration is realized locally for each predicted value. A useful prediction model must also be sharp, i.e., prediction intervals should be as tight as possible around any predicted value 46 . There is no formal definition for sharpness as there is for calibration.

Example
Let us consider a toy model where the predictive CDF F (V ) is learned from an ensemble of errors E i = R i − V i with uncorrected bias and without accounting for the linear dependence of E on V ( Fig. 1(a)). Assuming that the test set has the same properties as the learning set, using the system-independent F −1 p function will provide correct prediction intervals over the test set, and Eqs. 1-2 will be satisfied. However, if F −1 p is used to design predictive intervals for a local subset of the test set, for instance the positive values of V (( Fig. 1(d)), the intervals will clearly be oversized and ill-centered: sharpness will not be ensured. In the first case, calibration tests such as the PIT histogram presented below (Sect. III A 2) do not detect a problem ( Fig. 1(b)), while the same test on the local subset of the data is diagnostic of a calibration problem ( Fig. 1(e)).

Calibration tests and metrics
Several visual checks and statistics have been proposed to estimate calibration and sharpness of probabilistic predictions 43,44,46,79 . I review below the most pertinent ones for CC-UQ applications and extend the toolbox with a graphical representation enabling to assess simultaneously calibration and sharpness. Some of these tools have been used used in a ML setup, where statistical uncertainty on the validation statistics is negligible; this is not the case in a typical CC-UQ problem, so I added confidence intervals to all estimators. Probability integral transform histogram. The probability integral transform (PIT) is the value that the predictive CDF attains at a test value R i

Example (continued)
PIT histograms are shown in Fig. 1(b,e). For the first one, the binned PIT values stand in the confidence range for a uniform histogram, assessing the calibration of the predictions. A contrario, the second histogram presents significant deviations from the uniform, with an excess of small values pointing to a negative bias in the predictions.
Calibration curve. Calibration can also be checked by comparing the estimated success rate in the left hand side of Eq. 1 A calibration error (CalErr) is also proposed in the literature, as the sum of squared differences over the probabilities in the calibration curve 46 , or its square root 22 . The latter is retained here for its homogeneity with an error. As MisCal, it can be used to compare several methods on their calibration level.
Example (continued) Fig. 1(c,f) shows the calibration curves corresponding to both scenarios in Fig. 1(a,b).
The one corresponding to the full calibration dataset does not deviate notably from the identity line, and the MisCal statistic is much smaller that its upper limit MisCalUp.
In the second case, there is no ambiguity about miscalibration, either for the curve or the statistics. As a confirmation, the calibration error in the second case (CalErr= 1.3) is much larger than in the first case (0.048).
PICP testing for a series of target coverage probabilities. The effective coverage of a 100p % prediction interval I p is estimated by its prediction interval coverage probability (PICP) 23,41 , as the ratio of the number of successes S p , to the size M of the validation set where To test if a PICP value ν p is compatible with the intended coverage probability p, one checks if p lies within a 95 % confidence interval around ν p (see Section III D 1 for implementation details).
Direct application of Eq. 2 would lead to test PICP values for a regular set of coverage probabilities in [0, 1]. However, one is typically interested in large coverage values, and there is little practical interest to consider low-probability intervals. In the following, I will consider p = 0.5 (inter-quartile range), 0.75 and 0.95.

Sharpness tests and metrics
Local Coverage Probability (LCP). In order to evaluate sharpness, one can check local calibration by performing PICP tests on subsets of the validation data. The simplest scenario is to split the dataset into contiguous areas 82 of the predictor variable (typically the QoI V ), but other splitting schemes can be considered. It might also be interesting to use the prediction uncertainty, or a proxy, as predictor variable to check its impact on local PICP estimates.
Two constraints have to be considered: (i) the subsets sizes have to remain large enough for PICP testing to have some power, and (ii) using equi-sized subsets could simplify the appreciation of multiple tests. The first constraint limits the resolution of the LCP sharpness evaluation, and for small datasets, rather than splitting the dataset into small sets, one might use overlapping subsets. Whatever the splitting scheme, the LCP method proceeds as follows: 1. Choose a set of equi-sized areas of the predictor variable (contiguous or overlapping).

Repeat for several values of probability p:
(a) within each area (indexed by k), estimate the local PICP, ν p,k = S p,k /M k (b) plot ν p,k and its 95 % confidence interval at the center of each area.
The LCP plot enables to check (1) where Var(F i ) is the variance associated with the CDF at point i, F i (V ). sha has the dimension of an uncertainty on the QoI and it corresponds to the mean prediction uncertainty (MPU) used in an earlier study by Pernot et al. 6 . The mean predictive variance (sha 2 ) has also been used as a model performance metric by Proppe at al. 29 , and as a calibration statistic for the BEEF methods 4,40,67 .
It is considered that sha should be small, but no threshold value is available to distinguish between sharp and unsharp methods. sha is therefore a convenient metric to compare several methods, but not for self-standing sharpness assessment.

B. Validation of expanded uncertainty
The majority of CC-UQ methods does not provide a prediction CDF (F i ), but limited summary statistics, such as sets of predicted values V i ± u V,i . In order to use validation methods based on Eqs. 1-2, one would have to make assumptions on the underlying CDFs. As we have seen that there is no typical shape for the CC errors distributions (Sect. II A), it is preferable to avoid such assumptions and derive distribution-free validation methods. Let us see how the elementary test in Eq. 2 can be implemented in such cases.
Assuming the symmetry of the prediction interval, one can write where U p,i is the expanded prediction uncertainty at the p level. Eq. 8 is therefore equivalent to Note the U p should account for both calculation and reference errors.
So, when an expanded prediction uncertainty for the errors -typically U 95,i -is available, on can consider Eq. 10 without hypothesis on the errors distribution other than symmetry and the ones made to estimate U p,i . This enables us to use PICP testing (Sect. III A 2) and the LCP analysis for sharpness.

C. Validation of standard uncertainty using z-scores
When a standard uncertainty u E,i is available, Eq. 10 can be rewritten as where f p,i is a coverage factor associated with the errors CDF F i , at point i. Considering that the errors distribution depends strongly of the dominant error source in a calculation, it is often difficult to define the errors CDF on which to estimate coverage factors. In such cases, it might be interesting to test directly the consistency of the errors and uncertainties through their ratio, as shown below.
Let us assume that we have unbiased errors with unknown distribution, but known variance Then the z-scores, z i = E i /σ i , are unit-scaled and zero-centered variables, i.e., Therefore, if the errors are unbiased and the σ i values are correctly estimated by u E,i , the distribution of z-scores should be unbiased with unit variance: In absence of information on the z-scores distribution, testing the value of Var(Z) cannot be done using the popular chis-squared test, which is not robust with respect to deviations from the underlying normality assumption (Appendix B). As for the PICP, the z-scores variance is tested against a 95 % confidence interval (estimated by bootstrapping; see Sect. III D 2).

Local z-scores variance (LZV) analysis
Similar to the LCP analysis for coverage probabilities (Sect. III A 3), I propose a sharpness estimation analysis based on the variance of z-scores, estimated in a series of contiguous or overlapping areas (LZV analysis). For each subset, Var(Z) and its 95% confidence interval are estimated and plotted on a graph. Intervals that do not cross the unity variance line might be considered as problematic, but here also, trends are of diagnostic interest.

Example
A dataset of errors has been generated in the hypothesis of a Student's-t errors distribution (T (ν)) and non-uniform random variances issued from a chi-squared distribution (χ 2 (ν)) scaled to have unit mean The errors and z-scores distributions for a sample of N = 1000 z-scores are shown in Fig. 3(a,b). The z-score statistics are E(Z) = −0.045 (30) and Var(Z) = 0.894 (53), with a 95 % probability interval equal to [0.8, 1.0] for the variance, which is validated. The LZV analysis has been performed along the calculated values (Fig. 3(c)) and along the prediction uncertainty ( Fig. 3(d)), using five areas of 200 points. Both graphs show the compatibility of the z-scores with the unity variance requirement, and no remarkable trend is observed.
D. Implementation

Binomial proportions confidence intervals
As presented in Appendix A about PICP testing, the discreteness of binomial proportions ν p and the asymmetry of the associated confidence intervals make the estimation of binomial proportions confidence intervals (BPCI) a complex problem. Numerous methods are available, from which I retained the continuity corrected Wilson method for this study 83 . 15. If the number of subsets is smaller than 5, a sliding window of size M/n is used.

Testing z-scores variance
In order to test the variance of z-scores with respect to the target value (1), Appendix B concluded on the use of a bootstrapping method, with the limits of a 95 % confidence intervals estimated by the BC a method 84 . The design parameters for the LZV analysis are the same as for the LCP analysis.

IV. APPLICATIONS
The tools and statistics presented above for the validation of calibration and sharpness of PU estimates are applied below to several datasets extracted from the literature or provided by kind colleagues. The choice of datasets is subjective, mostly guided by availability and complementarity in methods and properties. They cover the major embedded and a-posteriori CC-UQ methods.

A. Datasets with expanded prediction uncertainty
Despite the recommendations of Ruscic 5 for the systematic use of U 95 expanded uncertainty in reference databases 63 , I did not find many CC-UQ studies providing them explicitly for a reasonable dataset size. I selected two recent studies in which I found considerations about prediction sharpness.

Expanded uncertainty of errors
The ideal scenario is to have sets of predictions and reference data with their expanded un- It is then directly possible to test the PICP value leading to a single assessment of calibration at the p = 0.95 level.
When the uncertainties of the reference data are negligible, this approach tests directly the quality of the expanded uncertainties for the predicted values U (V ) 95 . Otherwise, the quality of the reference set uncertainties U (R) 95 will also affect the results. They did not identify any concern about the predicted prediction intervals and conclude that their method "provides a fair estimate of 95 % confidence" 21 . There is no uncertainty on the reference values for this dataset. The errors distribution is shown in Fig. 4(a).

The BAK2021 dataset
Considering that expanded U 95 uncertainties are reported, I tested only this probability level.
The global PICP ν 0.95 = 0.92(3) is not significantly different from the target value.
Three alternative versions of the LCP analysis are tested. The LCP areas have been chosen along the calculated ZPE values (Fig. 4(b)), the estimated prediction uncertainty (Fig. 4(c)) and the heteroatoms content of the molecules (Fig. 4(d)). Because of the small dataset, a sliding window of width M = N/2 was used in all cases. Despite the large error bars, a few points indicate a statistical rejection of local PICP values. In ZPE space (Fig. 4(b)), one notices a global trend, with PICP values increasing from 0.89 to 0.96 as the ZPE value increases. This might reveal an underestimation of the uncertainties for small ZPE values. An analog trend is observed in PU space (Fig. 4(c)), with PICP values increasing from 0.84 to 1. This would indicate that small prediction uncertainties are too small, and large ones too large. The third representation is along the property that was used to design the PU model (Fig. 4(d)), and the positive trend, although still present, is much weaker than in the ZPE or PU space.
Overall, one can agree with Bakowies and von Lilienfeld that their UQ method globally provides "fair" U 95 values, but with a caveat about a sharpness issue, as revealed in the PU-space LCP analysis.

The PRO2021 dataset
In a recent study, Proppe and Kircher 42 estimated prediction uncertainty for reaction rates In such conditions, the coverage statistics cannot be used to differentiate the models, and due to the saturation of the PICP values (very close to the upper limit), the LCP analysis does not provide useful information (not shown).
Considering that the expanded uncertainties were derived by Proppe and Kircher from standard These tests bring evidence that the polynomial mode (b) is a notable improvement over the uniform model (a). However there remains a global overestimation of the uncertainties, and the lack of sharpness is not fully resolved.

B. Datasets with standard prediction uncertainty
I consider here datasets where the coverage level of PU is not explicitly provided.

Standard uncertainty of errors
The most common scenario in the CC-UQ literature is based on standard uncertainties of the predictions and reference data: V i , u V,i , R i and u R,i . The uncertainty on the errors is directly accessible through the combination of variances enabling the computation of z-scores without further hypothesis, to be tested by their mean and variance.
Again, non-negligible uncertainties on the reference data might affect the results. If the reference uncertainty is missing, it can be taken as null (for instance in the case of high-accuracy calculated values) or as constant, if a typical value can be found.

The FEL2008 dataset
This dataset has been extracted manually from For the reference data, in absence of specific information, I assumed that expanded uncertainties were used (U (R) 95 ). The process of PU estimation in the FPD method has been summarized above (Section II B 1), and it is difficult to infer its nature, beyond a possible over-estimation implied by the worst-case scenario strategy.
As an initial assessment, I estimated that FPD PUs are close to standard uncertainties and computed u E values as u E,i = u 2 V,i + U The distribution of errors with their uncertainty is shown in Fig. 6(a). SiH stands as a strong outlier, with an uncertainty too small to cover its large error level. The case is not discussed in the original article, although it is stated that the heat of formation of silicon was not well established at the time of publication. However, there does not seem to be a systematic error, as other Sicontaining systems are not visibly affected. Other outliers such as BN or B 2 have error bars large enough to cover their deviation.
A plot of the z-scores distribution after removal of the SiH system ( Fig. 6(b)) reveals two points with outstanding values: the atomization energies of S 2 and CF 2 O have apparently underestimated uncertainties. The histogram of the z-scores (Fig. 6(b), left panel) shows a strong concentration of small z-scores values and a strong departure from a normal shape. Nevertheless, the variance of the z-scores is close to 1: Var(Z) = 1.04 (33), which does not invalidate the derivation of u E . The LZV analysis against prediction uncertainty is shown in Fig. 6(c) for a sliding window of width M/2, showing no significant deviation from the target, considering the large uncertainties. There is a slight positive trend, with a step around 0.5 kcal/mol, which is linked to the prevalence of negative z-score values for larger PUs (Fig. 6(b)).
Unless my treatment of the original data is unduly lucky, it appears that the "crude and hopefully conservative" appreciation of the original authors 18 , might be self-deprecating, but more data would be necessary to conclude. I propose here to check how the proposed validation methods perform with this dataset.

The PAN2015 dataset
The error set and the corresponding z-scores are plotted in Fig. 7. Both distributions deviate from normality. Besides a small trend in both quantities, the main feature in these plots is the heterogeneity of the z-scores distribution along the predicted values, with a noticeable negative  Fig. 7(c)). The LZV analysis with respect to the prediction uncertainty shows the opposite trend, decreasing from a value above 3 to about 0.5 for the larger PUs (Fig. 7(d)). Briefly, small uncertainties are underestimated, while large ones are overestimated.

The PAR2019 dataset
A dataset of 35 harmonic vibrational frequencies (PAR2019) has been extracted from an article by Parks et al. 87 (Table I, columns µ and σ), as another example of BEEF-generated uncertainties. This is a small dataset, which puts the validation methods to their limits. The errors are plotted in Fig. 8(a), and the corresponding z-scores in Fig. 8(b).
On Fig. 8 This clearly questions the calibration and more essentially the sharpness of BEEF-derived PUs.

C. A-posteriori UQ methods
A-posteriori prediction uncertainty has to be estimated by a statistical model using the differences between model predictions and reference data, and complementary uncertainty information, when available. The main target is to provide reliable PU in the shape of standard and/or expanded uncertainty.
The first step in PU estimation is the correction of trends in the errors dataset. The estimation of trends relies typically on the calculated value, V , as a predictor variable 6,24,25,29 , but more complex scenarios can be considered, as might be done in QSAR 88 and ML 22,78 methods.
The workhorse model of trend correction is the low-degree polynomial. I note below PTCn (Polynomial Trend Correction of order n) the correction of the errors trend vs. V by a polynomial of degree n. For this study, the PTCn model's parameters and their uncertainty are estimated by standard least-squares, which is known to be robust to non-normal error distributions 89 . More complex PU models could be used for datasets presenting complex trends, requiring weighting schemes (e.g., in the case of heterogeneous reference values uncertainties 6 ) or the consideration of heteroscedasticity 21,42 . More specific trend correction models are defined in their application case.
Once trends have been corrected, one is left with a set of residual errors which can be considered as unpredictable and can be treated as random variables.
I have considered in this study three methods (DIST, PRED, EQ) to estimate either prediction uncertainty or the limits of 100p % prediction intervals. All these methods assume the homoscedasticity of the errors: • DIST: based on mean (MSE) and standard deviation of the errors set (RMSD) provides a standard uncertainty The interest of this model is to enable the user to infer a PU from statistics often reported in benchmark tables (MSE and RMSD). Note that for well corrected trends, one should have µ 0. In order to avoid hypotheses on the errors distribution, the use of z-scores based validation methods is best suited to this case. Note that, by construction, one will get Var(Z) = 1, and the focus should be on the sharpness assessment by LZV analysis.
• PRED: based on the least-squares based statistical predictions of the trend correction model 6,28,29 . The least-squares formalism provides standard prediction uncertainty or expanded uncertainty at any probability level. The prediction intervals account for the parametric uncertainty of the correction polynomial, but are symmetrical and derive from the expansion of the prediction uncertainty by Student's-t factors. This approach can be validated for calibration and sharpness by PICP-based and z-scores-based methods.
• EQ: based on Empirical Quantiles q p estimated from the errors set. This model makes no distribution hypothesis, except that the errors ECDF is a good proxy for the CDF of prediction errors F . It ignores the parametric uncertainty due to the trend correction model, but enables to treat non-normal errors distributions and to estimate 100p % prediction intervals from quantiles of the errors distribution This approach guarantees a good calibration, but not sharpness, as the global distribution might not be locally optimal in presence of heteroscedasticity. The EQ method can be validated by the full arsenal based on PICP estimation. A symmetrized expanded uncertainty, noted SEQ p can be defined as the half range of I p .
Cross-validation. For a-posteriori methods, one can benefit from the ability to estimate a PU repeatedly to enrich the LCP analysis with cross-validation. Schematically, one proceeds as follows: 1. Cross-validation (repeat until all points have been tested) (a) Split randomly the data into learning and test sets.
(b) Use the learning set to estimate prediction intervals.
(c) Test these intervals on the test set (does the interval contain the test value or not ?).
2. Dispatch the test results into the chosen areas of the predictor variable, and estimate the PICPs as the percentage of positive tests in each interval ν i .
Leave-One-Out (LOO) or k-fold cross-validation can be used in step 1. To improve the estimation of ν, k-fold cross-validation can also be repeated several times on randomly reordered datasets.

The PER2017 dataset
The most famous a-posteriori method might be the scaling of harmonic frequencies. A few years ago, I showed that it had the notable interest to enable the estimation of a prediction uncertainty 86 .
As an example, I used a set of 2278 frequencies calculated at the CCD/6-31G* level, extracted from the CCCBDB 90 . To check calibration, I then estimated the 2-σ PICP, which ideally was 95 %, and much better than for a concurrent method. I propose here to revisit this dataset through the calibration/sharpness prism.
The distribution of errors after scaling is presented in Fig. 9(a), from which one can make three observations: (1) in the absence of the constant term in the correction model, the scaling does not fully correct the bias (there is still a very small trend in the errors, which is probably irrelevant); (2) the histogram is far from being normal; and (3) its symmetry is not perfect. More problematic, the dispersion of the errors does not seem to be uniform along the predictor axis. One might thus expect less than optimal sharpness...
Let us consider first the z-scores analysis with an PU derived by the DIST approach ( Fig. 9(b)).
By construction, the variance of the z-scores is 1, but the LZV analysis reveals a sharpness issue, with several areas with small variance values deviating significantly from the target.
Not surprisingly, considering the shape of the errors distribution in Fig. 9(a) A much better calibration is automatically obtained by using prediction intervals based on the EQ approach (Fig. 9(c)), but the issue of sharpness becomes prominent, as some local PICP values deviate significantly from their targets. Note however that the situation is not so bad at the 0.95 level, and one would still get a reasonable estimation of a uniform U 95 value (70 cm −1 ), depending on the intended application.
Nevertheless, this questions the reliability of the PU derived by the scaling procedure for this dataset, which features multiple underlying trends. More elaborate, bond-based, scaling methods 91,92 might enable to derive sharper predictions.

The PAN2015 dataset: an alternative view
For comparison with the original PU estimation, I applied an a-posteriori analysis to the PAN2015 dataset analyzed above (Sections IV B 3). Considering the trend observed in the errors, I used a linear trend correction (PTC1). The distribution of residual errors is shown in Fig. 10(a). The distribution is still not perfectly normal.
Using the PRED method, the PICP values 0.55(3)/0.77(3)/0.93 (2) are nevertheless consistent with their targets Fig. 10(b). Although the LCP analysis does not provide evidence of sharpness issues, there seems to be a small systematic bias on the local PICP values that could be compensated for by using the EQ approach Fig. 10(c).
The mean prediction uncertainty estimated by the PRED model is about 0.13 kcal/mol. By comparison, the BEEF-based PUs range from 1/4 th to 3 times this value, with about 67 % of them in excess.

ML and ∆-ML datasets
The data are issued from a study by Zaspel et al. 75 , as used and analyzed by Pernot et al. 93 .
They contain the effective atomization energies (EAE) for the QM7b dataset 71 , for molecules up  It has been observed previously that the MP2 errors have a quasi-normal distribution, which was not the case for the SLATM-L2 dataset. Assuming a normal CDF of prediction errors through the use of a PTC1/PRED scenario, I built calibration curves to illustrate their diagnostic features.
Note that, considering the large number of points, the parametric uncertainty for the linear correction is negligible, and one is basically testing here the normality of the corrected error sets.
The plots in Fig. 11, confirm the calibration for MP2 PUs, while the notable heavy tails of the SLATM-L2 errors do not benefit from a trend correction, leading to a distorted calibration curve. The PTC1/PRED method provides us with standard PUs and confidence intervals. Their sharpness can be tested by LZV and LCP analysis, respectively. The LZV plots are presented in Fig. 12, and the LCP plots in Fig. 13.
The uniform standard PU estimated by the PTC1/PRED method seems to have severe sharpness issues, notably for MP2 and SLATM-L2 (Fig. 12). For MP2, the PU overestimates the dispersion of the errors for negative energies, while it underestimates it for positive ones. For SLATM-L2,  The PRED method enables also to generate prediction intervals for a LCP analysis. Sharpness of the PCT1/PRED scenario is evaluated in Fig. 13. One sees that despite its quasi-normal errors distribution and good calibration statistics, MP2 presents a sharpness problem at all coverage levels ( Fig. 13(a)). At the 0.5 and 0.75 levels, the LCP curves confirm the observations made through the LZV analysis (PU overestimation for the negative energy values and underestimation for the positive ones). This effect is partially lost at the 0.95 level. For SLATM-L2 ( Fig. 13(b)), the poor calibration of the PRED intervals is clearly visible, notably at the 0.5 and 0.75 levels, and sharpness problems are also detected, notably at the smaller coverage probabilities. Using a quantile-based calibration (EQ) would improve calibration but not sharpness. In the ∆-ML/CM-L1 case (Fig. 13(c)) calibration and sharpness are globally very good, but one still observes a deviation for the small negative energy values, as seen in the LZV analysis.
Considering the acceptable-to-good sharpness of the ML methods at the 95 % level, one can consider the possibility to define a uniform prediction uncertainty. In a recent article comparing MP2 (normal error distribution) with SLATM-L2 (non-normal) 93 , the authors argued that the former would probably be a wiser choice for reliable predictions. This suggestion does not stand against the present results, as MP2 is worse than SLATM-L2 at the sharpness test. In fact, it would even be reasonable to use SLATM-L2 to define an expanded uncertainty U 95 , which is not the case for MP2.

V. CONCLUSION
The use of computational chemistry in practical applications is a strong incentive to establish levels of confidence on the calculated properties. In this paper, I considered the concepts of calibration and sharpness, and the associated metrics and graphical checks, developed for probabilistic forecasters, notably in meteorology and more recently in machine learning. I adapted these methods to CC uncertainty quantification scenarios and I applied them to a series of datasets covering both embedded and a-posteriori CC-UQ methods.
The validation methods were adapted to the level of available information, notably to avoid uncontrolled hypotheses on the errors distributions. In particular, while the calibration of expanded uncertainties can be directly tested using prediction interval coverage probabilities (PICP), validation of prediction uncertainty sets was based on the analysis of z-scores variance. A methodological novelty introduced in this study is to treat sharpness as local calibration along a predictor quantity, leading to local versions of the PICP (LCP analysis) and z-scores variance (LZV analysis) tests, which proved to be very convenient and useful in our context.
The main conclusions arising from the application cases are the following: predictions and then a subsequent model to make uncertainty estimates were more calibrated than models that attempted to make value and uncertainty predictions simultaneously". We saw that a-posteriori methods can be used to ensure calibration, however, they cannot always provide sharpness.
3. The sample sizes required to test confidently either PICP values or z-scores variances often exceed the size of CC benchmark datasets. This is even more limiting when testing sharpness.
In such cases, the trend in curves of the LCP/LZV analysis provide an interesting diagnostic, even when statistical uncertainties are large.  interval (BPCI) 96 .
As gathered from the rich literature on the best method choice for the estimation of BPCIs 83,96-102 , this is a complex problem, mainly because of the discreteness of the binomial distribution. The most striking symptom is that the coverage of all BPCI estimations oscillates sharply when either ν or M varies (see 96,99 for illustrations and Fig. 15(a)

Coverage
As in the present study the main target value is p = 0.95, I compared the five selected methods for different values of the PICP ν around p and varying sample size M (Fig. 14). For all practical purposes, on sees that the Wilson and Clopper-Pearson methods produce indistinguishable intervals in these conditions. Similarly, the Agresti-Coull and continuity corrected Wilson methods (Wcc) are very closely behaved. By contrast, the Jeffreys method consistently produces shorter intervals.
All methods converge as M increases.
The coverage property of the intervals provided the continuity-corrected Wilson and Clopper-Pearson methods have been estimated by Monte-Carlo on 10 6 values of binomial proportions for a series of M values ( Fig. 15(a)). The oscillatory nature of the coverage probability is striking, and it is not evident to assess the quality of a method. CP seems to be overshooting in average, while Wcc seems more balanced, notably for M over 200. However, Wcc presents larger excursion to the shorter coverages.
A particular care has to be taken with the extreme ν = 1 value. There is an infinity of ways to realize this PICP, from barely over-scaled prediction intervals to outrageously over-scaled ones.
Using (Fig. 14), the minimal sample size required to exclude p = 0.95 from I (ν=1) 0.95 is 100. For smaller sets, one should not take the non-rejection at face value. It is necessary to examine the data to conclude: if the range of the I 95 prediction interval is much larger than the range of the test points, the positive result of the statistical test can be considered as meaningless.

Power study
The sample size required to reject confidently that a PICP value ν is compatible with the target probability p can be estimated through the power of the ν = p one-sample proportion test (probability to reject correctly this hypothesis in favor of the alternative ν = p). The power was estimated using the Clopper-Pearson confidence interval for a range of M values and several ν values, with targets p = 0.5 ( Fig. 15(b)) and p = 0.95 ( Fig. 15(c)).
From these curves, one sees that the sample size to achieve a chosen power (e.g. 0.8) increases when the |ν − p| difference gets smaller. However, a notable asymmetry is visible for p = 0. where µ k = 1/M M i=1 (x i − µ) k and µ the arithmetic mean. In absence of further information on the distribution of the variance, confidence intervals might be estimated by using the normal quantiles, leading to a test where t α,ν is the 100α % quantile of the Student's-t distribution with ν degrees of freedom. The problem here is the symmetry of the confidence interval that might be problematic for small samples.
The third option considered here is to use bootstrapping to estimate the confidence interval.
This is more computer intensive, and does not come without its own problems, as the estimation of confidence intervals from a bootstrap sample requires corrections for which several options exist 84 .
I tested here three options for the estimation of CIs: the basic, percentile and BC a methods.

CI limits estimation
To illustrate the performances of the CI estimation methods for the variance of a sample presented above, I sampled repeatedly (1000 times) from various unit-variance distribution with sample sizes between M = 50 and 1000 points. The CI limits were then averaged over the repeats and compared. Fig. 16 presents the limits for a 95 % confidence interval on Var(Z), where Z is sampled from a standard normal distribution ( Fig. 16(a)), and the effective coverage of this interval (PICP; Fig. 16(b)), for increasing sample sizes N .   The same analysis has been performed for power exponential (also known as generalized normal or error ) 107,108 distributions p(x|µ, σ, λ) = 1 2k 1/λ Γ(1 + 1/λ)σ exp −|x − µ| λ σ λ (B4) with powers λ = 1 (large kurtosis; Fig. 17) and 4 (small kurtosis; Fig. 18), scaled to unit variance.
The leptokurtic distribution (λ = 1) poses problems to all the estimation methods, and all CI Overall, the bootstrapping estimation of CI limits to the variance of z-scores seems to be a sensible approach, and I retained the BC a method with 1500 repeats as a default for this study.
One has to keep in mind that for leptokurtic z-scores distributions, the risk of false negative is slightly larger than 5 %, probably smaller than 10 % for samples sizes above N = 50.

Power study
In order to have an estimation of minimum sample sizes for reliable testing, I estimated the power of the chi-square test for a range of sample sizes and sample variances (Fig. 19). If, as for PICP testing, one targets a power of 0.8, one sees that default of variance are easier to tackle than excesses, the difference decreasing when the variances get closer to 1. With samples of 100 points, variance differences of 0.5 are testable with good confidence, while 500 points enable to go down to differences of about 0.2. Differences of 0.1 would require more than 2000 points. Of course, these data are not strictly valid for non-normal samples, but they give us a global view that variance testing is not much easier than PICP testing.

Appendix C: Atlas of calibration checks
We present calibration graphical checks for various normal z-score distributions z ∼ N (µ, σ), where µ is the mean value and σ the standard deviation. The z-scores are generated from synthetic error datasets of N = 10 4 points: