Another look at the Grubbs estimators

Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

Contents lists available at SciVerse ScienceDirect

Chemometrics and Intelligent Laboratory Systems journal homepage: www.elsevier.com/locate/chemolab

Another look at the Grubbs estimators F. Lombard a,⁎, C.J. Potgieter b, c a b c

Centre for Business Mathematics and Informatics, Northwest University, Potchefstroom, South Africa Department of Statistics, University of Johannesburg, South Africa Institute for Applied Mathematics and Computational Science, Texas A&M University, United States

a r t i c l e

i n f o

Article history: Received 1 July 2011 Received in revised form 22 September 2011 Accepted 28 September 2011 Available online 8 October 2011 Keywords: Grubbs' estimators Instrument precision Online analyzers

a b s t r a c t We consider estimation of the precision of a measuring instrument without the benefit of replicate observations on heterogeneous sampling units. Grubbs (1948) proposed an estimator which involves the use of a second measuring instrument, resulting in a pair of observations on each sampling unit. Since the precisions of the two measuring instruments are generally different, these observations cannot be treated as replicates. Very large sample sizes are often required if the standard error of the estimate is to be within reasonable bounds and if negative precision estimates are to be avoided. We show that the two instrument Grubbs estimator can be improved considerably if fairly reliable preliminary information regarding the ratio of sampling unit variance to instrument variance is available. Our results are presented in the context of the evaluation of on-line analyzers. A data set from an analyzer evaluation is used to illustrate the methodology. © 2011 Elsevier B.V. All rights reserved.

1. Introduction The problem considered in this paper, first discussed by [5,6], involves estimation of the precision of a measuring instrument without the benefit of replicate observations on heterogeneous sampling units. In our applications, the instrument is an on-line analyzer that measures coal quality characteristics in real time, as opposed to a rather time consuming method consisting of mechanical and/or manual sampling from a coal stream followed by chemical analysis in a laboratory. Replicate measurements cannot be made because under normal operating conditions coal on a conveyor belt can pass through the analyzer only once. Furthermore, there is substantial fluctuation in the average daily values of the measured quality characteristic. Nonetheless, some form of repeated measurement is possible in that the characteristic of interest, e.g. ash content, can be measured simultaneously by one or more additional instruments. For instance, it is possible to obtain every day two samples of coal that have passed through the analyzer (the “test instrument” which gives a reading Y), one sample coming from a mechanical sampler (“reference instrument 1”) and another by manual sampling that involves stopping the conveyor belt at regular time intervals (“reference instrument 2”). The two samples are then submitted for chemical analysis, yielding respective readings X1 and X2. The main focus of interest is in assessing the measurement precision of the analyzer.

If n independent sets of readings involving n different units (daily coal samples) are obtained, the data consist of a n × 3 matrix with jth row (Yj, Xj1, Xj2), the n rows being statistically independent. An underlying structure is assumed, namely: Yj ¼ Tj þ ej ; X1j ¼ μ1 þ Tj þ e1j ; X2j ¼ μ2 þ Tj þ e2j

ð1Þ

where Tj has mean ν and variance τ2, μ1 and μ2 are fixed (nonrandom, but unknown) numbers and where the measurement errors ej, e1j and e2j are independently distributed with zero means and variances σ 2, σ12 and σ22 respectively. The random variable Tj represents the “true” value of the characteristic which varies from sample to sample. The variance of Tj, denoted by τ2, is commonly referred to as the “product variance”. The parameters μ1 and μ2 represent unknown additive measurement biases while the random variables ej, e1j and e2j represent instrument measurement errors. [5] proposed the sample covariance SY − X1, Y − X2 between the data vectors Y − X1 and Y − X2 as an unbiased estimator of σ 2. Since this sample covariance can take negative values, often with quite large probability, the estimator n ˜ 2 ¼ max SY−X ; σ 1

o

Y−X2 ; 0

ð2Þ

is used instead. This estimator is asymptotically unbiased and normally distributed with mean σ 2 and variance ⁎ Corresponding author. E-mail address: [email protected] (F. Lombard). 0169-7439/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.chemolab.2011.09.013

  n   o ˜ 2 ¼ n−1 2σ 4 þ σ 2 σ12 þ σ22 þ σ12 σ22 : var σ

ð3Þ

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

