{'fraction_non_alphanumeric': 0.029709105139577356, 'fraction_numerical': 0.010827028437255413, 'mean_word_length': 4.361013986013986, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multivariate Fitting and the Error Matrix in Global Analysis of Data

When a large body of data from diverse experiments is analyzed using a theoretical model with many parameters, the standard error matrix method and the general tools for evaluating errors may become inadequate. We present an iterative method that significantly improves the reliability of the error matrix calculation. To obtain even better estimates of the uncertainties on predictions of physical observables, we also present a Lagrange multiplier method that explores the entire parameter space and avoids the linear approximations assumed in conventional error propagation calculations. These methods are illustrated by an example from the global analysis of parton distribution functions.


Introduction
The subject of this paper is a problem that arises when a large body of data from diverse experiments is analyzed according to a theoretical model that has many adjustable parameters. Consider a generic data fitting problem based on experimental measurements {D I , I = 1, . . . , N} with errors {σ I }. The data are to be compared to predictions {T I } from a theoretical model with unknown parameters {a i , i = 1, . . . , n}. A common technique for comparing data with theory is to compute the χ 2 function defined by or a generalization of that formula if correlations between the errors are known in terms of a set of correlation matrices. The physics objectives are (i) to find the best estimate of the parameters {a i } and their uncertainties, and (ii) to predict the values and uncertainties of physical quantities {X (α) , α = 1, 2, . . .} that are functions of the {a i }.
If the errors are randomly distributed, and the correlations well determined, then standard statistical methods of χ 2 minimization [1,2] apply, and established fitting tools like the CERN Library program MINUIT [3] can be employed. However, real problems are often more complex. This is particularly true in a "global analysis," where the large number of data points {D I } do not come from a uniform set of measurements, but instead consist of a collection of results from many experiments on a variety of physical processes, with diverse characteristics and errors. The difficulties are compounded if there are unquantified theoretical uncertainties, if the number of theoretical parameters n is large, or if the best parametrization cannot be uniquely defined a priori. All of these difficulties arise in the global analysis of hadronic parton distribution functions (PDFs) [4,5,6], which originally motivated this investigation. Several groups have addressed the question of estimating errors for the PDF determinations [7,8,9,10]. But the problem is clearly more general than that application.
Of the many issues that confront a global analysis, we address in this paper two, for which we have been able to significantly improve on the traditional treatment. The improvements allow a more reliable determination of the uncertainties of {a i } and {X (α) } in complex systems for which conventional methods may fail. To define these problems, we assume the system can be described by a global fitting function χ 2 global , or χ 2 for short, that characterizes the goodness-of-fit for a given set of theory parameters {a i }. This χ 2 distills all available information on the theory and on the global experimental data sets, including their errors and correlations. One finds the minimum value χ 2 0 of χ 2 , and the best estimate of the theory parameters are the values {a 0 i } that produce that minimum. The dependence of χ 2 on {a i } near the minimum provides information on the uncertainties in the {a i }. These are usually characterized by the error matrix and its inverse, the Hessian matrix H ij , where one assumes that χ 2 can be approximated by a quadratic expansion in {a i } around {a 0 i }. Once the Hessian is known, one can estimate not only the uncertainties of {a i }, but also the uncertainty in the theoretical prediction for any physical quantity X, provided the dependence of X on {a i } can be approximated by a linear expansion around {a 0 i }, and is thus characterized by its gradient at {a 0 i } (cf. Sec. 2).
The first problem we address is a technical one that is important in practice: If the uncertainties are very disparate for different directions in the n-dimensional parameter space {a i }, i.e., if the eigenvalues of H ij span many orders of magnitude, how can one calculate the matrix H ij with sufficient accuracy that reliable predictions are obtained for all directions? To solve this problem, we have developed an iterative procedure that adapts the step sizes used in the numerical calculation of the Hessian to the uncertainties in each eigenvector direction. We demonstrate the effectiveness of this procedure in our specific application, where the standard tool fails to yield reliable results.
The second problem we address concerns the reliability of estimating the uncertainty ∆X in the prediction for some physical variable X that is a function of the {a i }: How can one estimate ∆X in a way that takes into account the variation of χ 2 over the entire parameter space {a i }, without assuming the quadratic approximation to χ 2 and the linear approximation to X that are a part of the error matrix approach? We solve this problem by using Lagrange's method of the undetermined multiplier to make constrained fits that derive the dependence of χ 2 on X. Because this method is more robust, it can be used by itself or to check the reliability of the Hessian method.
Section 2 summarizes the error matrix formalism and establishes our notation. Section 3 describes the iterative method for calculating the Hessian, and demonstrates its superiority in a concrete example. Section 4 introduces the Lagrange multiplier method and compares its results with the Hessian approach to the same application. Section 5 concludes.

Error Matrix and Hessian
First we review the well-known connection between the error matrix and the Hessian matrix of second derivatives. We emphasize the eigenvector representations of those matrices, which are used extensively later in the paper.
The basic assumption of the error matrix approach is that χ 2 can be approximated by a quadratic expansion in the fit parameters {a i } near the global minimum. This assumption will be true if the variation of the theory values T I with {a i } is approximately linear near the minimum. Defining y i = a i − a 0 i as the displacement of parameter a i from its value a 0 i at the minimum, we have where the derivatives are evaluated at the minimum point y i = 0 and H ij are the elements of the Hessian matrix. 1 There are no linear terms in y i in (2), because the first derivatives of χ 2 are zero at the minimum. Being a symmetric matrix, H ij has a complete set of n orthonormal eigenvectors V v ik with eigenvalues ǫ k : These eigenvectors provide a natural basis to express arbitrary variations around the minimum: we replace {y i } by a new set of parameters {z i } defined by These parameters have the simple property that In other words, the surfaces of constant χ 2 are spheres in {z i } space, with ∆χ 2 the squared distance from the minimum.
The orthonormality of v ij can be used to invert the transformation (6): The Hessian and its inverse, which is the error matrix, are easily expressed in terms of the eigenvalues and eigenvector components: Now consider any physical quantity X that can be calculated according to the theory as a function of the parameters {a i }. The best estimate of X is the value at the minimum X 0 = X(a 0 i ). In the neighborhood of the minimum, assuming the first term of the Taylorseries expansion of X gives an adequate approximation, the deviation of X from its best estimate is given by where are the components of the z-gradient evaluated at the global minimum, i.e., at the origin in z-space.
Since χ 2 increases uniformly in all directions in z-space, the gradient vector X i gives the direction in which the physical observable X varies fastest with increasing χ 2 . The maximum deviation in X for a given increase in χ 2 is therefore obtained by the dot product of the gradient vector X i and a displacement vector Z i in the same direction with length ∆χ 2 , i.e., Z i = X i ∆χ 2 /Σ j X 2 j . For the square of the deviation, we therefore obtain the simpler formula The traditional formula for the error estimate (∆X) 2 in terms of the original coordinates {y i } can be derived by substituting X i = ∂X ∂z i = j ∂X ∂y j ∂y j ∂z i in (13) and using (6) and (10).
The result is This standard result can of course also be derived directly by minimizing χ 2 in (2) with respect to {a i }, subject to a constraint on X.
Equations (13) and (14) are equivalent if the assumptions of a linear approximation for X and a quadratic approximation for χ 2 are exact. But in practice, the numerical accuracy of the two can differ considerably if these conditions are not well met over the relevant region of parameter space. To calculate the error estimate ∆X, we prefer to use Eq. (13) using derivatives X i calculated by finite differences of X at the points z i = ± 1 2 ∆χ 2 (with z j = 0 for j = i). This is generally more accurate, because it estimates the necessary derivatives using an appropriate step size, and thus reduces the effect of higher order terms and numerical noise.
In a complex problem such as a global analysis, the region of applicability of the approximations is generally unknown beforehand. A situation of particular concern is when the various eigenvalues {ǫ i } have very different orders of magnitude-signaling that the function χ 2 varies slowly in some directions of a i space, and rapidly in others. The iterative method described in the next section is designed to deal effectively with this situation.

Iterative Procedure
In practical applications, the Hessian matrix H ij is calculated using finite differences to estimate the second derivatives in (3). A balance must be maintained in choosing the step sizes for this, since higher-order terms will contribute if the intervals are too large, while numerical noise will dominate if the intervals are too small. This noise problem may arise more often than is generally realized, since the theory values {T I } that enter the χ 2 calculation may not be the ideally smooth functions of the fit parameters that one would associate with analytic formulas. For in complex theoretical models, the {T I } may be computed from multiple integrals that have small discontinuities as functions of {a i } induced by adaptive integration methods. These numerical errors forbid the use of a very small step size in the finite difference calculations of derivatives. Furthermore, as noted above, the eigenvalues of H ij may span a wide range, so excellent accuracy is needed especially to get the smaller ones right.

The Procedure
We want to evaluate H ij by sampling the values of χ 2 in a region of parameter space where Eq. (2) is a good approximation. In principle, the parameters {z i } are the natural choice for exploring this space; but of course they are not known in advance. We therefore adopt the following iterative procedure: 1. Define a new set of coordinates {ξ i } by where u ij is an orthogonal matrix and {t i } are scale factors. In the first iteration, these are chosen as u ij = δ ij and t i = 1, so that ξ i = y i . This makes the first round of iteration similar to the usual procedure of taking derivatives with respect to a i . The iterative method is designed such that with successive iterations, u ij , t i , and ξ i converge to v ij , 1/ǫ i , and z i respectively.
2. Calculate the effective second derivative matrix Φ ij defined by using finite differences of the ξ i . The step size in ξ i is chosen to make the increase in χ 2 due to the diagonal element Φ ii ξ 2 i equal to a certain value δχ 2 . The choice of δχ 2 is determined by the particular physics application at hand. Naively, one might expect δχ 2 ≃ 1 to be the right choice. That would indeed be appropriate for a χ 2 function obeying ideal statistical requirements. But when the input to χ 2 global is imperfect, a reasonable choice of δχ 2 must be based on a physics judgement of the appropriate range of that particular χ 2 function. We therefore leave the choice of δχ 2 open in this general discussion. 2 In any case, if the final results are to be trustworthy, they must not be sensitive to that choice.
We calculate each off-diagonal second derivative by evaluating χ 2 at the four corners of the rectangle (+δ i , +δ j ), This is a modification of the technique used in MINUIT [3]. For the sake of efficiency, the MINUIT subroutine HESSE estimates off-diagonal elements using only one of those corners, together with values at (δ i , 0) and (0, δ j ) that are already known from calculating the diagonal elements of the Hessian. Our method is slower by a factor of 4, but is more accurate because it fully or partly cancels some of the contributions from higher derivatives. The first derivatives ∂χ 2 /∂ξ i are also calculated at this stage of the iteration and used to refine the estimate of the location of the minimum.

Compute the Hessian according to Φ ij ,
4. Find the normalized eigenvectors of the Hessian, as defined by Eqs. (4) and (5).
5. Replace u ij by v ij , t j by 1/ǫ j , and go back to step 1. The steps are repeated typically 10-20 times, until the changes become small and Φ ij converges to δ ij .
This iterative procedure improves the estimate of the Hessian matrix, and hence of the error matrix, because in the later iterations it calculates the Hessian based on points that sample the region where ∆χ 2 has the magnitude of physical interest.

Results from a Sample Application
As an example, we apply the iterative procedure to the application that motivated this study-the global analysis of PDFs [7]-and compare the results with those obtained from MINUIT. The experimental input for this problem consists of N = 1295 data points from 15 different experimental data sets involving four distinct physical processes. All the potential complexities mentioned earlier are present in this system. The theory is the quark parton model, based on next-to-leading order perturbative Quantum Chromodynamics (QCD). The model contains n = 16 parameters a i that characterize the quark and gluon distributions in the proton at some low momentum scale Q 0 . From a calculational point of view, the theoretical model consists of the numerical integration of an integro-differential equation and multiple convolution integrals that are evaluated mostly by adaptive algorithms. The fitting function χ 2 global in this case combines the published statistical and systematic errors of the data points in quadrature. The only correlated errors incorporated are the published overall normalization uncertainties of the individual experiments. The fitting program is the same as that used to generate the CTEQ parton distributions [4,6]. The global χ 2 minimum for this system defines the CTEQ5M1 set of PDFs, for which χ 2 0 ≈ 1200 [6]. We find that the eigenvalues {ǫ i } of the Hessian for this system range over 5-6 orders of magnitude (distributed approximately exponentially).
The value of ∆χ 2 that corresponds to a given confidence level is well defined for an ideal experiment: e.g., ∆χ 2 < 1 defines the 68% confidence region. But in a real-world global analysis, the experimental and theoretical values in Eq. (1) include systematic errors, and the uncertainties σ I include subjective estimates of those errors, so the relation between ∆χ 2 and confidence level requires further analysis. From independent detailed studies of the uncertainties [11,12], we estimate that an appropriate choice of δχ 2 for the iterative calculation is around 10 in our application, and only the region ∆χ 2 > 100 can be ruled out for the final fits.
The error matrix approach relies on a quadratic approximation to χ 2 in the neighborhood of the minimum. To test whether that approximation is valid, we plot χ 2 as a function of distance along a particular direction in {a i } space, as shown in Fig. 1. The direction chosen is a typical one-specifically it is the direction of the eigenvector with median eigenvalue. The dotted curve in Fig. 1 is the exact χ 2 and the solid curve is the quadratic approximation (2). The approximation is seen to provide a rather good description of the function. Even at points where χ 2 has increased by 50, the quadratic approximation reproduces the increase to 20% accuracy. Figure 1: Variation of χ 2 with distance along a typical direction in parameter space. The dotted curve is the exact χ 2 and the solid curve is the quadratic approximation based on the Hessian. The quadratic form is seen to be a rather good approximation over the range shown.
To correctly measure the curvature of the quadratic approximation, it is important to fit points that are displaced by an appropriate distance. This can be seen from Fig. 2, which shows the difference between the two curves in Fig. 1 in the central region. The difference displays a small cubic contribution to χ 2 . It also reveals contributions that vary erratically with a magnitude on the order of 0.03 . These fluctuations come from the noise associated with switching of intervals in the adaptive integration routines. Because the fluctuations are small, they do not affect our results in principle. But they do require care in estimating the derivatives. In particular, they would make finite-difference estimates based on small intervals extremely unreliable. The iterative method avoids this problem by choosing a suitable scale for each eigenvector direction when evaluating the Hessian. Figures 1 and 2 show the behavior of χ 2 along a single typical direction in the 16 dimensional parameter space. Fig. 3 shows a complementary test of the iterative method for all possible directions. We have chosen 1000 directions at random in {z i } space. We displace the parameters away from the minimum in each of these directions by a distance that makes ∆χ 2 = 5. We then compute the value of ∆χ 2 predicted by the quadratic approximation (2), using the Hessian calculated by the iterative method and, for comparison, by the routine HESSE within the MINUIT package. The results are displayed in Fig. 3 as histograms, with ∆χ 2 on the horizontal axis and the number of counts on the vertical axis. If χ 2 were quadratic in {a i }, then a perfect computational method would yield a delta function at ∆χ 2 = 5. Fig.  3 shows that: Figure 2: Difference between χ 2 and its quadratic approximation (2), both of which are shown in Fig. 1. A cubic contribution can be seen, along with a noticeable amount of numerical noise. The fine structure revealed here is small compared to the main variation of χ 2 itself, which rises by 20 over the region shown, as can be seen in Fig. 1.
• For the solid histogram-the result of the iterative procedure-the quadratic approximation is close to the exact result in all directions, and hence Eq. (2) is a pretty good representation of χ 2 . Quantitatively, the middle 68% of the distribution is contained in the region 5.4 ± 0.6 .
• For the dotted histogram-based on the general purpose program MINUIT-the distribution is also spread around the expected value of 5, but it is very broadly distributed. This estimate of the Hessian is therefore unsatisfactory, because we might be interested in a quantity whose gradient direction is one for which the Hessian computed by MINUIT is widely off the mark. A major source of this problem is the numerical noise visible in Fig. 2: MINUIT uses a small step size to calculate the derivatives, and gets misled by the small-scale discontinuities in χ 2 . For some directions, ∆χ 2 even becomes negative because the errors in one or more of the small eigenvalues are big enough to allow their calculated values to become negative. (Within MINUIT, this circumstance elicits a warning message, and a constant is added to all the eigenvalues, which in the context of Fig. 3 corresponds to shifting the dotted distribution to the right.) Figure 4 shows the results of a similar study, in which the 1000 random directions are chosen only from the subspace of {z i } that is spanned by the 10 directions with the largest eigenvalues ǫ i . The larger eigenvalues correspond to directions in which χ 2 rises most rapidly, or in other words, directions in which the parameters are more strongly constrained by data. Because the distance moved away from the minimum in {a i } space is smaller in this case, the quadratic approximation is generally better, so it is not surprising that the histograms are more sharply peaked than in Fig. 3. But the advantage of the iterative method remains apparent.

Comment
Information from the iteratively-improved Hessian provides a useful tool for refining the choice of functional forms used to parametrize a continuous degree of freedom in the theoretical model. Explicitly, the relevant formulas are as follows.
The length squared of the displacement vector in the space of fit parameters is while ∆χ 2 = i z 2 i by (7). Hence the directions in which the parameters are well determined (the steep directions) correspond to eigenvectors of the Hessian with large eigenvalues, while the shallow directions in which they are weakly determined correspond to small eigenvalues.
The extreme values for any particular a i are where Figure 4: Same as Fig. 3, except that the displacements are restricted to the parameter subspace spanned by the 10 steepest directions.
Equation (21) can be used to see if each parameter is appropriately well constrained. Furthermore, the individual terms in the sum show the contributions to ∆a i from the various eigenvectors, so if a parametrization leads to a poorly defined minimum because it allows too much freedom-which is indicated by a failure of the iteration to converge for the smallest eigenvalues of the Hessian-it is easy to see which of the parameters are most responsible for the too-shallow directions.

Lagrange Multiplier Method
The Hessian, via its inverse which is the error matrix, provides a general way to propagate the uncertainties of experimental and theoretical input to the fit parameters {a i }, and thence on to a measurable quantity X({a i }) by Eqs. (13) or (14). But these equations are based on assuming that χ 2 and X can be treated as quadratic and linear functions of {a i } respectively. In this section we describe a different approach, based on the mathematical method of the Lagrange undetermined multiplier, which avoids those assumptions.

The Procedure
Let X 0 be the value of X at the χ 2 minimum, which is the best estimate of X. For a fixed value of λ, called the Lagrange multiplier, one performs a new minimization with respect to the fit parameters {a i }, this time on the quantity to obtain a pair of values (χ 2 (λ), X(λ)). (The constant term −λX 0 here is not necessary, because it does not affect the minimization; but it makes the minimum value of F easier to interpret.) At this new minimum, χ 2 (λ) is the lowest possible χ 2 for the corresponding value X(λ) of the physical variable X. Thus one achieves a constrained fit in which χ 2 is minimized for a particular value of X. By repeating the minimization for many values of λ, one maps out the parametricallydefined curve (χ 2 (λ), X(λ)). Since λ is just the parameter for this curve, its value is of no particular physical significance. The relevant range for λ can be found by trial-and-error; or it can be estimated using the Hessian approximation, which predicts that λ ≈ −2 ∆χ 2 /∆X. In that approximation, F goes down by the same amount that χ 2 goes up.
One way to understand the Lagrange Multiplier method is to imagine that the quantity X is simply one of the fitting parameters, say a 1 . The variation of χ 2 with a 1 could be mapped out by minimizing χ 2 with respect to {a 2 , . . . , a n } for a sequence of values of a 1 . (That operation is indeed so useful that MINUIT provides a procedure MINOS to carry it out.) In the more general case that X is a function of all {a i }, one wants to similarly minimize χ 2 for fixed values of X. That is exactly what the Lagrange Multiplier method does, since including the undetermined multiplier term in (22) renders the {a i } independent in spite of the constraint on X.
A more phenomenological way to understand the Lagrange Multiplier method is to imagine that X has just been measured, with result X new ± σ new . To decide whether this hypothetical new measurement is consistent with the old body of data, one would add a term [(X new − X)/σ new ] 2 to χ 2 global of Eq. (1) and redo the minimization. The added contribution to χ 2 consists of a constant, a linear term in X, and a quadratic term in X. This is equivalent to Eq. (22), because a constraint on X 2 is equivalent to a constraint on X itself.
The essential feature of the Lagrange Multiplier method is that, for a given ∆χ 2 , it finds the largest range of X allowed by the global data set and the theoretical model, independent of any approximations. The full parameter space {a i } is explored in the minimization procedure, not just the immediate neighborhood of the original χ 2 minimum as in the Hessian method, and no approximations based on a small deviation from the original minimum are needed.
The only drawback to the Lagrange Multiplier method is that it can be slow computationally, since it requires a separate series of minimizations for each observable X that is of interest.

Example
We now look at an example of the Lagrange Multiplier method from our application, the uncertainty of parton distribution functions. For the physical quantity X, we consider the cross section σ W for W ± production in pp collisions at the energy √ s = 1.8 TeV of the Tevatron collider at Fermilab. We want to estimate σ W and the uncertainty on that estimate, based on the global analysis of parton distributions. The points in Fig. 5 show χ 2 global as a function of σ W B in nanobarns, where B = 0.106 is the branching ratio assumed for W → e ν. These points are obtained by the Lagrange Multiplier method using λ = 0, ±1000, ±2000, ±3000, ±4000. They are thus discrete cases of χ 2 global versus σ W , without approximations. The smooth curve in Fig. 5 is the parabola given in Eq. (14), using the Hessian computed by the iterative method and treating σ W in the linear approximation. The comparison between this curve and the discrete points from the Lagrange Multiplier calculation tests the quality of the quadratic and linear approximations and the reliability of the iterative calculation of H ij . For this application, we conclude that the improved Hessian method works very well, since the difference between the points and the curve is small, and indicates only a small cubic contribution. If the two results did not agree, the correct result would be the one given by the Lagrange Multiplier method.
To estimate ∆X, the uncertainty of X consistent with the global analysis of existing data, one needs to specify what range of ∆χ 2 global is allowed. As discussed earlier, the acceptable limit of χ 2 global depends on the nature of the original definition of this fitting function, in the context of the specific system under study. For the case of σ W , Ref. [12] estimates ∆χ 2 global ≈ 100, which translates into ∆σ W /σ W ≈ ±3% according to Fig. 5.

Conclusion
We have addressed some computational problems that arise in a global phenomenological analysis, in which a complex theory with many parameters confronts a large number of data points from diverse experiments. The traditional error-matrix analysis is based on a quadratic approximation to the function χ 2 that measures the quality of the fit, in the neighborhood of the minimum that defines the best fit. The iterative method proposed in Sec. 3 improves the calculation of the Hessian matrix which expresses that quadratic approximation, for a complex system in which general-purpose programs may fall short. The inverse of this improved version of the Hessian matrix is an improved version of the error matrix. It can be used to estimate the uncertainty of predictions using standard error matrix formulas.
Our iterative procedure for calculating the Hessian is implemented as an extension to the widely-used CERN fortran library routine MINUIT [3]. The code is available from http://www.pa.msu.edu/∼pumplin/iterate/, or it can be requested by e-mail from pumplin@pa.msu.edu. Included with this code is a test example which demonstrates that the iterative method is superior to standard MINUIT even for a χ 2 function that has no numerical noise of the type encountered in Fig. 2.
The Lagrange Multiplier method proposed in Sec. 4 calculates the uncertainty on a given physical observable directly, without going through the error matrix. It thus avoids the assumption that the theoretical quantities can be approximated by linear functions of the search parameters, which is intrinsic to the Hessian approach.
For simplicity, we have discussed only the problem of obtaining error estimates on a single quantity X. It is straightforward to generalize our methods to find the region allowed simultaneously for two or more variables by a given ∆χ 2 . For example, in the case of two variables X (1) and X (2) , the allowed region according to the Hessian method is the interior of an ellipse. The Lagrange multiplier method can be generalized for this case by adding two terms, λ 1 X (1) + λ 2 X (2) , to χ 2 .
Although the Lagrange Multiplier procedure is conceptually simple and straightforward to implement, it is slow computationally because it requires many full minimizations to map out χ 2 as a function of X, and this must be done separately for each quantity X whose error limits are of interest. In contrast, once the Hessian has been determined from the global analysis, it can be applied to any physical observable. One needs only to compute the gradient ∂X/∂a i of the observable X and substitute into Eq. (14); or better, to compute the gradient X i = ∂X/∂z i and substitute into Eq. (13). For computational efficiency, the iteratively calculated Hessian is therefore the method of choice, provided its linear approximations are sufficiently accurate. Whether or not that is the case can be determined by comparing Hessian and Lagrange Multiplier results. We use both methods in a detailed study of the uncertainties in the CTEQ5 parton distribution functions [7,11,12] that is based on the work presented here.