- Email: [email protected]
Generalized numerical error analysis with applications to geochronology and thermodynamics
Gmhrmlccl 6 Rrpmon
001a7037/87/53.00
CTCo.nmxhrm~a Acta Vol. 51. pp. 2129-2135 Jotuda Lad. 1987. Rioted in U.S.A.
+ .oa
Generalized numerical error analysis with applications to geochronology and thermodynamics* J. C. RODDICK Geological Surveyof Canada, 601 Booth Street,Ottawa, Ontario, Canada KlA 0E8
(Received January 7, 1987; accepted in revised form May 1 I, 1987) Abstrwt-A generalnumericalerrorpropagation procedure is developed to calculate the errorin a quantity derived from measurements which are subject to errors. It is an alternative to computer intensive techniques such as Monte Carlo and can be applied to quite complex analytical problems where correlation among the measurement errors and among the final errors in results are present. Previous formulations have usually ignored correlated errors. Some of the assumptions involved in error propagation can also be checked numerically. This is an advantage over analytical formulations which cannot assessthe validity of the calculated errors.Formulation as computer subroutines permits the analysis to be added to existing prcgmms. Examples from the fieldsof geochronology and thermodynamics are used to highlight the advantages and the flexibility of the method. INTRODUCTION THE ASSIGNMENT OF a realistic uncertainty in any measured quantity is important in all physical sciences. In 1939 R. T. BIRGE commented on the paucity of error analyses and outlined mathematical procedures for rectifying this. Twenty-seven years later Ku ( 1966) still considered that error propagation was neither adequately defined nor applied, and presented an expository paper on the subject. While the situation perhaps has improved in some fields of science the geological sciences often do not consider the treatment of errors in a rigorous manner. A number of texts are available on the application of statistics in geology (e.g. DAVIS, 1986) but these works usually relate to problems where the required results are overdetermined, i.e., there are more data than are required to solve the problem. For these cases least squares regression type procedures are employed to provide error estimates. In contrast, many problems contain data just sufficient to provide the required results. These types of problems, common in the fields of physics and chemistry, are also present in geology as a result of the techniques of those fields being applied to precise measurements of parameters of geochemical interest. Error propagation is an important technique that can provide estimates of the uncertainties of these parameters. In computing any quantity by combining a number of measurements in some mathematical relation, the uncertainties in the measurements are passed on or propagated into the final quantity. Calculations using the principles of error propagation to obtain the final error in a parameter are often very tedious to derive, as part of the procedure involves determining partial derivatives for each measurement variable. Additional complications arise if there are correlations among the measurements. Despite these difficulties, error prop
* Geological Survey of Canada Contribution No. 44286.
agation calculations have been used in a number of geological applications (e.g. DALRYMPLE and LANPHERE, 197 1; LUDWIG, 1980) but the algebraic expressions derived apply to specific formulations of the problems and often contain approximations. In addition, any change in a formula requires a new error equation derivation. In order to propagate errors in a general way ANDERSON (1976) suggested the use of a Monte Carlo technique. One disadvantage of this ap preach is its intensive use of computer time as the technique involves repeated calculations with randomly-varied errors in the measurements REES ( 1984) suggested that numerical differentiation procedures could obtain more accurate results without extensive computer calculations. A more important limitation of both the Monte Carlo technique, as formulated by Anderson, and the numerical procedures suggested by Rees, is that correlation of errors is not accommodated. Correlation both among the measured variables and among the final results can be extremely important and is more common than is usually acknowledged. Recently, POWELL and HOLLAND (1985) and HOLLANDand POWELL( 1985) have applied numerical error propagation with correlation to thermodynamic problems but do not describe the error analysis in detail. The presentation that follows attempts to rectify most of the problems of error analysis discussed above. The basic equations of error propagation are rewritten in a form which can be easily solved by numerical means and applied to a variety of situations. Correlation is accounted for in the numerical calculations. Computer program subroutines have been developed which can be incorporated into most programs. They are self contained and do not require additional routines such as normally-distributed random number generators as for Monte Carlo methods. With these subroutines it is easy to assess the effects of errors on some or all the measurement variables in a problema particularly useful technique in designing new experiments. It is often not appreciated that error prop agation is an approximate method with a number of
2129
J. C. Roddick
2130
assumptions that should be fulfilled if the results are to he valid. These assumptions are discussed and procedures to test the validity of some of them are incorporated in computer subroutines. Finally several ap plications are presented which highlight the value and the problems of numericaf error analysis. Emphasis is placed on the practical application of the technique. ANALYTICAL
the standard deviations of the means in percent as Ci (termed % error), where Ci = IOOUJXi.
(5)
Then the variance of the mean of parameter G (in 4%’units) is:
FORMULATION AND NUMERICAL CALCULATION
Many problems in physical science involve the determination of some pururnefer G that is a function of several measured variabtes which are subject to uncertainty in their values. The relationship of the variables to G can be expressed as a known mathemati~ function G = _&x1,x2, - * * x.) vahd over a cettain range of the x1 to x. variables. The variance of G can be approximated by using the geneml equation for the propagation of errors in the form (PUGH and Wmsr,ow, 1969):
4i!&~+($q+.
Provided the partial derivatives can be determined this equation can be used to derive an analytical expression for the propagation of errors to a given parameter. A slight rearrangement of the terms and replacement of the partial differentials with the ratios of small increments @f; ax*)emphasizes a numerica solution. Equation (6) can be expressed as:
. . +2($9(j3m +2f
( I(221 dxg
ax,
u13+ * - *
(1)
or in summation notation (HAHN and SHAPIRO,1967):
i
correlationco&icienr: (3) The division of the covariance by the associated standard deviations results in a value of pij that is independent of the units of measure. The de~nitions of the variance and covariance of Xi and Xj limit the correIation coefficient to a range of - 1 to +I with p0 = 1 for perfect positive correlation of xi and x,, pij = -1 for perfect negative correlation and pii = 0 (i.e., q = 0) for uncorrelated variables (PUGH and WINSLOW, 1969). The value of pil is estimated from m pairs of observations; i.e., for variables u (‘Xi) and D (=X1)with mean values iiand C respectively: m
V/f The terms of the form represent ratios of smalt fracaxj/xt tional changes in the function resulting from small ti-actional changes in each variabIe hdding the remaining variables constant, and for sufficiently small changes closely approximate first partial derivatives. These terms may be determined numerically by sequentially deviating each variable from its mean value by a small fractional amount, while keeping all other variables at their mean values, and determining the cormsponding change in the functionj: In practice a deviation of 0.01% from the mean value of the variable is adequate to approximate first partial derivatives of a range of functions. Iffis a linear function of a variable this calculation of the ftrst partial derivative is exact As expfained in the next section, ftmctions should be a~~xirna~iy linear in the region of evaluation for the error propagation equation to be valid, thus the use of more complex interpolation formulae such as Sterling’s method to calculate fitst derivatives is unneomsary. Because of their association with the % error of each variable these terms represent an error mugn@cution or EMAG factor related to that particular variable. Thus a tabulation of these
EMAG factorscan kc used to assessthe potential contribution of each variable to the final error in the parameter G. Alternatively if the errors are known or estimated, variance contributions can be calculated as the square of the product of the % error and the EMAG factor. Each pair of variables have a covariancc, expressed as a correlation coefficient, hut in practice many of the variables are not correlated and py = 0. Thus for the double summation in Eon. (7) only variables with significant correlation coefficients (known or estimated) need be evaluated. Having obtained the necessary variance and covariance terms the summations are carried out to calculate a total variance for the parameter. The vahres of other parameters which are also fun~o~ly related to the variables in the G function (i.e., B = h(x, , x2, * - * x,,) may be required. As G and B are functions of the same variables they are correlated and have a covariance given by:
Estimated pii = ruv=
Equation (2) may be written in a form that is independent of the units of measure by substituting Eon. (3) and expressing
It is more convenient to express thii covatiance as a correlation coefficient for parameters G and B. If we also express the errors in % using Eqn. (5) and express the partial differentials as ratios of increments as was done in Eqn. (7) then the correlation coefficient is:
Numerical error analysis applications
2131
surements. Also, if a normal distribution is an adequate approximation most of the higher moments are zero. A detailed study of the significance of the omitted higher order terms concluded that the approximation “turns out to be better than expected” (TUICEY, 1967, This equation is of a similar form to Eon. (7) and the EMAG factors,variance and covariance components may be evaluated as outlined above. To evaluate any number of parametersand their correlation coefficients it is necessaq to repeat the complete calculation of the parameters n + I times where n is the number of variables. In each cycle the EMAG factor for each variable is determined. These calculations can be expressed and solved using matrix notation (CLIFFORD,1973) but the presentation given here clearly outlines the approach to the numerical calculations. A general algorithm of this procedure has been written in Hewlett Packard (200/300 series) Basic and allows the evaluation of up to 8 parameters having up to 25 variables, though this may be increased if required. The algorithm can be added to existing programs with a minimum of modifications.
The propagation of errors, fluctuations and tolerance. Unpubl. tech. reports 10, I 1, 12. Princeton University
Press). Ku ( 1966) suggests in practice, that the higher terms in the error propagation equation can be neglected if the second and higher order partial derivatives, when evaluated at the means, are small and that the standard deviations of the variables are also small. For small standard deviations Ku suggests values of 0.1% or less, though examples given below suggest standard deviations of several percent are acceptable in most cases. For the requirement of small higher partial derivatives Ku does not suggest any values but the use of numerical calculations permits further tests to assess the significance of the higher order terms, and also large variable ACCURACY OF THE ERROR PROPAGATION errors, to the accuracy of the calculated parameter erEQUATION rors. The mathematical derivation of the error propagaThese tests involve comparing two methods of caltion equation involves a number of assumptions and culating EMAG factors. The previously described calapproximations, which affect the validity of the final culation of the EMAG factors involves successive .O1% error estimate in calculated parameters, and of which deviations in each variable. To obtain a variance conthe user should be aware. In addition to the numerical tribution the squared product of a variable’s EMAG approximations to determine the &t partial derivatives factor and estimated % error is calculated (Eqn. 7). discussed in the last section the other main assumptions The estimated % errors of the variables may be as large are: as several percent in contrast to the small deviation used for the EMAG calculation. An alternate value for (1) The mean values of the variables and calculated the EMAG factor can be obtained by replacing the parameters am derived from approximately normally .Ol% deviation with the estimated SJ error of the distributed populations. variable. These new alternate EMAG factors can then (2) The error propagation equation is an approxibe compared with the original factors to assess if the mation from a Taylor series expansion about the point function is approximately linear (i.e., constant EMAG at which each of the component variables take on their factors or partial derivatives) at the point of evaluation mean value. Terms in the series higher than the first and over the range of the error limits. If an EMAG partial derivatives are considered negligible and thus factor is not constant then the linear approximation is omitted from the analysis. inadequate and higher order terms may be significant. This procedure can be carried out for all variables with Empirical observations have shown that the normal distribution is a good representation for many popua check on the linearity of the calculations to some lations in physical science. Even if it is only an ap arbitrary level. The computer subroutine flags differences greaterthan 5%in the two methods of calculation proximation the central limit theorem states that the sampling distribution of the mean of n independent of EMAG factors and thus highlights potential probobservations of any distribution with finite mean and lems. Analytical formulations of error propagation are variance approaches a normal distribution as n ap incapable of assessing these non-linearity effects. Calculations with alternate EMAG factors can prc+ proaches infinity. For a parameter dependent on four or more variables from different populations the pa- vide some estimate of the uncertainty of calculated parameter errors. The algorithm first calculates alterrameter mean tends to be normally distributed even for small n (EISENHART, 1963; HAHN and SHAPIRO, nate EMAG factors by increasing the mean of each variable by its estimated error limit. Variance contri1967). butions are then calculated and summed to give one The omission of the higher order terms in the Taylor error limit on the parameter. The procedure is then series is usually a valid approximation for most functions. In the case of linear functions (e.g. G = Ui + bX_) repeated but with alternate EMAG factors calculated by decreasing the mean of each variable by its estimated the higher terms vanish and Eon. (2) is exact. inclusion of the higher order terms is usually not very practical error limit. This provides another limit ofthe parameter error. If the function is linear in all variables the two as they contain the higher moments of a distribution (the variance is the second moment; HAHN and SHA- limits will be equal to the original calculation and the PIRO, 1967) which are not available for many meaerror calculation is exact. Non-linear behavior results
J. C. Roddick
2132
in a range in the error limits, however, these limits are only approximations as the method of calculation with alternate EMAG factors presumes a positive correlation in the determination of the Eh4AG factors. This arises from the fact that the calculation sums the variance contributions derived from alternate EMAG factors calculated by changing (increasing or decreasing) each variable in the same way. If several variables cause the function to behave in a non-linear fashion the calculated parameter errors will be poor approximations as evaluation also holds all other variables at their mean values. The net effect of simultaneous random variation of several variables could be to enhance or cancel the non-linear effects. If the linear approximation is inadequate the correct variance would be obtained only if a parameter were a Bmction of a single variable. For situations where the linear appro~mation is inadequate Monte Carlo simulation may be a better approach to the error evaluation though it would be easiest with functions with uncorrelated variables. EXAMPLES OF ERROR ANALYSIS In this section three examples of error propagation
analysis from thermodynamics and geochronology are given. The examples range in increasing order of complexity from the determination of parameters of linear equations to a relatively complex example where the equations are not linear, some variables are correlated and a correlation coefficient for two parameters is required. In applying the computer routines it is necessary to define a subroutine containing the calculations of interest as they would normally be computed with measured variables with known errors and any correlations among variables. The particular variables which have error estimates are defined by equating them to a working array, while results to be calculated are defined in another array. Correlated variables are specified as array element pairs. For any given set of variables, errors, and correlations a call to the error calculation routine will produce the calculated results with their associated errors and correlation coefficients between results for that particular data set. A table of indi~d~ EMAG factors or variance contributions for each variable may be listed, and any significant non-liner EMAG calculations noted, along with the uncertainties in the errors due to non-linearity. The examples which follow display typical tables of error analysis.
2 chlorite + 8 kyanite -t 1 I quartz = 5 con&rite + 5.5 Hz0
the equation for determining the Gibbs free energy of cordierite is (BIRD and ANDERSON,1973): 5 AGt,, = AG,, + A&J T- 25) -Al’,@‘-l)-5SCPHrO
(11)
where: AGf = standard
P= T= G*H20 = AT, =
A&$ =
molar Gibbs free energy of formation from the elements at 298”K, 1 bar (standard state); r = reactants, equilibrium pressure, equilibrium temperature, a form of Gibbs free energy of I&O at P and T (USHER and ZEN, 197 I), molar volume change of the solids in a reaction, third law molar entropy change of the solids and for Eqn. (10) is:
asf,~=5~f.~-2hs~chl-s~f.Ly-
ll~fmz.
(12)
An analogous equation may be written for AV,. Data for this reaction, with 2 standard errors of the means are listed in BIRD and ANDERSON( 1973). Table I details a numerical error analysis of these data at 7000 bars and 65O”C, the same P, T conditions as used by ANDERSON(1976) in his Monte Carlo analysis The table presents all 14 variables with the error prop agation of each variable expressed both in terms of the EMAG factor and the variance. EMAG factors may be positive or negative and reflect the resultant fractional change in AG- for positive change of a variable. Many of the EMAG factors are small (~0.1) and when coupled with the errors on the variables contribute essentially nothing to the variance of AG,_. This ability to separate error components highlights the influence of each variable on the total variance of a parameter. Equation (11) is very nearly linear in this case (7000 bars, 65O’C) as the P and T variables contribute little to the total variance in Table 1. This means that the error propagation calculation can be regarded as es-
VARIABLE NAME
PARAMETERI VALUE: 96 ERRORc
-ai ,219
% ERROR’
EMAct
96’ VAWANCE
-.02?? a022 -.0164 .I760 .w73 -.0626 .ZlSS
.OOLI .OoOO .oocQ .OZU .0018 .0143 .0019
GZ,O
1.34 I.&3 :f’
6=ci11 A=b 2”s
.OWZ 1.91 .2
TmP
Free energy of formation of cordierite Tlte free energies of formation of mineral phases are required for the understanding of chemical reactions in rock-forming systems. BIRD and ANDERSON( 1973) used hydrothermal equilibrium data to calculate the free energy of cordierite and later ANDERSON(1976) applied a Monte Carlo technique to assess the errors in the calculation. For the reaction:
(10)
A&W
.221
.I223
.0007
igz
.7la :YF .0336 .I61 .a036
-.0287 - .0577 -.OOa .OL87 -.Oo% -.OO@O
.0021 .OWO .OOOfl .OOOO .OOOO .cmO
6 Vch, AVcw 6Vky AVsxt~
’ = Errtmattwostmdudemm t
= Error rrugificaianof
nra
in cakulati~
ACf,m
Numerical
sentially exact. The calculated value (-2089.15 Kcall mole) is the same as that obtained by BIRD and ANDERSON (1973) while the error of .219% (2 sigma) is slightly larger than the value of .206’S obtained by hand calculations (BIRD and ANDERSON, 1973) and the value of .2 12% by the Monte Carlo method (ANDERSON, 1976). The numerical calculation has thus obtained an exact result in 15 calculation cycles where the Monte Carlo employed 862 cycles. The chlorite variables AGfchl and A&,,,, dominate the variance in the final error in AGr.,,. Table 1 thus shows that improved measurements in these variables would contribute most to an improvement in the free energy calculation of cordierite. Correlation in thermodynamic calculations
ANDERSON(1977) discussed the general characteristics of error propagation in thermochemical calculations of mineral equilibria such as Eqn. ( 11). He provided a clear example, with the calculation of enthalpies of sanidine and muscovite, of how data are often interdependent with highly correlated relations that unfortunately are usually ignored when applying the data to mineral equilibria calculations. In solution calorimetry the calculation of the enthalpy of formation of a single compound (Mr) involves the measurement of a number of heats of reaction. For muscovite 12 independent measurements are required while for sanidine 11 measurements of heats of reaction plus its heat of solution are required. Eleven of these measurements are common to the A?& of both minerals. A thermochemical calculation may require the enthalpy of the reaction muscovite = sanidine + corundum + water.
2133
error analysis applications
(13)
In this case the algebraic sum of enthalpies from the four components are required. In calculating the uncertainty in this enthalpy of reaction the covariance in the errors of muscovite and sanidine, as a result of the common source of their enthalpies, should be considered Anderson’s procedure was to calculate the variance of (Ai&*, - A&,,,) from the original 13 independent reactions. He emphasized that, while providing the correct variance for these two minerals in this reaction, this procedure cannot be used generally. In calculating the difference of AHr six of the 13 reactions cancel out and do not contribute to the variance of the difference. In contrast if the sum of enthalpies of muscovite and sanidine were required the individual heats of reaction do not cancel and thus result in a greater variance. Numerical error propagation, however, can provide a general solution to this problem. Using the 13 independent reactions and the linear summation equations for AHCW and AHcmur given in Appendix 2 of ANDERSON(1977) the errors and the correlation coefficient for these enthalpies may be calculated.ThevaluesareA&,,,=-1421.23+-1.192 Kcal and AHr_ = -946.453 + -1.098 (2 sigma errors), the Same as determined by Anderson. In addition
the correlation coefficient between these two enthalpy errors is .9 19 and confirms that the errors are highly correlated. The determination of the correlation coefficient provides a value which can be used in other calculations of mineral reactions regardless of whether and AHr,, are added or subtracted. For the Mr.,,, example, for Eqn. ( 13) the errors and correlation coefficient can be included in the numerical error analysis to calculate the enthalpy of reaction. This can be done by specifying which variables are correlated and the value of the correlation coefficient. The resulting covariance will then be propagated into the calculated error on AZ& This example shows how numerical error propagation may be used in several stages. The resulting parameters with their errors and correlation from one calculation may become the variables in a subsequent error analysis. U/Pb ratios for concordia diagrams
In U-Pb geochronology in order to calculate the age of mineral phases that contain U but very little Pb when they are formed (termed common Pb) the amounts of mdiogenic 207Pb(20’Pb*) and x”%%(%Pb*) in a sample are required. Often the ratio of m7Pb*/ 206Pb* provides an estimate of the age, but in most cases there is a loss of radiogenic Pb subsequent to crystallization of a mineral. In these cases it is necessary to calculate ratios of zo7Pb*/238U and 2osPb*/23sU of several sub-samples of a mineral and assess their significance on a concordia diagram (see FAURE (1977) for principles). Thus these U/PC, ratios and their errors for several sub-samples are required for subsequent assessment, often with two-error regression treatments (e.g. YORK, 1969). In general the total Pb in a mineral analysis consists of three components: radiogenic Pb, common Pb and blank Pb added during sample processing. If samples are analyzed by the isotope dilution technique with a mixed isotopic spike containing (for simplicity in this example) pure *05Pb and pure 235U, then the concentrations of the Pb and U isotopes are (abbreviating the isotopes to their masses): 206Pb* = (206/205),(205,)
- 206u - (206/204),(204),,
207Pb* = (207/205),(205,)
( 14)
- (207’204)bi(206u) (206/204), - (207/204),(204,)
(15)
238U = (238/235),/ [ 1 - (238/235),/Ku](235,)
- 238,,,
(16)
where the subscripts sp, m, bl, cm refer to spike, measured, blank and common Pb or U isotopic compositions respectively. Ku = 137.88 = (238/235) in the sample. 205,, 238, = the amounts of spike *05Pband 238Urespectively in the mixture.
J. C. Roddick
2134
PARAMETER: VALUE: % ERROR: VARIABLE NAME
% ERROR
2071206,
.2
3 4 5 6 7
;;;g& 238l235,,, pb Frac. 2061205sp ;3;g$P
.3 .5 .I
9 LO II
Fo Blank 206/2wbl 207/20$,~ 206/2Wc,,, 207/2%,
8
f:
30 .03 30 .05 30 3 .7 2 .4
2imt231 .2019 .a25
207*1235 5.598 1.21
%’ VARIANCE .0640 .I439 .0148 .0103 .0030 .woo .0063 .0367 .2251 .0181 2 .I579 Z
207*/206* .2011 .783
COR
COEFf .7660
CONTRIBUTIONS 2 .4940 .4404 .0103 .0484 2 .0063 .0367 .2136 .0021 .0330 .A76
.O#O .I046 -2936 .0:72 .OOw Z . of2 .0326 .0330 .I579 .I876
Z .2666 .0808 .0103 .0121 .:063 .0367 .2192 -.0062-
-.oow :
COR. COEF. v ( NO’S) 10, 11 12, 13
% ERRORS 3, .7 2, .4
.; .7
TOTAL VARIANCE:
t Z +, +
= = = =
.6802
206*/238 ERROR -.I0 ERROR .I1
IN PARAMETER
= (204/205),205,
.0171 .I205
.6137
.7636
ERROR3 FOR a/- 1NPUT ERROR3
207’1235 -.I4 .I4
207+1206* -.37 .39
206bl = the amount of mPb in a separate blank analysis processed at the same time as the sample. 204, = the amount of 204Pbinitially in the sample, where: 204,
I.4608
-.wM) -.2409
Correlation Coefficient for first tvm parameters EMAC factor and variance wle zero fa this variable EMAG factor subject to > 5% variation in calculation far +/- Input Errors EhtAG factor increase (- = EMAC decrease) for input Erra incrure
96 UNCERTAINTY
+ -
%’ COVARIANCE -.Mw)6 -.0118 . 0000 0.0000
- 206u/(206/204)u.
( 17)
Equations ( 14 to 17) and the value of Ku may be combined to obtain the required U/Pb ratios. Table 2 presents a numerical error analysis of the U-Pb ratios and m7Pb*/206Pb* using all the variables required in Eqns. (14, 15 and 16). In addition 3 other variables are listed: the (238/235), of U is normalized for mass fractionation’ using a double 233-235 spike and so (233/235), is listed; measured Pb ratios, subject to isotopic fractionation, are normalized thus Pb fractionation (Pb. Frac.) is a variable; and the 2osPb spike contains about 5% 206Pb so the spike 206Pb/205Pbis listed though itscontribution is insignificant in this example. Only variance contributions are listed in the table but a distinction is made between very small variances (.OOOO)and those that are zero because EMAG factors are zero (Z). LUDWIG (1980) presented an analytical solution for error propagation in U-Pb ratios similar to that given here and a worked example. The data in Table 2 are the same as Ludwig’s with about 10 ng of sample Pb and 1 ng of blank Pb. The parameter error estimates and correlation coefficient are in good agreement with Ludwig’s with slight differences related to different
’ The (238/235), is arbitrarily nomalized to a fixed value of 233/235 during U analysis. Consequently error in the analysis is expressed only by error in (238/235),.
starting points; Ludwig assumes that the concentrations of 206, and U were previously determined with known errors. Here the error analysis used mainly isotopic ratios derived directly from mass spectrometer analyses. A complete analytical solution to Eqn. (6) using Eqns. ( 14 to 17) and mass spectrometer data would be tedious to derive and subject to algebraic errors while the numerical procedures can overcome these difficulties. In this example correlations between two pairs of variables are included. The common Pb isotopic compositions and the blank isotopic compositions have been given the same correlation coefficient of 0.7 as in Ludwig’s example. These correlations have little intluence on the variance of the U-W ratios but &nihcantly reduce the variance on the 207*/206* as is reflected by the negative covariance components. This is important in the case of a sample which does not have significant Pb loss since its age and error may be calculated directly from the 207*/206* ratio and its error. The uncertainties in spike U and Pb concentrations have not been included in the errOranalysis for two reasons. Because a mixed spike is employed only its U-Pb ratio is required to determine sample U-Pb ratios. While the quantity of spike used is required to make blank corrections, the contribution of uncertainty in this factor (weight X concentration) to the errors in UPb ratios is not significant. Secondly, different samples regressed on a coneordia plot are analyzed with the same mixed spike, so should the calibration of the spike U-l% ratio be biased (within calibration errors) all samples will have the same bias. Thus the uncertainty in the spike U-Pb ratio should be included in the regression procedure and then, only for comparison of ages with data from other laboratories using a different U-Pb spike.
2135
Numerical error analysis applications
These calculations of U/W are derived from equations that are not linear and so significant approximations might be expected especially in the case of error limits as high as +30% as listed for the blanks and Pb fractionation. Using the procedure outlined in the previous section on accuracy all EMAG factors are linear to better than 3% except for the (206/204A, variance contribution (>5%) to the correlation coefficient (Table 2). However, this contribution to the total variance is small while the uncertainty in the propagated error estimates of the three parameters ranges from
use of these parameters in geological and geochemical models. Acknowledgements-M. J. Dalwitz is thanked for unknowingly enlightening me, more than 10 years ago, on the numerical approach to error analysis. Critical comments from Terry Gordon and Rob Berman greatly improved the text and ensured that erroneous ideas did not creep into the mathematical development. An anonymous reviewer is thanked for correcting some misconceptions in the U-Pb example. Editorial handling:D. M. Shaw
only 0.10% to .39% (Table 2). Thus non-linearity is very minor and the error calculations are considered
to be valid. These examples demonstrate the general application of the technique. More complex problems involving additional parameters and correlations are treated in a similar way. While non-linearity in these examples did not produce problems in the error analyses other applications may have significant non-linearity. An example of this sort of problem in 40Ar-39Ardating will be given elsewhere (RODDICK, in prep.). CONCLUSIONS The application of rigorous error analysis to many analytical problems in geology is seldom a usual procedure. There are several reasons for this but the main one is a lack of a technique that can be applied to a range of problems. Numerical error propagation can provide the required general technique. It can be ap plied to complex problems where analytical procedures are difficult and is faster and more informative than a Monte Carlo method. In addition, the procedure can also test whether some of the assumptions involved in the error propagation equation have been fulfilled and thus if the error estimates are- valid. The equations required for numerical error propagation, including correlations are given and the method of calculation discussed. Examples from thermodynamics and geochronology outline the procedures and the advantages of the numerical technique. They also highlight the importance of correlation among errors of variables in the problem and among errors in the final results that are calculated. An algorithm written as subroutines in HP Basic was employed in the calculations and listings of the program for each example are available upon request. Wider application of rigorous error analysis to geological problems will provide a better appreciation of the uncertainties on parameters derived from complex problems. This will permit a more realistic assessment of the conclusions to be drawn from the
RJCFJSRJZNCE!S ANDERSONG. M. (1976) Error propagation by the Monte
Carlo method in geochemical calculations. Geochim. Cosmochim. Acta 48,2309-23 I 1. ANDERSONG. M. (1977) The accuracy and precision of Calculated mineral dehydration equilibria In Thermodynamics in Geology,pp. I 15-135. D. Reidel, Dordrecht-Holland. BIRDG. W. and ANDERSONG. M. (1973) The free energy of formation of magnesian cordietite and phlogopite. Amer. J. Sci. 273, 84-9 I. BIRGER. T. (1939) The propagation of errors. Amer. Phys.
Teacher 7,351-357. CLFFXXDA. A. (1973) Multivariate ErrorAnalysis. J. Wiley and Sons, 112~~. DALRYMPLEG. B. and LANPHERE M. A. (1971) *Ar/=Ar
technique of K-Ar dating: a comparison with the convcntional technique. Earth Planet Sci. L&t. 12,300-308. DAVISJ.
C. (1986) Statistics and Data Analysis in Geology.
J. Wiley and Sons, 55Opp. FAUREG. (1977) Principles of Isotope Geology.J. Wiley and Son%464pp. E~SENHART C. (1963) Realistic evaluation of the precision and accumcy of instrument calibration systems. J. Res. NBS
C. Eng. Instrum. 67C, 161-188. FISHERJ. R. and ZEN E-AN (197 I) Thermochemical calculations from hydrothermal phase equilibrium data and the bee energy of formation of HrO. Amer. J. Sci. 270,297-
314. HAHN G. J. and SHAPIROS. S. (1967)
Statistical Models in Engineering.J. Wiley and Sons, 355pp.
HOLLANDT. J. B. and POWELL R. (1985) An internally consistent thermodynamic dataset with uncertainties and carrelations: 2. Data and results. J. Meta. Geol. 3.343-370. KU H. H. (1966) Notes on the use of propagation of error formulas J. Res. NBS C. Eng. Instrum. 7oC, 263-273. LUDWIGK. R. (1980) Calculation of uncertainties of U-Ph isotope data. Earth Planet. Sci. Lett. 46.2 12-220. POWELL R. and HOLUND T. J. B. (1985) An internally consistent thermodynamic dataset with uncertainties and correlations: I. Methods and a worked example. J. Meta. Geol.
3,327-342. PUGH E. M. and WINSLOWG. H. (1969) The analysis of physical measurements. Addison-Wesley, 246~~. REEKC. E. (1984) Error propagation calculations. Geochim.
Cosmochim. Acta 48,2309-23 I I. YORK D. (1969) Least squares fitting of a straight line with
correlated errors. Earth Planet. Sci. Lett. 5, 320-324.
001a7037/87/53.00
CTCo.nmxhrm~a Acta Vol. 51. pp. 2129-2135 Jotuda Lad. 1987. Rioted in U.S.A.
+ .oa
Generalized numerical error analysis with applications to geochronology and thermodynamics* J. C. RODDICK Geological Surveyof Canada, 601 Booth Street,Ottawa, Ontario, Canada KlA 0E8
(Received January 7, 1987; accepted in revised form May 1 I, 1987) Abstrwt-A generalnumericalerrorpropagation procedure is developed to calculate the errorin a quantity derived from measurements which are subject to errors. It is an alternative to computer intensive techniques such as Monte Carlo and can be applied to quite complex analytical problems where correlation among the measurement errors and among the final errors in results are present. Previous formulations have usually ignored correlated errors. Some of the assumptions involved in error propagation can also be checked numerically. This is an advantage over analytical formulations which cannot assessthe validity of the calculated errors.Formulation as computer subroutines permits the analysis to be added to existing prcgmms. Examples from the fieldsof geochronology and thermodynamics are used to highlight the advantages and the flexibility of the method. INTRODUCTION THE ASSIGNMENT OF a realistic uncertainty in any measured quantity is important in all physical sciences. In 1939 R. T. BIRGE commented on the paucity of error analyses and outlined mathematical procedures for rectifying this. Twenty-seven years later Ku ( 1966) still considered that error propagation was neither adequately defined nor applied, and presented an expository paper on the subject. While the situation perhaps has improved in some fields of science the geological sciences often do not consider the treatment of errors in a rigorous manner. A number of texts are available on the application of statistics in geology (e.g. DAVIS, 1986) but these works usually relate to problems where the required results are overdetermined, i.e., there are more data than are required to solve the problem. For these cases least squares regression type procedures are employed to provide error estimates. In contrast, many problems contain data just sufficient to provide the required results. These types of problems, common in the fields of physics and chemistry, are also present in geology as a result of the techniques of those fields being applied to precise measurements of parameters of geochemical interest. Error propagation is an important technique that can provide estimates of the uncertainties of these parameters. In computing any quantity by combining a number of measurements in some mathematical relation, the uncertainties in the measurements are passed on or propagated into the final quantity. Calculations using the principles of error propagation to obtain the final error in a parameter are often very tedious to derive, as part of the procedure involves determining partial derivatives for each measurement variable. Additional complications arise if there are correlations among the measurements. Despite these difficulties, error prop
* Geological Survey of Canada Contribution No. 44286.
agation calculations have been used in a number of geological applications (e.g. DALRYMPLE and LANPHERE, 197 1; LUDWIG, 1980) but the algebraic expressions derived apply to specific formulations of the problems and often contain approximations. In addition, any change in a formula requires a new error equation derivation. In order to propagate errors in a general way ANDERSON (1976) suggested the use of a Monte Carlo technique. One disadvantage of this ap preach is its intensive use of computer time as the technique involves repeated calculations with randomly-varied errors in the measurements REES ( 1984) suggested that numerical differentiation procedures could obtain more accurate results without extensive computer calculations. A more important limitation of both the Monte Carlo technique, as formulated by Anderson, and the numerical procedures suggested by Rees, is that correlation of errors is not accommodated. Correlation both among the measured variables and among the final results can be extremely important and is more common than is usually acknowledged. Recently, POWELL and HOLLAND (1985) and HOLLANDand POWELL( 1985) have applied numerical error propagation with correlation to thermodynamic problems but do not describe the error analysis in detail. The presentation that follows attempts to rectify most of the problems of error analysis discussed above. The basic equations of error propagation are rewritten in a form which can be easily solved by numerical means and applied to a variety of situations. Correlation is accounted for in the numerical calculations. Computer program subroutines have been developed which can be incorporated into most programs. They are self contained and do not require additional routines such as normally-distributed random number generators as for Monte Carlo methods. With these subroutines it is easy to assess the effects of errors on some or all the measurement variables in a problema particularly useful technique in designing new experiments. It is often not appreciated that error prop agation is an approximate method with a number of
2129
J. C. Roddick
2130
assumptions that should be fulfilled if the results are to he valid. These assumptions are discussed and procedures to test the validity of some of them are incorporated in computer subroutines. Finally several ap plications are presented which highlight the value and the problems of numericaf error analysis. Emphasis is placed on the practical application of the technique. ANALYTICAL
the standard deviations of the means in percent as Ci (termed % error), where Ci = IOOUJXi.
(5)
Then the variance of the mean of parameter G (in 4%’units) is:
FORMULATION AND NUMERICAL CALCULATION
Many problems in physical science involve the determination of some pururnefer G that is a function of several measured variabtes which are subject to uncertainty in their values. The relationship of the variables to G can be expressed as a known mathemati~ function G = _&x1,x2, - * * x.) vahd over a cettain range of the x1 to x. variables. The variance of G can be approximated by using the geneml equation for the propagation of errors in the form (PUGH and Wmsr,ow, 1969):
4i!&~+($q+.
Provided the partial derivatives can be determined this equation can be used to derive an analytical expression for the propagation of errors to a given parameter. A slight rearrangement of the terms and replacement of the partial differentials with the ratios of small increments @f; ax*)emphasizes a numerica solution. Equation (6) can be expressed as:
. . +2($9(j3m +2f
( I(221 dxg
ax,
u13+ * - *
(1)
or in summation notation (HAHN and SHAPIRO,1967):
i
correlationco&icienr: (3) The division of the covariance by the associated standard deviations results in a value of pij that is independent of the units of measure. The de~nitions of the variance and covariance of Xi and Xj limit the correIation coefficient to a range of - 1 to +I with p0 = 1 for perfect positive correlation of xi and x,, pij = -1 for perfect negative correlation and pii = 0 (i.e., q = 0) for uncorrelated variables (PUGH and WINSLOW, 1969). The value of pil is estimated from m pairs of observations; i.e., for variables u (‘Xi) and D (=X1)with mean values iiand C respectively: m
V/f The terms of the form represent ratios of smalt fracaxj/xt tional changes in the function resulting from small ti-actional changes in each variabIe hdding the remaining variables constant, and for sufficiently small changes closely approximate first partial derivatives. These terms may be determined numerically by sequentially deviating each variable from its mean value by a small fractional amount, while keeping all other variables at their mean values, and determining the cormsponding change in the functionj: In practice a deviation of 0.01% from the mean value of the variable is adequate to approximate first partial derivatives of a range of functions. Iffis a linear function of a variable this calculation of the ftrst partial derivative is exact As expfained in the next section, ftmctions should be a~~xirna~iy linear in the region of evaluation for the error propagation equation to be valid, thus the use of more complex interpolation formulae such as Sterling’s method to calculate fitst derivatives is unneomsary. Because of their association with the % error of each variable these terms represent an error mugn@cution or EMAG factor related to that particular variable. Thus a tabulation of these
EMAG factorscan kc used to assessthe potential contribution of each variable to the final error in the parameter G. Alternatively if the errors are known or estimated, variance contributions can be calculated as the square of the product of the % error and the EMAG factor. Each pair of variables have a covariancc, expressed as a correlation coefficient, hut in practice many of the variables are not correlated and py = 0. Thus for the double summation in Eon. (7) only variables with significant correlation coefficients (known or estimated) need be evaluated. Having obtained the necessary variance and covariance terms the summations are carried out to calculate a total variance for the parameter. The vahres of other parameters which are also fun~o~ly related to the variables in the G function (i.e., B = h(x, , x2, * - * x,,) may be required. As G and B are functions of the same variables they are correlated and have a covariance given by:
Estimated pii = ruv=
Equation (2) may be written in a form that is independent of the units of measure by substituting Eon. (3) and expressing
It is more convenient to express thii covatiance as a correlation coefficient for parameters G and B. If we also express the errors in % using Eqn. (5) and express the partial differentials as ratios of increments as was done in Eqn. (7) then the correlation coefficient is:
Numerical error analysis applications
2131
surements. Also, if a normal distribution is an adequate approximation most of the higher moments are zero. A detailed study of the significance of the omitted higher order terms concluded that the approximation “turns out to be better than expected” (TUICEY, 1967, This equation is of a similar form to Eon. (7) and the EMAG factors,variance and covariance components may be evaluated as outlined above. To evaluate any number of parametersand their correlation coefficients it is necessaq to repeat the complete calculation of the parameters n + I times where n is the number of variables. In each cycle the EMAG factor for each variable is determined. These calculations can be expressed and solved using matrix notation (CLIFFORD,1973) but the presentation given here clearly outlines the approach to the numerical calculations. A general algorithm of this procedure has been written in Hewlett Packard (200/300 series) Basic and allows the evaluation of up to 8 parameters having up to 25 variables, though this may be increased if required. The algorithm can be added to existing programs with a minimum of modifications.
The propagation of errors, fluctuations and tolerance. Unpubl. tech. reports 10, I 1, 12. Princeton University
Press). Ku ( 1966) suggests in practice, that the higher terms in the error propagation equation can be neglected if the second and higher order partial derivatives, when evaluated at the means, are small and that the standard deviations of the variables are also small. For small standard deviations Ku suggests values of 0.1% or less, though examples given below suggest standard deviations of several percent are acceptable in most cases. For the requirement of small higher partial derivatives Ku does not suggest any values but the use of numerical calculations permits further tests to assess the significance of the higher order terms, and also large variable ACCURACY OF THE ERROR PROPAGATION errors, to the accuracy of the calculated parameter erEQUATION rors. The mathematical derivation of the error propagaThese tests involve comparing two methods of caltion equation involves a number of assumptions and culating EMAG factors. The previously described calapproximations, which affect the validity of the final culation of the EMAG factors involves successive .O1% error estimate in calculated parameters, and of which deviations in each variable. To obtain a variance conthe user should be aware. In addition to the numerical tribution the squared product of a variable’s EMAG approximations to determine the &t partial derivatives factor and estimated % error is calculated (Eqn. 7). discussed in the last section the other main assumptions The estimated % errors of the variables may be as large are: as several percent in contrast to the small deviation used for the EMAG calculation. An alternate value for (1) The mean values of the variables and calculated the EMAG factor can be obtained by replacing the parameters am derived from approximately normally .Ol% deviation with the estimated SJ error of the distributed populations. variable. These new alternate EMAG factors can then (2) The error propagation equation is an approxibe compared with the original factors to assess if the mation from a Taylor series expansion about the point function is approximately linear (i.e., constant EMAG at which each of the component variables take on their factors or partial derivatives) at the point of evaluation mean value. Terms in the series higher than the first and over the range of the error limits. If an EMAG partial derivatives are considered negligible and thus factor is not constant then the linear approximation is omitted from the analysis. inadequate and higher order terms may be significant. This procedure can be carried out for all variables with Empirical observations have shown that the normal distribution is a good representation for many popua check on the linearity of the calculations to some lations in physical science. Even if it is only an ap arbitrary level. The computer subroutine flags differences greaterthan 5%in the two methods of calculation proximation the central limit theorem states that the sampling distribution of the mean of n independent of EMAG factors and thus highlights potential probobservations of any distribution with finite mean and lems. Analytical formulations of error propagation are variance approaches a normal distribution as n ap incapable of assessing these non-linearity effects. Calculations with alternate EMAG factors can prc+ proaches infinity. For a parameter dependent on four or more variables from different populations the pa- vide some estimate of the uncertainty of calculated parameter errors. The algorithm first calculates alterrameter mean tends to be normally distributed even for small n (EISENHART, 1963; HAHN and SHAPIRO, nate EMAG factors by increasing the mean of each variable by its estimated error limit. Variance contri1967). butions are then calculated and summed to give one The omission of the higher order terms in the Taylor error limit on the parameter. The procedure is then series is usually a valid approximation for most functions. In the case of linear functions (e.g. G = Ui + bX_) repeated but with alternate EMAG factors calculated by decreasing the mean of each variable by its estimated the higher terms vanish and Eon. (2) is exact. inclusion of the higher order terms is usually not very practical error limit. This provides another limit ofthe parameter error. If the function is linear in all variables the two as they contain the higher moments of a distribution (the variance is the second moment; HAHN and SHA- limits will be equal to the original calculation and the PIRO, 1967) which are not available for many meaerror calculation is exact. Non-linear behavior results
J. C. Roddick
2132
in a range in the error limits, however, these limits are only approximations as the method of calculation with alternate EMAG factors presumes a positive correlation in the determination of the Eh4AG factors. This arises from the fact that the calculation sums the variance contributions derived from alternate EMAG factors calculated by changing (increasing or decreasing) each variable in the same way. If several variables cause the function to behave in a non-linear fashion the calculated parameter errors will be poor approximations as evaluation also holds all other variables at their mean values. The net effect of simultaneous random variation of several variables could be to enhance or cancel the non-linear effects. If the linear approximation is inadequate the correct variance would be obtained only if a parameter were a Bmction of a single variable. For situations where the linear appro~mation is inadequate Monte Carlo simulation may be a better approach to the error evaluation though it would be easiest with functions with uncorrelated variables. EXAMPLES OF ERROR ANALYSIS In this section three examples of error propagation
analysis from thermodynamics and geochronology are given. The examples range in increasing order of complexity from the determination of parameters of linear equations to a relatively complex example where the equations are not linear, some variables are correlated and a correlation coefficient for two parameters is required. In applying the computer routines it is necessary to define a subroutine containing the calculations of interest as they would normally be computed with measured variables with known errors and any correlations among variables. The particular variables which have error estimates are defined by equating them to a working array, while results to be calculated are defined in another array. Correlated variables are specified as array element pairs. For any given set of variables, errors, and correlations a call to the error calculation routine will produce the calculated results with their associated errors and correlation coefficients between results for that particular data set. A table of indi~d~ EMAG factors or variance contributions for each variable may be listed, and any significant non-liner EMAG calculations noted, along with the uncertainties in the errors due to non-linearity. The examples which follow display typical tables of error analysis.
2 chlorite + 8 kyanite -t 1 I quartz = 5 con&rite + 5.5 Hz0
the equation for determining the Gibbs free energy of cordierite is (BIRD and ANDERSON,1973): 5 AGt,, = AG,, + A&J T- 25) -Al’,@‘-l)-5SCPHrO
(11)
where: AGf = standard
P= T= G*H20 = AT, =
A&$ =
molar Gibbs free energy of formation from the elements at 298”K, 1 bar (standard state); r = reactants, equilibrium pressure, equilibrium temperature, a form of Gibbs free energy of I&O at P and T (USHER and ZEN, 197 I), molar volume change of the solids in a reaction, third law molar entropy change of the solids and for Eqn. (10) is:
asf,~=5~f.~-2hs~chl-s~f.Ly-
ll~fmz.
(12)
An analogous equation may be written for AV,. Data for this reaction, with 2 standard errors of the means are listed in BIRD and ANDERSON( 1973). Table I details a numerical error analysis of these data at 7000 bars and 65O”C, the same P, T conditions as used by ANDERSON(1976) in his Monte Carlo analysis The table presents all 14 variables with the error prop agation of each variable expressed both in terms of the EMAG factor and the variance. EMAG factors may be positive or negative and reflect the resultant fractional change in AG- for positive change of a variable. Many of the EMAG factors are small (~0.1) and when coupled with the errors on the variables contribute essentially nothing to the variance of AG,_. This ability to separate error components highlights the influence of each variable on the total variance of a parameter. Equation (11) is very nearly linear in this case (7000 bars, 65O’C) as the P and T variables contribute little to the total variance in Table 1. This means that the error propagation calculation can be regarded as es-
VARIABLE NAME
PARAMETERI VALUE: 96 ERRORc
-ai ,219
% ERROR’
EMAct
96’ VAWANCE
-.02?? a022 -.0164 .I760 .w73 -.0626 .ZlSS
.OOLI .OoOO .oocQ .OZU .0018 .0143 .0019
GZ,O
1.34 I.&3 :f’
6=ci11 A=b 2”s
.OWZ 1.91 .2
TmP
Free energy of formation of cordierite Tlte free energies of formation of mineral phases are required for the understanding of chemical reactions in rock-forming systems. BIRD and ANDERSON( 1973) used hydrothermal equilibrium data to calculate the free energy of cordierite and later ANDERSON(1976) applied a Monte Carlo technique to assess the errors in the calculation. For the reaction:
(10)
A&W
.221
.I223
.0007
igz
.7la :YF .0336 .I61 .a036
-.0287 - .0577 -.OOa .OL87 -.Oo% -.OO@O
.0021 .OWO .OOOfl .OOOO .OOOO .cmO
6 Vch, AVcw 6Vky AVsxt~
’ = Errtmattwostmdudemm t
= Error rrugificaianof
nra
in cakulati~
ACf,m
Numerical
sentially exact. The calculated value (-2089.15 Kcall mole) is the same as that obtained by BIRD and ANDERSON (1973) while the error of .219% (2 sigma) is slightly larger than the value of .206’S obtained by hand calculations (BIRD and ANDERSON, 1973) and the value of .2 12% by the Monte Carlo method (ANDERSON, 1976). The numerical calculation has thus obtained an exact result in 15 calculation cycles where the Monte Carlo employed 862 cycles. The chlorite variables AGfchl and A&,,,, dominate the variance in the final error in AGr.,,. Table 1 thus shows that improved measurements in these variables would contribute most to an improvement in the free energy calculation of cordierite. Correlation in thermodynamic calculations
ANDERSON(1977) discussed the general characteristics of error propagation in thermochemical calculations of mineral equilibria such as Eqn. ( 11). He provided a clear example, with the calculation of enthalpies of sanidine and muscovite, of how data are often interdependent with highly correlated relations that unfortunately are usually ignored when applying the data to mineral equilibria calculations. In solution calorimetry the calculation of the enthalpy of formation of a single compound (Mr) involves the measurement of a number of heats of reaction. For muscovite 12 independent measurements are required while for sanidine 11 measurements of heats of reaction plus its heat of solution are required. Eleven of these measurements are common to the A?& of both minerals. A thermochemical calculation may require the enthalpy of the reaction muscovite = sanidine + corundum + water.
2133
error analysis applications
(13)
In this case the algebraic sum of enthalpies from the four components are required. In calculating the uncertainty in this enthalpy of reaction the covariance in the errors of muscovite and sanidine, as a result of the common source of their enthalpies, should be considered Anderson’s procedure was to calculate the variance of (Ai&*, - A&,,,) from the original 13 independent reactions. He emphasized that, while providing the correct variance for these two minerals in this reaction, this procedure cannot be used generally. In calculating the difference of AHr six of the 13 reactions cancel out and do not contribute to the variance of the difference. In contrast if the sum of enthalpies of muscovite and sanidine were required the individual heats of reaction do not cancel and thus result in a greater variance. Numerical error propagation, however, can provide a general solution to this problem. Using the 13 independent reactions and the linear summation equations for AHCW and AHcmur given in Appendix 2 of ANDERSON(1977) the errors and the correlation coefficient for these enthalpies may be calculated.ThevaluesareA&,,,=-1421.23+-1.192 Kcal and AHr_ = -946.453 + -1.098 (2 sigma errors), the Same as determined by Anderson. In addition
the correlation coefficient between these two enthalpy errors is .9 19 and confirms that the errors are highly correlated. The determination of the correlation coefficient provides a value which can be used in other calculations of mineral reactions regardless of whether and AHr,, are added or subtracted. For the Mr.,,, example, for Eqn. ( 13) the errors and correlation coefficient can be included in the numerical error analysis to calculate the enthalpy of reaction. This can be done by specifying which variables are correlated and the value of the correlation coefficient. The resulting covariance will then be propagated into the calculated error on AZ& This example shows how numerical error propagation may be used in several stages. The resulting parameters with their errors and correlation from one calculation may become the variables in a subsequent error analysis. U/Pb ratios for concordia diagrams
In U-Pb geochronology in order to calculate the age of mineral phases that contain U but very little Pb when they are formed (termed common Pb) the amounts of mdiogenic 207Pb(20’Pb*) and x”%%(%Pb*) in a sample are required. Often the ratio of m7Pb*/ 206Pb* provides an estimate of the age, but in most cases there is a loss of radiogenic Pb subsequent to crystallization of a mineral. In these cases it is necessary to calculate ratios of zo7Pb*/238U and 2osPb*/23sU of several sub-samples of a mineral and assess their significance on a concordia diagram (see FAURE (1977) for principles). Thus these U/PC, ratios and their errors for several sub-samples are required for subsequent assessment, often with two-error regression treatments (e.g. YORK, 1969). In general the total Pb in a mineral analysis consists of three components: radiogenic Pb, common Pb and blank Pb added during sample processing. If samples are analyzed by the isotope dilution technique with a mixed isotopic spike containing (for simplicity in this example) pure *05Pb and pure 235U, then the concentrations of the Pb and U isotopes are (abbreviating the isotopes to their masses): 206Pb* = (206/205),(205,)
- 206u - (206/204),(204),,
207Pb* = (207/205),(205,)
( 14)
- (207’204)bi(206u) (206/204), - (207/204),(204,)
(15)
238U = (238/235),/ [ 1 - (238/235),/Ku](235,)
- 238,,,
(16)
where the subscripts sp, m, bl, cm refer to spike, measured, blank and common Pb or U isotopic compositions respectively. Ku = 137.88 = (238/235) in the sample. 205,, 238, = the amounts of spike *05Pband 238Urespectively in the mixture.
J. C. Roddick
2134
PARAMETER: VALUE: % ERROR: VARIABLE NAME
% ERROR
2071206,
.2
3 4 5 6 7
;;;g& 238l235,,, pb Frac. 2061205sp ;3;g$P
.3 .5 .I
9 LO II
Fo Blank 206/2wbl 207/20$,~ 206/2Wc,,, 207/2%,
8
f:
30 .03 30 .05 30 3 .7 2 .4
2imt231 .2019 .a25
207*1235 5.598 1.21
%’ VARIANCE .0640 .I439 .0148 .0103 .0030 .woo .0063 .0367 .2251 .0181 2 .I579 Z
207*/206* .2011 .783
COR
COEFf .7660
CONTRIBUTIONS 2 .4940 .4404 .0103 .0484 2 .0063 .0367 .2136 .0021 .0330 .A76
.O#O .I046 -2936 .0:72 .OOw Z . of2 .0326 .0330 .I579 .I876
Z .2666 .0808 .0103 .0121 .:063 .0367 .2192 -.0062-
-.oow :
COR. COEF. v ( NO’S) 10, 11 12, 13
% ERRORS 3, .7 2, .4
.; .7
TOTAL VARIANCE:
t Z +, +
= = = =
.6802
206*/238 ERROR -.I0 ERROR .I1
IN PARAMETER
= (204/205),205,
.0171 .I205
.6137
.7636
ERROR3 FOR a/- 1NPUT ERROR3
207’1235 -.I4 .I4
207+1206* -.37 .39
206bl = the amount of mPb in a separate blank analysis processed at the same time as the sample. 204, = the amount of 204Pbinitially in the sample, where: 204,
I.4608
-.wM) -.2409
Correlation Coefficient for first tvm parameters EMAC factor and variance wle zero fa this variable EMAG factor subject to > 5% variation in calculation far +/- Input Errors EhtAG factor increase (- = EMAC decrease) for input Erra incrure
96 UNCERTAINTY
+ -
%’ COVARIANCE -.Mw)6 -.0118 . 0000 0.0000
- 206u/(206/204)u.
( 17)
Equations ( 14 to 17) and the value of Ku may be combined to obtain the required U/Pb ratios. Table 2 presents a numerical error analysis of the U-Pb ratios and m7Pb*/206Pb* using all the variables required in Eqns. (14, 15 and 16). In addition 3 other variables are listed: the (238/235), of U is normalized for mass fractionation’ using a double 233-235 spike and so (233/235), is listed; measured Pb ratios, subject to isotopic fractionation, are normalized thus Pb fractionation (Pb. Frac.) is a variable; and the 2osPb spike contains about 5% 206Pb so the spike 206Pb/205Pbis listed though itscontribution is insignificant in this example. Only variance contributions are listed in the table but a distinction is made between very small variances (.OOOO)and those that are zero because EMAG factors are zero (Z). LUDWIG (1980) presented an analytical solution for error propagation in U-Pb ratios similar to that given here and a worked example. The data in Table 2 are the same as Ludwig’s with about 10 ng of sample Pb and 1 ng of blank Pb. The parameter error estimates and correlation coefficient are in good agreement with Ludwig’s with slight differences related to different
’ The (238/235), is arbitrarily nomalized to a fixed value of 233/235 during U analysis. Consequently error in the analysis is expressed only by error in (238/235),.
starting points; Ludwig assumes that the concentrations of 206, and U were previously determined with known errors. Here the error analysis used mainly isotopic ratios derived directly from mass spectrometer analyses. A complete analytical solution to Eqn. (6) using Eqns. ( 14 to 17) and mass spectrometer data would be tedious to derive and subject to algebraic errors while the numerical procedures can overcome these difficulties. In this example correlations between two pairs of variables are included. The common Pb isotopic compositions and the blank isotopic compositions have been given the same correlation coefficient of 0.7 as in Ludwig’s example. These correlations have little intluence on the variance of the U-W ratios but &nihcantly reduce the variance on the 207*/206* as is reflected by the negative covariance components. This is important in the case of a sample which does not have significant Pb loss since its age and error may be calculated directly from the 207*/206* ratio and its error. The uncertainties in spike U and Pb concentrations have not been included in the errOranalysis for two reasons. Because a mixed spike is employed only its U-Pb ratio is required to determine sample U-Pb ratios. While the quantity of spike used is required to make blank corrections, the contribution of uncertainty in this factor (weight X concentration) to the errors in UPb ratios is not significant. Secondly, different samples regressed on a coneordia plot are analyzed with the same mixed spike, so should the calibration of the spike U-l% ratio be biased (within calibration errors) all samples will have the same bias. Thus the uncertainty in the spike U-Pb ratio should be included in the regression procedure and then, only for comparison of ages with data from other laboratories using a different U-Pb spike.
2135
Numerical error analysis applications
These calculations of U/W are derived from equations that are not linear and so significant approximations might be expected especially in the case of error limits as high as +30% as listed for the blanks and Pb fractionation. Using the procedure outlined in the previous section on accuracy all EMAG factors are linear to better than 3% except for the (206/204A, variance contribution (>5%) to the correlation coefficient (Table 2). However, this contribution to the total variance is small while the uncertainty in the propagated error estimates of the three parameters ranges from
use of these parameters in geological and geochemical models. Acknowledgements-M. J. Dalwitz is thanked for unknowingly enlightening me, more than 10 years ago, on the numerical approach to error analysis. Critical comments from Terry Gordon and Rob Berman greatly improved the text and ensured that erroneous ideas did not creep into the mathematical development. An anonymous reviewer is thanked for correcting some misconceptions in the U-Pb example. Editorial handling:D. M. Shaw
only 0.10% to .39% (Table 2). Thus non-linearity is very minor and the error calculations are considered
to be valid. These examples demonstrate the general application of the technique. More complex problems involving additional parameters and correlations are treated in a similar way. While non-linearity in these examples did not produce problems in the error analyses other applications may have significant non-linearity. An example of this sort of problem in 40Ar-39Ardating will be given elsewhere (RODDICK, in prep.). CONCLUSIONS The application of rigorous error analysis to many analytical problems in geology is seldom a usual procedure. There are several reasons for this but the main one is a lack of a technique that can be applied to a range of problems. Numerical error propagation can provide the required general technique. It can be ap plied to complex problems where analytical procedures are difficult and is faster and more informative than a Monte Carlo method. In addition, the procedure can also test whether some of the assumptions involved in the error propagation equation have been fulfilled and thus if the error estimates are- valid. The equations required for numerical error propagation, including correlations are given and the method of calculation discussed. Examples from thermodynamics and geochronology outline the procedures and the advantages of the numerical technique. They also highlight the importance of correlation among errors of variables in the problem and among errors in the final results that are calculated. An algorithm written as subroutines in HP Basic was employed in the calculations and listings of the program for each example are available upon request. Wider application of rigorous error analysis to geological problems will provide a better appreciation of the uncertainties on parameters derived from complex problems. This will permit a more realistic assessment of the conclusions to be drawn from the
RJCFJSRJZNCE!S ANDERSONG. M. (1976) Error propagation by the Monte
Carlo method in geochemical calculations. Geochim. Cosmochim. Acta 48,2309-23 I 1. ANDERSONG. M. (1977) The accuracy and precision of Calculated mineral dehydration equilibria In Thermodynamics in Geology,pp. I 15-135. D. Reidel, Dordrecht-Holland. BIRDG. W. and ANDERSONG. M. (1973) The free energy of formation of magnesian cordietite and phlogopite. Amer. J. Sci. 273, 84-9 I. BIRGER. T. (1939) The propagation of errors. Amer. Phys.
Teacher 7,351-357. CLFFXXDA. A. (1973) Multivariate ErrorAnalysis. J. Wiley and Sons, 112~~. DALRYMPLEG. B. and LANPHERE M. A. (1971) *Ar/=Ar
technique of K-Ar dating: a comparison with the convcntional technique. Earth Planet Sci. L&t. 12,300-308. DAVISJ.
C. (1986) Statistics and Data Analysis in Geology.
J. Wiley and Sons, 55Opp. FAUREG. (1977) Principles of Isotope Geology.J. Wiley and Son%464pp. E~SENHART C. (1963) Realistic evaluation of the precision and accumcy of instrument calibration systems. J. Res. NBS
C. Eng. Instrum. 67C, 161-188. FISHERJ. R. and ZEN E-AN (197 I) Thermochemical calculations from hydrothermal phase equilibrium data and the bee energy of formation of HrO. Amer. J. Sci. 270,297-
314. HAHN G. J. and SHAPIROS. S. (1967)
Statistical Models in Engineering.J. Wiley and Sons, 355pp.
HOLLANDT. J. B. and POWELL R. (1985) An internally consistent thermodynamic dataset with uncertainties and carrelations: 2. Data and results. J. Meta. Geol. 3.343-370. KU H. H. (1966) Notes on the use of propagation of error formulas J. Res. NBS C. Eng. Instrum. 7oC, 263-273. LUDWIGK. R. (1980) Calculation of uncertainties of U-Ph isotope data. Earth Planet. Sci. Lett. 46.2 12-220. POWELL R. and HOLUND T. J. B. (1985) An internally consistent thermodynamic dataset with uncertainties and correlations: I. Methods and a worked example. J. Meta. Geol.
3,327-342. PUGH E. M. and WINSLOWG. H. (1969) The analysis of physical measurements. Addison-Wesley, 246~~. REEKC. E. (1984) Error propagation calculations. Geochim.
Cosmochim. Acta 48,2309-23 I I. YORK D. (1969) Least squares fitting of a straight line with
correlated errors. Earth Planet. Sci. Lett. 5, 320-324.