However, in small samples the estimator can be quite severely biased and then the asymptotic standard deviation provides a poor indicator of its accuracy. The estimator (Eq. (2)) has been applied in a variety of fields such as climatology [3], agriculture [9] and on-line coal quality analysis [10]. In fact, the “Grubbs method” seems to be firmly established in the coal industry as the preferred method for estimating analyzer precision and has been incorporated into at least one international standard — see [7]. Substantial generalizations of the methodology have been made since Grubbs' original work. In this regard, the work of [1,8] is particularly noteworthy. These works focus for the most part on questions of hypothesis testing, such as whether or not the precisions of multiple instruments may be regarded as being identical or whether one or more of a number of instruments are biased relative to the others. However, these methods are not of direct relevance in the evaluation of online analyzers. First, we typically know a priori that the analyzer delivers substantially greater precision than a process consisting of mechanical or manual sampling followed by chemical analysis. The question is rather whether the increased precision is of an order that would justify the cost of acquiring and maintaining such an instrument. Second, sampling system bias is typically not problematic because evaluation of an analyzer invariably consists of two phases: calibration followed by an assessment of precision. If an analyzer cannot be calibrated to produce results that are unbiased relative to those produced by the reference instruments, the evaluation process terminates in the first phase. Needless to say, if the instrument does not hold its calibration during the precision assessment phase, and this is not noticed, then a whole new data analysis problem presents itself. We forego discussion of such complications in the present paper. Our data set consists of n = 94 triples (Y, X1, X2) of ash contents (% of mass) of low grade coal collected in an analyzer evaluation. (The data are available from the authors upon request.) For these 2 2 2 data, σ˜ ¼ 0:013; σ˜ 1 ¼ 0:101 and σ˜ 2 ¼ 0:105 and Eq. (3) give a nom2 inal estimated standard error of 0.012 for σ˜ , hence an approximate nominal of 95% confidence interval (max{− 0.011, 0,}, 0.037) = (0, 0.037) for σ 2. Given the cost involved in using two reference instruments, it is somewhat disappointing that 94 triples of observations cannot deliver a standard error substantially smaller than 0.012, which could lead to a confidence interval with a positive lower limit. Notice, however, that both reference instruments are evidently producing considerably less precise readings than the analyzer. Thus, it is hardly reasonable to expect a precise estimate of σ to arise from a comparison of the analyzer with these reference instruments. In fact, in large samples the probability that SY − X1, Y − X2 is 2 negative, hence σ˜ is zero, is approximately 1 − Φ(λn− 1) where n   o1=2 2 2 2 2 2 2 2 2 σ2 =σ λn ¼ 2 þ σ1 =σ þ σ2 =σ þ σ1 =σ

o

Y−X1 ; 0

ð4Þ

which has large-sample variance   n   o ˆ 2 ¼ n−1 2σ 4 þ σ 2 σ12 þ τ2 þ σ12 τ2 : var σ

Eq. (5), we see that in the latter the product variance τ2 figures as the measurement error variance of a fictitious “reference instrument”. Thus if τ2 is substantially larger than either σ 2 or σ12, which is typically the case with coal and many other minerals, its presence in formula (5) has the consequence that an even larger n will be needed to produce an acceptably small standard error. In fact, the two instrument variance estimator (Eq. (4)) is very frequently equal to zero in such circumstances. Thus, an important question is whether the precision of the twoinstrument estimator can be improved without increasing the sample size. Clearly, since Eq. (4) is the maximum likelihood estimator in the two instrument setup, an improvement will be possible only if some external information is available about the parameters in the model. We will show that a reasonably precise numerical specification of the ratio τ 2/σ 2 can be used to improve substantially, at no extra cost, the two-instrument assessment of analyzer error variance. In fact, in relatively small samples the proposed estimator generally performs much better than the three-instrument Grubbs estimator. The objective of the present paper is to explore the properties and application of the new estimator. In Section 2 of the paper we define the new estimator and discuss the calculation of its standard error. We illustrate the methodology using data sets collected in the evaluation of online coal analyzers. In Section 3 we use Monte Carlo simulation to compare the sampling distribution of the new estimator to those of the two-instrument and three-instrument Grubbs estimators. The results suggest that in realistic circumstances the new estimator performs substantially better than either of the Grubbs estimators. Some technical results are collected in Appendix A. 2. A new two instrument estimator Because the focus is on a two-instrument method, we henceforth denote X1j by Xj, X2j being redundant. Construction of the new estimator requires ideally that the ratio w0 = τ 2/σ 2 be known. This requirement is akin to a requirement in the well known errors in variables problem [4] that the reliability ratio be known. However, rather than assuming a known value of w0 we will construct our estimator using a “working” value w of w0 and then consider the sensitivity of the estimator to discrepancies between w and the true value w0. The Monte Carlo simulation results reported in Section 3 are mainly aimed at assessing the impact of incorrect specification of w0 on the properties of the estimator that we now propose. The data are n pairs of analyzer — reference instrument readings (Yj, Xj), j = 1, …, n. Set dj ¼ ð1 þ wÞXj −wYj

