- Email: [email protected]
Estimation of error components in a multi-error linear regression model, with an application to track fitting
NUCLEAR INSTRUMENTS & METHODS IN PHYSICS RESEARCH Section A
Nuclear Instruments and Methods in Physics Research A324 (1993) 317-319 North-Holland
Estimation of error components in a multi-error linear regression model, with an application to track fitting Rudolf Frühwirth '
Institut für Hochenergiephysik der Österreichischen Akademie der Wissenschaften, Vienna, Austria
Received 20 July 1992 We present an estimation procedure of the error components in a linear regression model with multiple independent stochastic error contributions. After solving the general problem, we apply the results to the estimation of the actual trajectory in track fitting with multiple scattering . 1. A multi-error linear regression model Let us consider a linear regression model in which the stochastic deviation from the linear relationship is a superposition of n + 1 independent random variables: y=Ap+e o + - . . + En . Here y denotes the vector of observations, A the matrix of the known linear model, and p a vector of unknown parameters . We assume that dim(p) is less than dim(y), and that the rank of A is equal to dim(p). In section 2 we shall interpret y as the measurement vector belonging to a track and p as the track parameters to be estimated. The measurement errors of the detector will play the role of eo , and the stochastic deviations due to multiple scattering the role of E, . For a given set of observations y, we wish to estimate the unknown parameters p and the actual values of e i = 0, - - - , n . This can be done if there is prior knowledge on the distribution of the e, . We assume that we know the first two moments (expectation and covariance matrix) of all e, : E(e,)=0, var(e,)=V,=Gl', i=0,---,n . The assumption that E(e,) = 0 implies no loss of generality. Let us denote the actual values of e, as e, . Then our prior knowledge can be expressed by the following extended linear model in which we consider e en as additional parameters to be estimated: A I .. . I (y 0) (0 1
The asymmetry between e O and e, (i > 0) is only an apparent one; we shall show below that the estimator is of the same form for all i >_ 0. The least-squares estimators p and i, (i > 0) are given by the well-known formula which is valid for arbitrary linear models :
Iy~
p~
it I
-N-
10
0 with (ATGOA G OA G OA M=
and
N
ATGO Go Go + GZ
... ... ...
Go Go
Go
Go
...
Go + Gn )
GOA ATG, Go Go
AT GO
ATGO Go + G1 Go
o G,
.. . .. .
0 p
p
.. .
Gn I
Here we use the fact that the e, are independent . p and i, are the linear unbiassed estimators with minimum variance, and M- ' is their covariance matrix [1, p. 871. The matrix C = M - ' is symmetric and consists of the following submatrices :
Permanent address: Nikolsdorfer Gasse 18, A-1050 Wien, Austria. 0168-9002/93/$06 .00 (0 1993 - Elsevier Science Publishers B.V . All rights reserved
318
R Friihwtrth / Estimation of error components
where
The symmetry of the estimator is obvious. Therefore the previous results hold also for i = 0.
Ccm = (ATGA) Cot = -(AT GA) -I ATGV, (t, c,1 =8 ,~V, - V,G BVj with
(i > 0), > 0),
2. Application to track fitting with multiple scattering
var(i, - t,) =V,'-V,'GBV,',
The estimation of track parameters in the presence of multiple scattering can either be done in the framework of a regression model or in the framework of a state-space model. The first approach leads to the classical least-squares estimator ("global fit"), the second one to the Kalman filter/smoother. Both methods have been discussed extensively in the literature [2, and references therein] . In particular, it has been shown that in a linear model the smoother yields minimum variance estimates of the actual track parameters at every point of measurement (and possibly at other points along the trajectory). In contrast, a global fit yields an estimate of the parameter vector at a single fixed reference surface . If there is multiple scattering, the extrapolation of this global estimate does not coincide with the actual scattered track. In case of strong multiple scattering, the extrapolated track may be very far indeed from the actual track, thus decreasing the chances of finding outliers or kinks [3]. Hence we may ask whether estimates of the actual track parameters like those produced by the smoother can also be computed in a global fit. This can indeed be done, if we apply the multi-error regression model developed in section 1, for the case n = 1 . In fact, the difference of the measured track and the track model is a superposition of detector errors and multiple scattering, and these are in most cases independent, at least locally, i.e . in a small portion of the phase space containing a specific track. Therefore we set up the following correspondence : p: true track parameters at a reference surface; y: vector of measurements, measured track; A: linear of linearized track model; eo: vector of detector measurement errors ; e, : vector of deviations from the true track model due to multiple scattering; Ap : ideal unscattered track in measurement space; tl : Ap+e,, actual scattered ("true") track in measurement space. We are interested in an estimate of the actual track t, . We have shown in section 1 that i, is the linear minimum variance estimator of the actual track t, :
vary-i,) =V,'GBV,' .
il = (I - VOGB)y,
n
V = r_
f=U
VI,
G=V - ', B = I - A(ATGA) -' ATG G B = BTG = GB = BTGB = G - GA(AGA) -I ATG. It can readily be verified that indeed C - M = I. Note that BA =ATBT =O, VG B = B,
GBV = BT,
B'= B. Now we can write down the estimated parameter vector p and the estimated error components i, (i > 0) : (ATGA)-IATGy,
P=
i, = VG By . The variances and covariances of these estimates are the respective submatrices of C: var( fi -p) = (ATGA) -', cov(p-p, ë,-e,) = -(ATGA)-IATGV, var(é, -e,) = V, - V,GBV cov(é,-e é,-eJ = -V,GBV1 . In section 2 we shall need a few more quantities and their covariance matrices which we compute by linear error propagation : t,=Ap+e t,=Ap+é,=(I-V,'GB)y, V,'=V-V, ; var(é,) = VG V var(i,-Ap) =A(ATGA) -I AT +VGBV
Finally, we obtain an estimate of eo in the following way: ë o =y-Ap-
- . ..
-en = VoGBy.
var( i, -Ap) =A(ATGA) -I AT +V,G BV,, var( i,-t,) =Vo-VOGBV,
R. Frühwirth / Estimation of error components
In practice it is frequently useful to check this estimator and its covariance matrix by computation of the "pulls" q, i.e . the studentized residuals r: r=y - ii = VOGBy, var(r) = VOGBV0,
qk = rkl var( rk ) . All q k should have mean 0 and a standard deviation of 1.
The equivalence of the estimate il with the one produced by the smoother can be proved by a simple statistical argument : by virtue of the Gauss-Markov theorem, both estimators are minimum variance estimators in the same linear model. Therefore they are identical almost certainly (with probability 1), up to a linear transformation [1, p. 18]. It should be noted, however, that the smoother yields estimates in track parameter space (state space), whereas the above algorithm yields an estimate in measurement space. In many cases, however, a detector does not measure the full parameter vector, but only a subset . Therefore an estimate in parameter space usually can be projected into measurement space, but the parameter vector
319
cannot always be computed from the measurement. If, as for instance in a time projection chamber, the basic measurements of the track are space points, the estimator il yields only the position of the actual track in the points of measurements, but neither the track direction nor the curvature is immediately available. In the filter/smoother, on the other hand, the full estimated parameter vector is available everywhere . Only if the basic measurements are genuine track elements, comprising position, direction, and momentum, the global fit and the filter/smoother give strictly equivalent results. References [1] M.G. Kendall and A. Stuart, The Advanced Theory of Statistics, vol . 2, 2nd ed . (Griffin, London, 1967). [2] M. Regler (ed.), Data Analysis Techniques for High-energy Physics Experiments (Cambridge University Press, Cambridge, 1990) . [3] R. Friihwirth, Application of Filter Methods to the Reconstruction of Tracks and Vertices in Events of Experimental High Energy Physics, HEPHY-PUB 516/88, Vienna (1988) unpublished .
Nuclear Instruments and Methods in Physics Research A324 (1993) 317-319 North-Holland
Estimation of error components in a multi-error linear regression model, with an application to track fitting Rudolf Frühwirth '
Institut für Hochenergiephysik der Österreichischen Akademie der Wissenschaften, Vienna, Austria
Received 20 July 1992 We present an estimation procedure of the error components in a linear regression model with multiple independent stochastic error contributions. After solving the general problem, we apply the results to the estimation of the actual trajectory in track fitting with multiple scattering . 1. A multi-error linear regression model Let us consider a linear regression model in which the stochastic deviation from the linear relationship is a superposition of n + 1 independent random variables: y=Ap+e o + - . . + En . Here y denotes the vector of observations, A the matrix of the known linear model, and p a vector of unknown parameters . We assume that dim(p) is less than dim(y), and that the rank of A is equal to dim(p). In section 2 we shall interpret y as the measurement vector belonging to a track and p as the track parameters to be estimated. The measurement errors of the detector will play the role of eo , and the stochastic deviations due to multiple scattering the role of E, . For a given set of observations y, we wish to estimate the unknown parameters p and the actual values of e i = 0, - - - , n . This can be done if there is prior knowledge on the distribution of the e, . We assume that we know the first two moments (expectation and covariance matrix) of all e, : E(e,)=0, var(e,)=V,=Gl', i=0,---,n . The assumption that E(e,) = 0 implies no loss of generality. Let us denote the actual values of e, as e, . Then our prior knowledge can be expressed by the following extended linear model in which we consider e en as additional parameters to be estimated: A I .. . I (y 0) (0 1
The asymmetry between e O and e, (i > 0) is only an apparent one; we shall show below that the estimator is of the same form for all i >_ 0. The least-squares estimators p and i, (i > 0) are given by the well-known formula which is valid for arbitrary linear models :
Iy~
p~
it I
-N-
10
0 with (ATGOA G OA G OA M=
and
N
ATGO Go Go + GZ
... ... ...
Go Go
Go
Go
...
Go + Gn )
GOA ATG, Go Go
AT GO
ATGO Go + G1 Go
o G,
.. . .. .
0 p
p
.. .
Gn I
Here we use the fact that the e, are independent . p and i, are the linear unbiassed estimators with minimum variance, and M- ' is their covariance matrix [1, p. 871. The matrix C = M - ' is symmetric and consists of the following submatrices :
Permanent address: Nikolsdorfer Gasse 18, A-1050 Wien, Austria. 0168-9002/93/$06 .00 (0 1993 - Elsevier Science Publishers B.V . All rights reserved
318
R Friihwtrth / Estimation of error components
where
The symmetry of the estimator is obvious. Therefore the previous results hold also for i = 0.
Ccm = (ATGA) Cot = -(AT GA) -I ATGV, (t, c,1 =8 ,~V, - V,G BVj with
(i > 0), > 0),
2. Application to track fitting with multiple scattering
var(i, - t,) =V,'-V,'GBV,',
The estimation of track parameters in the presence of multiple scattering can either be done in the framework of a regression model or in the framework of a state-space model. The first approach leads to the classical least-squares estimator ("global fit"), the second one to the Kalman filter/smoother. Both methods have been discussed extensively in the literature [2, and references therein] . In particular, it has been shown that in a linear model the smoother yields minimum variance estimates of the actual track parameters at every point of measurement (and possibly at other points along the trajectory). In contrast, a global fit yields an estimate of the parameter vector at a single fixed reference surface . If there is multiple scattering, the extrapolation of this global estimate does not coincide with the actual scattered track. In case of strong multiple scattering, the extrapolated track may be very far indeed from the actual track, thus decreasing the chances of finding outliers or kinks [3]. Hence we may ask whether estimates of the actual track parameters like those produced by the smoother can also be computed in a global fit. This can indeed be done, if we apply the multi-error regression model developed in section 1, for the case n = 1 . In fact, the difference of the measured track and the track model is a superposition of detector errors and multiple scattering, and these are in most cases independent, at least locally, i.e . in a small portion of the phase space containing a specific track. Therefore we set up the following correspondence : p: true track parameters at a reference surface; y: vector of measurements, measured track; A: linear of linearized track model; eo: vector of detector measurement errors ; e, : vector of deviations from the true track model due to multiple scattering; Ap : ideal unscattered track in measurement space; tl : Ap+e,, actual scattered ("true") track in measurement space. We are interested in an estimate of the actual track t, . We have shown in section 1 that i, is the linear minimum variance estimator of the actual track t, :
vary-i,) =V,'GBV,' .
il = (I - VOGB)y,
n
V = r_
f=U
VI,
G=V - ', B = I - A(ATGA) -' ATG G B = BTG = GB = BTGB = G - GA(AGA) -I ATG. It can readily be verified that indeed C - M = I. Note that BA =ATBT =O, VG B = B,
GBV = BT,
B'= B. Now we can write down the estimated parameter vector p and the estimated error components i, (i > 0) : (ATGA)-IATGy,
P=
i, = VG By . The variances and covariances of these estimates are the respective submatrices of C: var( fi -p) = (ATGA) -', cov(p-p, ë,-e,) = -(ATGA)-IATGV, var(é, -e,) = V, - V,GBV cov(é,-e é,-eJ = -V,GBV1 . In section 2 we shall need a few more quantities and their covariance matrices which we compute by linear error propagation : t,=Ap+e t,=Ap+é,=(I-V,'GB)y, V,'=V-V, ; var(é,) = VG V var(i,-Ap) =A(ATGA) -I AT +VGBV
Finally, we obtain an estimate of eo in the following way: ë o =y-Ap-
- . ..
-en = VoGBy.
var( i, -Ap) =A(ATGA) -I AT +V,G BV,, var( i,-t,) =Vo-VOGBV,
R. Frühwirth / Estimation of error components
In practice it is frequently useful to check this estimator and its covariance matrix by computation of the "pulls" q, i.e . the studentized residuals r: r=y - ii = VOGBy, var(r) = VOGBV0,
qk = rkl var( rk ) . All q k should have mean 0 and a standard deviation of 1.
The equivalence of the estimate il with the one produced by the smoother can be proved by a simple statistical argument : by virtue of the Gauss-Markov theorem, both estimators are minimum variance estimators in the same linear model. Therefore they are identical almost certainly (with probability 1), up to a linear transformation [1, p. 18]. It should be noted, however, that the smoother yields estimates in track parameter space (state space), whereas the above algorithm yields an estimate in measurement space. In many cases, however, a detector does not measure the full parameter vector, but only a subset . Therefore an estimate in parameter space usually can be projected into measurement space, but the parameter vector
319
cannot always be computed from the measurement. If, as for instance in a time projection chamber, the basic measurements of the track are space points, the estimator il yields only the position of the actual track in the points of measurements, but neither the track direction nor the curvature is immediately available. In the filter/smoother, on the other hand, the full estimated parameter vector is available everywhere . Only if the basic measurements are genuine track elements, comprising position, direction, and momentum, the global fit and the filter/smoother give strictly equivalent results. References [1] M.G. Kendall and A. Stuart, The Advanced Theory of Statistics, vol . 2, 2nd ed . (Griffin, London, 1967). [2] M. Regler (ed.), Data Analysis Techniques for High-energy Physics Experiments (Cambridge University Press, Cambridge, 1990) . [3] R. Friihwirth, Application of Filter Methods to the Reconstruction of Tracks and Vertices in Events of Experimental High Energy Physics, HEPHY-PUB 516/88, Vienna (1988) unpublished .