Analyst, November, 1973, Vol. 98, jy5. 777-781 777 Monte Carlo Simulation of Matrix Correction Effects BY R. J. HOWARTH (Applied Geochemistry Research Group, Department of Geology, Imperial College of Science and Technology, London, SW7 2BP) Monte Carlo simulation is useful for the precise evaluation of the effects of complex systems of matrix correction equations (such as occur in spectro- graphic analysis). If the error distributions for the interfering elements are experimentally determined, that induced by interaction in the correction equation system for the elements subject to interference can be predicted. IN spectrographic systems the observed concentrations of both major and trace elements may have been affected by errors that are dependent on the concentrations of other elements present in the sample.The magnitude ef these errors, and the form of the functional relation- ship between the affected element and the interfering elements, will vary with the nature of the sample, the system of correction equations being generally referred to as “matrix corrections.” This functional relationship can be represented in the general form where Ci is the corrected concentration of the j t h affected element; Cj the observed concen- tration of the j t h affected element; xi the observed concentration of the ith interfering element (i = 1, . . ., T Z ) ; andfi ( ) the generalised functional relationship for the j t h affected element. The form of fj ( ) is generally determined by observations made on known samples spiked with various concentrations of the interfering elements.Any element that has been corrected may, in its turn, be used in further’correction equations, so building up a complex system for the simultaneous correction of all the elements being determined. However, direct appraisal of the effects of such equation systems cannot be carried out in more than qualitative terms. The detailed matrix corrections will normally be specific to particular laboratories as they will depend on the instrumentation, the calibration techniques used and the nature of the samples being analysed. The general principles of emission-spectrographic and X-ray fluorescence analysis have been outlined by Ahrens and Taylor1 and Norrish and Chappell,2 respectively, and matrix correction equations have been used in both emission-spectrographic and, to a lesser extent, X-ray fluorescence analysis. The papers of Tennant and Sewell3 and Leake et aL4 may be cited as examples of typical laboratory applications.The purpose of this paper, however, is to indicate how Monte Carlo simulation can be used as a tool to assess quantitatively the effects of error propagation through a battery of matrix correction equations, an aspect that has not previously been investigated. It will be illustrated with examples drawn from an emission-spectrographic system, although the technique could be applied to any situation where complex inter-element interference occurs. MONTE CARL0 SIMULATION- A Monte Carlo simulation is based upon the development of a mathematical model (the matrix correction equation system) that accurately represents the real-world situation to be investigated. The concepts of Monte Carlo simulation are discussed in a number of texts.6-7 The sources of random error are represented in the model by pseudo-random number generators, random values being drawn from populations that have the same prob- ability distribution and parameters (in this instance the mean and standard deviation) as those in the real-world system.The model is translated into a computer program and simulation runs are conducted by the computer to represent randomly selected real-world trials. A suitable basis for the simulation is the Muller method8 of generating pairs of random uncorrelated values ( A and B) from a Gaussian parent population of mean m and standard deviation s, these values having been determined by observation of the real-world system it is desired to simulate.Let U and V be independent random variables uniformly distributed @ SAC and the author.778 HOWARTH : MONTE CARL0 SIMULATION OF [Analyst, VOl. 98 in the interval (0, 1); then, the probability distribution for U and I/ is described by Oify ( 0 { Oify > 1 P ( Y < y ) = yifOQ-3 . The majority of computer centres have library programs for generating uniform random variables of this type. Now .. * (2) .. * * (3) A = m + s d-2 InU cos (z~v> B = m + s 4 - 2 InU sin ( 2 n ~ ) .. .. Because the Monte Carlo simulation involves random values, the results obtained are subject to statistical fluctuations; thus, the larger the number of trials carried out, the more precise will be the final answer.SIMULATION OF AN EMISSION-SPECTROGRAPHIC SYSTEM- As an example of the application of this method to a real system, we will briefly consider simulation of the matrix correction equation system for the ARL 29000B direct-reading optical spectrograph that is being used by the Applied Geochemistry Research Group for rapid, low-precision analysis of approximately 50 000 stream sediment samples in order to compile a regional geochemical atlas of England and Wales. With such an instrument we have the following defects : spectral interference between the lines present for two (or more) elements; a background effect, principally caused by continuous radiation, scattering, or fine spectra due to molecular emission ; and thirdly, the arc effect, which is an intensification or diminution of the intensity of a given spectral line caused largely by variation in the temperature of the arc as a result of differing rock matrix composition.l It is necessary to evaluate the effect that variations in the major element determinations have upon the trace-element values, by acting through the system of matrix correction equations.Analytical control is based on a series of eight representative natural standards (stream sediment samples from streams draining known homogeneous rocks) and two synthetic standards spiked with either a low or a high trace-element concentration. From the initial period of operation of the spectrograph it was possible to obtain the mean and standard deviation values for the element concentrations in all of the standards, based on a large number of replicate determinations for the major elements (aluminium, calcium, iron, potassium, magnesium and silicon) and each of twenty-four trace elements.If we assume, on the basis of the observed behaviour of the element, that the error distribution for the major elements is Gaussian and that for the ith element it is distributed with mean m, and standard deviation si, then we can simulate the major element variation for any matrix type by substituting the appropriate values of mi and si into equations of the form of equations (2) and (3). The observed uncorrected mean trace-element values [Cj of equation ( l ) ] should also be recorded. For each standard, 1000 simulated sets of major element values were drawn at random by, using the Muller method, and the corresponding matrix-corrected trace-element values [C, of equation (l)] were evaluated.The change, or perturbation, of the initial values caused by the correction equations was then calculated as a percentage ratio. For the j t h trace element the perturbation is given by These data were output by the computer program in the form of histograms, together with the mean, standard deviation and maximum and minimum perturbation values for each trace element for each of the ten standards. As an additional check on the validity of the method, the total percentage of oxide was calculated for each set of simulated major element values. These results had means acceptably close to 100 per cent.for all standards. SIMULATION RESCLTS OBTAIKED- It is to be expected that any correction equation of the formNovember, 19731 MATRIX CORRECTION EFFECTS 779 14 1 \- (I 0 - 50 0 50 100 pu, per cent. Fig. 1. Histogram of percentage frequency ( j ) of perturb- ation values (pa, per cent.) for strontium in the limestone-derived stream sediment control sample (where the ai are coefficients) will yield a Gaussian distribution for C;.B That the frequency distributions for the perturbation values were Gaussian in all instances where the nature of the correction equations would indicate it, was confirmed by testing against fitted Gaussian distributions by using the Kolmogorov - Smirnov statistic.1° This reproductive property of the normal distribution does not apply to non-linear correction equations.For example, Fig. 1 shows the positively skewed perturbation frequency distri- + 1- 0 1 pu20, per cent. I 3 Fig. 2. Variation of pcr- centage perturbation at two standard deviations (pu 20, per cent. ,) with increasing aver- age nickel content in the control samples. Lithologies of the stream sediment source rocks are shown: shale (9); sand- stone (0) ; limestone (0) ; gran- ite (+) ; basic igneous (A) ; and synthetic standard ( x )780 HOWARTH: MONTE CARL0 SIMULATION OF [A~zaZyjst, Vol. 98 bution for strontium in the limestone standard, resulting from a correction equation of the form Sr* = (Sr + a,Ca)/(l + a,Al + a,Ca + a,Fe + a&%) The general reduction in the size of the relative perturbation with increase in the trace- element content is typified by the behaviour of nickel (Fig.2). This element is corrected solely for calcium and therefore shows a Gaussian distribution for the perturbation values. The graph shows the relative perturbation at two standard deviations and indicates that we can expect that only 5 per cent. of the nickel determinations would be perturbed by more than the values shown. While it has to be remembered that each analytical system will yield a unique matrix correction system (and hence perturbation effects), it is of interest to note how the behaviour predicted by simulation compared in this instance with the results obtained in practice. Table I shows that, for the majority of elements, the mean perturbation value for the synthetic standard spiked with low element concentrations is close to zero (per cent.).How- ever, there are considerable variations in the magnitudes of the perturbation standard deviation. The spread of values obtained in a typical set of replicate analyses (made on one day, very much later in the project than the data used for the simulation) is indicated by the precision, defined here as twice the element standard deviation divided by the mean, in terms of concentration. It is clear from Table I that for most elements the effect of the matrix corrections has been to increase the precision value to some extent, and that when this increase has been in excess of a factor of 1.3 compared with the uncorrected value (for silver, arsenic, beryllium, cobalt, lithium, molybdenum, phosphorus, tungsten and zirconium) it correlates very well with high-perturbation standard deviations obtained by the earlier simulation of the system.Some difficulty would therefore be expected in obtaining reliable low-level analyses for these elements. It was found necessary in practice to use alternative methods of analysis for arsenic, molybdenum, cadmium and zinc ; silver, phos- phorus, tungsten and zirconium have not been used in the preparation of geochemical maps. The simulation results have therefore been well borne out in day-to-day experience of the rapid, low-precision analytical system necessary to cope with the very large number of multi-element analyses necessary for a geochemical reconnaissance of this type. TABLE I COMPARISON OF SIMULATED AND ACTUAL PERTURBATION EFFECTS FOR Element (i) 2: * ' ..Ba .. Be .. Bi .. Cd .. G 3 .. Cr .. c u .. Ga .. Li .. Mn .. Mo . . Ni . . P .. sc .. Sn . . Sr .. Ti .. V .. W .. Zn .. Zr .. A SYNTHETIC STANDARD - * - c: PUi .. .. i 3 .. .. 719 0 .. .. 146 0 .. .. 26 -1 .. .. 34 0 .. .. 26 -1 .. .. 16 0 .. .. 190 0 .. .. 76 0 .. .. 18 0 .. .. 47 0 .. .. 60 0 .. .. 6 -2 .. .. 166 0 .. .. 496 2 .. .. 21 0 .. .. 262 0 .. .. 197 1 .. .. 312 0 .. .. 240 1 .. .. 164 -1 .. .. 27 2 .. . . 200 2 uDu$ 92 86 6 17 1 28 62 3 1 0 20 4 49 1 69 12 1 13 2 6 19 46 84 +cl4nc§ 3.32 1-81 1.06 1.44 1.00 1.03 2.23 1.10 1.00 1.00 1431 1.03 1.69 1.03 1-64 1.00 1.07 1.02 1-06 1.2 1 1.69 1.16 2.50 * Mean concentration (p.p.m.) of matrix corrected results (11 samples). t Mean simulated perturbation (per cent.) for matrix correction (1000 trials).; Standard deviation of perturbation values (per cent.) for matrix correction (1000 trials). 9 Ratio of precision [$ = 2 (standard deviationlmean percentage)] for corrected ($c) and uncorrected ($uc) results (11 samples).November, 19731 MATRIX CORRECTION EFFECTS 781 CONCLUSION Visual evaluation of the effects of multi-component correction e-quations cannot be carried out in more than qualitative terms. However, Monte Carlo simulation allows the over-all magnitude of the matrix corrections to be evaluated precisely, and affords a method of comparison between the various corrected elements that helps to assess the relative magnitude of the perturbation effects. The method is easily applied to the most complex of matrix correction equation systems.It is economical in computer time and methods are available for determining the optimum number of trials for the evaluation of a given model.ll It is not intended that the specific results reported here should be applied to other emission-spectrographic systems, but rather encourage the investigation of other analytical system interactions by using simulation techniques when it is of interest to separate the effects of errors in the determination of the uncorrected element values from changes induced through the application of multi-component correction equations. The project of which this paper forms a part has been supported by a Natural Environ- ment Research Council grant for an investigation, under the direction of Professor J. s.Webb, into the applicability of computer methods to the analysis of regional geochemical data. The regional geochemical atlas of England and Wales has been made possible by a major grant from the Wolfson Foundation. The experimental determinations on the ARL 29000B system were obtained by Dr. M. Thompson. Computer time was provided by the Imperial College Computer Service. 1. 2. 3. 4. 5. 6 . 7. 8. 9. 10. 11. REFERENCES Ahrens, L. H., and Taylor, S. R., “Spectrochemical Analysis,” Addison-Wesley Co. Inc., London, 1961. Norrish, K., and Chappell, B. W., “X-ray Fluorescence Spectrography,” in Zussman, J., Editor, “Physical Methods in Determinative Mineralogy,” Academic Press, London, 1967, pp. 161 to 214. Tennant, W. C., and Sewell, J. R., Geochim. Cosmochim. Acta, 1969, 33, 640. Leake, B. E., Hendry, G. L., Kemp, A., Plant, A. G., Harvey, P. K., Wilson, J , R., Coats, J. S., Aucott, J. W., Lunel, T., and Howarth, R. J., Chem. Geol., 1969-70, 5, 7. Wagner, H. M., “Principles of Operations Research,” Prentice Hall, Englewood Cliffs, New Jersey, 1969. Hammersley, J. M., and Handscomb, D. C., “Monte Carlo Methods,” Methuen and Co. Ltd., London, 1964. Naylor, T. H., Balintfy, J. L., Burdick, D. S., and Chu, K., “Computer Simulation Techniques,” John Wiley and Sons, New York, 1966. Dillard, G. M., I.E.E.E. Trans. Inf. Theory, 1967, IT-13, 616. Topping, J., “Errors of Observation and their Treatment,” Institute of Physics Monograph, Siegel, S., “Nonparametric Statistics for the Behavioural Sciences,” McGraw-Hill, New York, 1956. Hahn, G. J., I.E.E.E. Trans. Systems Man Cybernetics, 1972, SMC-2, 678. Received February 19th, 1973 Accepted July 2nd, 1973 Chapman and Hall, London, 1957.