ð6Þ

where w is the “working” value of the ratio w0 = τ 2/σ 2. Notice that

and where Φ denotes the standard normal distribution function. This follows from the approximate normality of SY − X1, Y − X2. The estimate of this probability in the present instance is 0.14, which follows upon substituting the estimates of σ 2, σ12 and σ22 into the expression for λn. It therefore comes as no surprise that the confidence interval based on SY − X1, Y − X2 has a zero lower limit. The two-instrument version (one reference instrument only) of the Grubbs estimator is: n ^ 2 ¼ max SY; σ

75

ð5Þ

Eliminating one of the reference instruments decreases substantially the cost of the evaluation operation. However, comparing Eq. (3) and

h  i  E dj Tj ¼ t ¼ t so that dj is a predictor of Tj. The idea is now to use the dj values to attempt a stratification of the data into subsets that are more or less homogeneous in respect of the unobserved Tj values. Then product variation within any subset is typically substantially smaller than overall product variation, hence the effect on the variance of the twoinstrument Grubbs estimator calculated on the subset would also be smaller. When the estimators from each of the subsets are averaged an estimator with greatly reduced variance results. To implement this idea we use subsets consisting of two observation pairs each. Arrange the dj values in increasing order of magnitude, d(1) b ⋯ b d(n), and denote the observation pair from which d( j) was calculated by (Y[j], X[j]). For integers j ≠ j′, set  n   o Y½ j −Y½j′  − X½ j −X½ j′  =2 ζj; j′ ¼ Y½ j −Y½ j′ 

76

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80 0.15

The standard error of ζ can be estimated by

0.1

seζ ¼

0.05

 1=2  2 m−2 ∑j¼3 ζj −ζ =ðm−4Þðm−5Þ :

ð10Þ

Formula (10) arises from the fact, established in Appendix A, that ζ1, …, ζm are approximately uncorrelated with a common mean

ζj,j−1

0

E½ζ  ¼ A=ð1 þ wÞ −0.05

and variance 2

−0.1

2

σζ ¼ 2fA=ð1 þ wÞg ; where

−0.15

2 2

−0.2 −5

A¼σ τ −4

−3

−2

−1

0

1

2

3

4

log(d(j)−d(j−1)) Fig. 1. Plot of log(d(2j) − d(2j − 1)) against the Grubbs estimator ζ2j, the two Y observations at d(2j) and d(2j − 1).

2j − 1

formed from

which is the two instrument Grubbs estimator calculated on the two pairs (X[ j], Y[ j]),(X[ j′], Y[ j′]). It is shown in Appendix A that var[ζj, j′] is an increasing function of E[(d( j) −d( j′))2]. To minimize the variance we therefore choose j and j′ such that |j−j′|=1. Thus, we form m/2 pairs {d(2j − 1), d(2j)}, j=1,…,m and eliminate the first two and last two of these from further consideration. The justification for this elimination is that these are often the pairs with the largest values of (d(2j) −d(2j − 1))2, hence the largest (positive or negative) values of ζ2j, 2j − 1. Setting ζj :¼ ζ2j;

2j−1

 n   oo Y½2j −Y½2j−1 − X½2j −X½2j−1 =2 ð7Þ ¼ Y½2j −Y½2j−1

for j=3,…,m−2, the proposed estimator is   2 ^new σ ¼ max 0; ζ ;

  2 1−ρY;d ;

ð8Þ

ρY, d denoting the correlation between Yj and dj. It is perhaps worthwhile illustrating these points with a numerical example. Fig. 1 shows for w = 112 a plot of the 47 pairs (log(d(2j) − d(2j − 1)), ζ2j, 2j − 1), j = 1 …, 47 from the analyzer data set referred to in the Introduction when the mechanical sampler is used as the only reference instrument. There is a large negative estimate at j = 47 which would not be included in the calculation of Eq. (9). The relatively small variability among the remaining 46 pairs is quite striking, especially in light of the fact that the two instrument Grubbs ^ ¼ 0 with a (formal) standard error of estimator gives an estimate σ 0.045. Notice that the majority of these 46 pairs produce ζ values less than 0.045 in absolute value. Hence their standard deviation would be substantially smaller than 0.045 and the standard error of their average (i.e. of the new estimate) will be even smaller. Fig. 2 shows the same plot for w = 89.6 and w = 140. Evidently, the pattern is similar over a range of w values. If w = w0, then ρY, d = 0 and A = σ 2(1 + w), so that E[ζ] = σ 2 and there is no estimation bias. In general, let w = δw0. It is shown in Appendix A that for “large” m and for δ “close to” 1, h i 2 E ζ ≈σ ð1−ðδ−1ÞÞ;

ð11Þ

h i 2 so that the relative bias E ζ =σ −1 of the estimator is proportional to

where ζ ¼ ðζ3 þ ⋯ þ ζm−2 Þ=ðm−4Þ:

ð9Þ

the relative error δ − 1 in the specification w: over (under) specification of w0 will lead to under (over) estimation of σ 2. Thus, the extent

0.15

0.15

0.1

0.1

0.05

0.05

0

0

−0.05

−0.05

−0.1

−0.1

w=89.6

w=140

−0.15

−0.2 −6

−0.15

−4

−2

0

2

4

−0.2 −4

−2

0

Fig. 2. The same plots as in Fig. 1 for w = 89.6 and w = 140.

2

4

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

77

of the bias exhibited by the new estimator depends largely upon the (unknown) numerical value of the ratio w/w0. The usefulness of the method will be greatest if a good “working value”w can be found. In respect of the evaluation of an online gauge, the best source of information is the data generated during the calibration phase, prior to the evaluation phase — see the Introduction above. A very reasonable value of w can be obtained from these data. The bias aspect of the new estimator is discussed in more detail in Section 2. It is also not difficult to show that if all parameters are fixed and the sample size increases without limit, the relative bias and MSE of Grubbs' estimators both converge to zero while those of the new estimator do not, unless δ = 1. Thus, it could be argued that the proposed estimator does not obey a fundamental requirement of statistical estimation theory, namely consistency. Indeed, it does not. However, the raison d'ê tre of the new estimator is the fact that sample sizes at which consistency considerations are relevant are out of the question. A numerical comparison via simulations of the biases and mean squared errors of the three estimators in “small” samples is given in Section 3 of the paper. These simulation results show that in the cases considered, the new estimates almost always lie inside the nominal 95% intervals produced by the two- and three instrument Grubbs estimators and that their mean squared errors are smaller by an order of magnitude than those of the Grubbs estimators, even if the “working value”w is off the true value w0 by as much as 30%. In small samples, for which the new estimator was specifically developed, the confidence intervals produced by the Grubbs estimators are typically inordinately wide.

three instrument Grubbs method would require approximately 188 (= 47 × (.012/.006) 2) triples of observations (2 × 188 = 376 laboratory analyses) to deliver the same standard error as the estimator ζ which uses 94 pairs of observations (94 laboratory analyses). Clearly, a potentially more cost effective approach would be to apply the newly proposed estimator using a mutually agreed upon (between supplier and purchaser) value of w. In practice one would typically be in possession of reasonably accurate estimates of σ12, σ22 and τ 2 based upon historical data. In general, then the sensitivity of the new estimator to various values of w can be explored using calculations such as those that led to Table 1. Alternatively, an overall picture can be obtained by assigning a nondegenerate distribution to the relative error δ = w/w0 to express the uncertainty involved in the estimate of w0. If this distribution is centered at δ = 1 one is in effect saying that w is the best available estimate of w0. From Eq. (11) and the central limit theorem, the distribution of ζ is approximately normal with mean σ 2{1+ (1− δ)w0/(1 + w0)} and vari 2 ance seζ . If a normal(1, ξ2) distribution is assigned to δ, then the

3. Data analysis

where a > 0 and 0 b β b 1/2 are specified numbers. Then ξ = a/z1 − β/2 where z1 − β/2 denotes the 100(1− β/2) quantile of the standard normal distribution. For the analyzer data, suppose it was considered unlikely (probability = 0.10) that the preliminary estimate w = 112 would be more than 20% off the true value w0. Then a = 0.2 and β = 0.10, so that ξ = 0.102, v = 0.0062 and a 95% confidence interval for σ 2 is (0.0018, 0.0262). Returning briefly to Table 1 we notice that for w ≥ 117.6 the new estimates are identical, as are their standard errors. This is easily explained upon writing dj in the form

In our data set the agreed upon working value of w0 was w = 112. Table 1 shows for a range of values of w, up to a 25% deviation from w = 112, the corresponding estimates ζ and Eq. (4) of σ 2 together with their estimated standard errors seζ from Eqs. (10) and (5). We assume in these calculations that the mechanical sampler is used as the sole reference instrument. The estimates and standard errors produced by the three instrument Grubbs estimator (Eq. (2)) applied to the first 47 observation triples, now including the manual sampling operation, are also shown for comparison. (The total cost involved in gathering and chemically analyzing 94 mechanical samples is more or less commensurate with that involved in gathering 47 pairs of mechanically and manually extracted samples.) Notice that the two instrument Grubbs estimator (one reference instrument) produces an estimate of 0 for σ 2 and that the three instrument Grubbs estimator (two reference instruments) has a standard error of the same order of magnitude as the estimate of σ 2. In fact, all the new method estimates shown in the table fall well inside nominal 95% confidence intervals for σ 2 constructed from the three instrument Grubbs estimate. Some calculation also shows that the

Table 1 Estimates of σ 2 and standard errors for the analyzer data. Method

w

Est. × 104

Std. err. × 104

New

84.0 89.6 95.2 100.8 106.4 112.0 117.6 123.2 128.8 134.4 140.0 1 ref. instr. 2 ref. instr.

170 166 166 159 159 144 126 126 126 126 126 0 130

41 41 41 40 40 40 39 39 39 39 39 450 120

Grubbs Grubbs

distribution of ζ is approximately normal with mean σ 2 and variance  2 n o2 v2 ¼ seζ þ ξζ w=ð1 þ wÞ . An approximate 100(1− α)% confidence interval for σ2 is then centered at ζ with limits ± z1 − α/2v. An appropriate value of ξ2 can in practice be found by stating the perceived precision of the estimate w in the form Pr ½jw=w0 −1j > a ¼ Pr½j1−δj > a ¼ β;

  dj ¼ Xj þ w Xj −Yj : Clearly, if w is so large that w  min fjXj −Yj jg > max Xj ; j

j

then the ordering of the dj, hence the numerical value of the new estimate, will remain unchanged if w is further increased. 4. Monte Carlo simulations The behavior of the new estimator compared to the standard twoand three instrument Grubbs estimators will now be illustrated using Monte Carlo simulation. The following parameter configuration was suggested by an analysis of the full analyzer data set: 2

2

2

2

σ ¼ 0:013; σ1 ¼ 0:101; σ2 ¼ 0:105; τ ¼ 1:622: 10,000 sets of n (=60, 94) trivariate normal observations (Y, X1, X2) following this parameter configuration were generated using the model in Eq. (1). The new estimate, ζ and its estimated standard error (Eq. (10)) were calculated for each of the 10,000 sets using only the first two components (Y, X1). The first two and the last two of the m/2 pairs {d(2j − 1), d(2j)}, j =1,…,m were not used in forming the estimate. Table 2 gives in each row the mean of 10,000 simulated ζ values, the “true” (i.e. simulated) standard error (se), the average of the standard errors estimated from Eq. (10) (in brackets) and the “true” root mean

78

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

Table 2 Comparison of the new estimator with two Grubbs estimators (σ 2 = 0.013, σ12 = 0.101, τ2 = 1.622, n = 60, 94. w02 = 124.7). Method

New

Grubbs Grubbs

w

87.3 106.0 112.2 118.5 124.7 130.9 137.2 143.4 162.1 1 ref. instr. 2 ref. instr.

n = 60

n = 94

Method

Mean

Se

Rmse

Mean

Se

Rmse

185 153 144 138 130 123 118 113 100 293 166

71 (68) 66 (62) 64 (60) 64 (60) 62 (58) 60 (57) 59 (56) 58 (56) 56 (54) 375 170

90 70 65 64 62 61 60 60 63 409 174

186 151 146 138 130 124 120 114 99 230 149

50 (47) 43 (42) 42 (41) 38 (39) 41 (39) 38 (37) 37 (37) 38 (36) 35 (34) 302 139

75 48 45 39 41 38 38 41 46 318 141

squared error (rmse). Each row in the table corresponds to a different value of w in Eq. (6), namely w ¼ 87:3; 106; 112:2; 118:5; 124:7; 130:9; 137:2; 143:4; 162:1: The first four and last four values correspond to incorrect specification of the ratio w0 = 124.7 by 5%, 10%, 15% and 30% respectively. Also shown in Table 2 are the simulation means, standard errors and rmse of the two- and three instrument Grubbs estimators. For the latter estimator we use only the first n/2 of the n generated triples because with two reference instruments a sample of n/2 (Y, X1, X2) values involves n laboratory analyses which are comparable with the one reference instrument setup involving n (Y, X1) pairs. Tables 3 and 4 give the results when the product variance τ 2 is decreased from 1.622 to 0.101 and then to 0.013 respectively. The values of w02 in the latter two instances are 7.77 and 1.0 respectively. In all three tables the means, standard errors and root mean square errors have been multiplied by 10 4. We consider first the simulation results obtained at n = 60. Notice first that the three configurations differ only in respect of the assumed magnitude of the product variance τ 2. Since the three instrument Grubbs estimator does not involve product variance in any way, we would expect its performance to be similar under all three configurations. This is borne out by the simulation results. The three instrument Grubbs estimator is badly biased and has an excessively large nominal standard error under all three configurations. In the first two configurations, where product variance exceeds analyzer variance by one and two orders of magnitude, the two instrument Grubbs estimator fares almost equally as badly as its three instrument counterpart. Only in the third configuration, where product variance equals analyzer variance, does the two instrument Grubbs estimator exhibit more or less satisfactory performance. We pointed out in the

Table 3 Comparison of the new estimator with two Grubbs estimators (σ 2 = 0.013, σ12 = 0.101, τ2 = 0.101, n = 60, 94.w02 = 7.77). Method

New

Grubbs Grubbs 2

w

5.44 6.60 6.99 7.38 7.77 8.16 8.55 8.94 10.10 1 ref. instr. 2 ref. instr.

n = 60

Table 4 Comparison of the new estimator with two Grubbs estimators (σ 2 = 0.013, σ12 = 0.101, τ2 = 0.0.013, n = 60, 94. w02 = 1.00).

n = 94

Mean

Se

Rmse

Mean

Se

Rmse

176 149 143 135 130 124 120 115 103 143 165

51 (48) 43 (41) 42 (40) 40 (38) 38 (36) 37 (35) 35 (34) 35 (33) 32 (30) 125 170

69 47 44 40 38 37 37 38 42 125 173

177 150 142 136 130 124 120 114 103 139 151

39 (38) 33 (32) 32 (31) 30 (29) 30 (28) 27 (27) 26 (26) 26 (25) 23 (23) 105 141

61 39 34 31 30 28 28 31 36 105 143

New

Grubbs Grubbs 2

w

.70 .85 .90 .95 1.00 1.05 1.10 1.15 1.30 1 ref. instr. 2 ref. instr.

n = 60

n = 94

Mean

Se

Rmse

Mean

Se

Rmse

152 141 136 133 130 127 124 121 113 131 165

42 (40) 39 (37) 38 (36) 38 (35) 37 (35) 35 (34) 35 (33) 34 (32) 31 (30) 71 169

47 41 39 38 37 36 35 35 36 71 173

153 139 137 134 130 128 123 121 113 130 130

33 (32) 30 (29) 31 (29) 28 (28) 29 (27) 28 (27) 26 (26) 25 (25) 24 (24) 58 172

40 31 32 28 29 28 27 27 30 58 172

Introduction that in the two-sample Grubbs method the product variance figures as the measurement error variance of an imaginary third “measuring instrument”. This is manifested in the three tables above. When τ 2 is substantially larger than σ12 (Table 2) the two-instrument Grubbs estimator has a much larger root mean squared error than the three-instrument estimator. But when τ 2 is equal to or less than σ12, (Table 4) the opposite holds. On the other hand, we see that a reasonably accurate specification of the weight w limits bias in the new estimator and that the latter has by far the smallest root mean squared error among the three estimators and in all three cases considered above. Furthermore, the estimates shown in the tables at all values of w lie well within 95% confidence intervals constructed from either the two- or three instrument Grubbs estimators. The overall conclusions are very similar when the sample size is n = 94. Here the two Grubbs estimators show improved performance in terms of bias because of the greater amount of data involved. However, the 50% increase in the amount of data is not sufficient to bring the standard errors of these estimators down to the same levels as those of the new estimator.

5. Summary We have considered the methods proposed by [5] that can be used to estimate the error variance of a measuring instrument when replicate observations cannot be made. A modified version of Grubbs' two instrument estimator has been proposed. Successful application of the method requires a reasonably accurate numerical specification of the ratio of product variance to instrument variance. If this requirement is met to a satisfactory degree the new estimator is superior to both the two and three instrument Grubbs estimators in that it requires a substantially smaller sample size to attain a pre-specified level of estimation precision. The method is particularly well suited to situations where the precision of a new instrument, such as an online coal analyzer, must be estimated and where (i) the available reference instrument has a substantially larger error variance than the analyzer and (ii) the product variance is substantially greater than the error variance of either instrument. In the latter two circumstances the two instrument Grubbs estimator often suggests an analyzer error variance of zero unless the amount of data is excessively large. The two instrument variance estimator proposed in this paper largely overcomes this difficulty. However, we emphasize that the new method proposed here is not a panacea. If the prior estimate of the ratio of product variance to instrument variance is far off the mark, the method could produce a grossly incorrect estimate of gauge error. However, it is seldom the case that reasonably accurate prior information is not available. Hence, it will generally be safe to implement the new method. The potential savings in respect of time and of additional sampling and analysis costs certainly make

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

the method worth considering as an alternative to the threeinstrument Grubbs estimator. Acknowledgment The first author's work was supported by the National Research Foundation of South Africa. The second author's work was supported by Award No. KUS-C1-016-04, made by King Abdullah University of Science and Technology (KAUST). The authors thank two referees for valuable comments that led to the correction of some errors and to an improved exposition of the work.

Here we derive the mean and variance of the random variables ζ2j, 2j − 1 that were defined in Section 1. For notational simplicity we write ζj for ζ2j, 2j − 1. We assume that the measurement errors and product variability are normally distributed. Let σY2, σd2 and σY, d denote respectively the variance of Y, the variance of d and the covariance between Y and d and set ρY, d = σY, d/σY σd. Then we may write Y½ j 

and h i var ζj; j′ h ii h ii ¼ E var½ζ jdð jÞ ; dð j′ Þ þ var E½ζ jdð jÞ ; dð j′ Þ

    2  2  ¼ 2fA=ð1 þ wÞg2 1 þ E dð jÞ −dð j′ Þ C 2 =4A þ var dð jÞ −dð j′ Þ B2 =8A2 :

To minimize the variance we choose j and j′ such that | j − j′| = 1. Then, for 0 b γ b 1/2 and ⌊mγ⌋ ≤ j ≤ ⌈m − mγ⌉, E

Appendix A. Calculations for Eq. (11)

  ¼ ρY;d σY =σd dðjÞ þ Zj

79

  2    2    −2 −4 dð2jÞ −dð2j−1Þ ; var dð2jÞ −dð2j−1Þ ¼O n ¼O n

so that h i n   o −2 B=2A E ζj ¼ A=ð1 þ wÞ 1 þ O n

ð13Þ

and h i n     o 2 −2 2 −4 2 2 C =4A þ O n B =8A : var ζj ¼ 2fA=ð1 þ wÞg 1 þ O n

2 where Z1,…,Zn are i.i.d. normal with zero mean and variance σY2(1−ρY, d) and are independent of d(1),…,d(n) — see [2]. Let j≠j′ be integers between 1 and n. The conditional distribution of Y[ j] −Y[ j′], given d( j) and d( j′)are normal with mean

   ρY;d σY =σd dð jÞ −dðj′ Þ

Furthermore,    2  2  2 cov ζjþ1 ; ζj ¼ B cov dð2jþ2Þ −dð2jþ1Þ ; dð2jÞ −dð2j−1Þ   −4 ¼ O ðδ−1Þn which, because var(ζj) = O(1), implies that

and variance

  −4 : ρζjþ1 ;ζj ¼ O ðδ−1Þn

  2 2 2σY 1−ρY;d :

The conclusion is that even for moderately large sample sizes the ζj are, to good approximation, uncorrelated and identically distributed with mean A/(1 + w) and variance 2fA=ð1 þ wÞg2 . Now set

Set   2 2 A ¼ σY 1−ρY;d ;   B ¼ ρY;d σY =σd ρY;d σY =σd −1 ; C ¼ 2ρY;d σY =σd −1

2

2

w ¼ w0 δ ¼ δτ =σ : Notice that the covariance between Yj and dj is τ 2 − wσ 2, which is zero when w = τ 2/σ 2 : = w0. Thus, when δ = 1 we have ρY, d = 0, hence B = 0, C = − 1 and A/(1 + w) = σ 2, so that E[ζj] = σ 2 and

and  n   o Y½ j −Y½ j′  − X½ j −X½ j′  =2: ζj; j′ ¼ Y½ j −Y½ j′ 

  h i 2    4 2 2 : var ζj ¼ 2σ 1 þ E dð2jÞ −dð2j−1Þ =4 σ þ τ

Now,   Y½ j −X½ j ¼ Y½ j −dð jÞ =ð1 þ wÞ;

In this case ζj estimates σ 2 unbiasedly with a variance that increases with the expected squared distance between d(2j) and d(2j − 1). Finally, we derive Eq. (11). From Eqs. (1) and (6),

so that straightforward calculation gives   i  2 h   E ζj; j′ dð jÞ ; dðj′ Þ ¼ fA=ð1 þ wÞg 1 þ dð jÞ −dðj′ Þ B=2A ;

σY;d ¼ τ −σ w ¼ τ ð1−δÞ ¼ Oð1−δÞ

and

2 while σd2 = O(1). Thus ρY, d = o(1) as δ → 1 follows and, consequently, B = o(1) and C = − 1 + o(1). Then,

2

  i  2 h   2 2 var ζj; j′ dð jÞ ; dðj′ Þ ¼ 2fA=ð1 þ wÞg 1 þ dðjÞ −dðj′ Þ C =4A :

2

σY2 =ð1 þ wÞ ¼ σ 2 ð1 þ w0 Þ=ð1 þ δw0 Þ ¼ σ 2 ð1 þ w0 Þ=ð1 þ w0 −ð1−δÞw0 Þ ¼ σ 2 f1−ð1−δÞw0 =ð1 þ w0 Þg−1   ¼ σ 2 f1 þ ð1−δÞw0 =ð1 þ w0 Þg þ O ð1−δÞ2 :

Hence,    2  h i B=2A E ζj; j′ ¼ fA=ð1 þ wÞg 1 þ E dð jÞ −dðj′ Þ

2

ð12Þ

Substitution into Eq. (13) yields Eq. (11).

80

F. Lombard, C.J. Potgieter / Chemometrics and Intelligent Laboratory Systems 110 (2012) 74–80

References [1] R. Christensen, L.G. Blackwood, Tests for precision and accuracy of multiple measuring devices, Technometrics 35 (1993) 411–420. [2] H.A. David, J. Galambos, The asymptotic theory of concomitants of order statistics, Journal of Applied Probability 11 (1974) 762–770. [3] B. de Smedt, B. Mohymont, G.R. Demarée, Grubbs revisited: a statistical technique to differentiate measurement methods of the degree of exposure of rain gauges in a network, Hydrological Sciences 48 (2003) 891–897. [4] W.A. Fuller, Measurement Error Models, New York, John Wiley, 1987. [5] F.E. Grubbs, On estimating precision of measuring instruments and product variability, Journal of the American Statistical Association 43 (1948) 243–264.

[6] F.E. Grubbs, Errors of measurement, precision, accuracy and the statistical comparison of measuring instruments, Technometrics 15 (1973) 53–66. [7] International Standards Organization, Solid Mineral Fuels — Evaluation of the Measurement Performance of On-line Analysers 15239:2005, Geneva, 2005. [8] J.L. Jaech, Statistical Analysis of Measurement Errors, John Wiley, New York, 1985. [9] H.-P. Piepho, Application of a generalized Grubbs' model in the analysis of genotype–environment interaction, Heredity 73 (1993) 113–116. [10] C.D. Rose, Methods for assessing the accuracy of on-line coal analyzers, Journal of Coal Quality 10 (1993) 19–28.