Regression with Gaussian Mixture ModelsApplied to Track Fitting

Abstract

This note describes the application of Gaussian mixture regression to track fitting with a Gaussian mixture model of the position errors. The mixture model is assumed to have two components with identical component means. Under the premise that the association of each measurement to a specific mixture component is known, the Gaussian mixture regression is shown to have consistently better resolution than weighted linear regression with equivalent homoskedastic errors. The improvement that can be achieved is systematically investigated over a wide range of mixture distributions. The results confirm that with constant homoskedastic variance the gain is larger for larger mixture weight of the narrow component and for smaller ratio of the width of the narrow component and the width of the wide component.

1. Introduction

In the data analysis of high-energy physics experiments, and particularly in track fitting, Gaussian models of measurements errors or Gaussian models of stochastic processes such as multiple Coulomb scattering and energy loss of electrons by bremsstrahlung may turn out to be too simplistic or downright invalid. In these circumstances Gaussian mixture models (GMMs) can be used in the data analysis to model outliers or tails in position measurements of tracks [1,2], to model the tails of the multiple Coulomb scattering distribution [3,4], or to approximate the distribution of the energy loss of electrons by bremsstrahlung [5]. In these applications it is not known to which component an observation or a physical process corresponds, and it is up to the track fit via the Deterministic Annealing Filter or the Gaussian-sum Filter to determine the association probabilities [6,7,8,9]. The standard method for concurrent maximum-likelihood estimation of the regression parameters and the association probabilities is the EM algorithm [10,11].

In a typical tracking detector such as a silicon strip tracker, the spatial resolution of the measured position is not necessarily uniform across the sensors and depends on the crossing angle of the track and on the type of the cluster created by the particle. By a careful calibration of the sensors in the tracking detector a GMM of the position measurements can be formulated, see [12,13,14]. In such a model the component from which an observation is drawn is known. This can significantly improve the precision of the estimated regression/track parameters with respect to weighted least-squares, as shown below.

This introduction is followed by a recapitulation of weighted regression in linear and nonlinear Gaussian models in Section 2, both for the homoskedastic and for the heteroskedastic case. Section 3 introduces the concept of Gaussian mixture regression (GMR) and derives some elementary results on the covariance matrices of the estimated regression parameters. Section 4 presents the results of simulation studies that apply GMR to track fitting in a simplified tracking detector. The gain in precision that can be obtained by using a GMM rather than a simple Gaussian model strongly depends on the shape of the mixture distribution of the position errors. This shape can be described concisely by two characteristic quantities: the weight p of the wider component, and the ratio r of the smaller over the larger standard deviation. The study covers a wide range of p and r and shows the gain as a function of the number of measurements within a fixed tracker volume. It also investigates the effect of multiple Coulomb scattering on the gain in precision. Finally, Section 5 gives a summary of the results and discusses them in context of the related work in [12,13,14].

2. Linear Regression in Gaussian Models

This section briefly recapitulates the properties of linear Gaussian regression models and the associated estimation procedures.

2.1. Linear Gaussian Regression Models

The general linear Gaussian regression model (LGRM) has the following form:

$$\begin{array}{c}\mathit{y}=\mathit{X}\mathit{\theta }+\mathit{c}+\mathit{ϵ},\mathrm{with}\mathit{V}=\mathrm{Var}\left[\mathit{ϵ}\right]\mathrm{and}\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathit{V}),\end{array}$$

(1)

where $\mathit{y}$ is the vector of observations, $\mathit{\theta }$ is the vector of unknown regression parameters, $\mathit{X}$ is the known model matrix, $\mathit{c}$ is a known vector of constants, and $\mathit{ϵ}$ is the vector of observation errors, distributed according to the multivariate normal distribution $\mathcal{N}(\mathbf{0},\mathit{V})$ with mean $\mathbf{0}$ and covariance matrix $\mathit{V}$. In the following it is assumed that $\mathsf{dim}\left(\mathit{y}\right)\ge \mathsf{dim}\left(\mathit{\theta }\right)$ and $\mathsf{rank}\left(\mathit{X}\right)=\mathsf{dim}\left(\mathit{\theta }\right)$.

If the covariance matrix $\mathit{V}$ is diagonal, i.e., the observation errors are independent, two cases are distinguished in the literature. The errors in Equation (1) and the model are called homoskedastic, if $\mathit{V}$ is a multiple of the identity matrix:

$$\begin{array}{c}\mathit{V}={\sigma }^2·\mathit{I}.\end{array}$$

(2)

Otherwise, the errors in Equation (1) and the model are called heteroskedastic. If $\mathit{V}$ is not diagonal, the observation errors are correlated, in which case the the model is called homoskedastic if all diagonal elements are identical, and heteroskedastic otherwise.

2.2. Estimation, Fisher Information Matrix and Efficiency

In least-squares (LS) estimation, the regression parameters $\mathit{\theta }$ are estimated my minimizing the following quadratic form with respect to $\mathit{\theta }$:

$$\begin{array}{c}\mathcal{S}\left(\mathit{\theta }\right)=(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }){}^{\mathrm{T}}\mathit{G}(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }),\mathrm{with}\mathit{G}={\mathit{V}}^{-1}.\end{array}$$

(3)

Setting the gradient of $\mathcal{S}\left(\mathit{\theta }\right)$ to zero gives the LS estimate $\widehat{\mathit{\theta }}$ of $\mathit{\theta }$, and linear error propagation yields its covariance matrix $\widehat{\mathit{C}}$:

$$\begin{array}{c}\widehat{\mathit{\theta }}={\left(\mathit{X}{}^{\mathrm{T}}\mathit{G}\mathit{X}\right)}^{-1}\mathit{X}{}^{\mathrm{T}}\mathit{G}(\mathit{y}-\mathit{c}),\widehat{\mathit{C}}=\mathsf{Var}[\widehat{\mathit{\theta }}]={\left(\mathit{X}{}^{\mathrm{T}}\mathit{G}\mathit{X}\right)}^{-1}.\end{array}$$

(4)

It is easy to see that $\widehat{\mathit{\theta }}$ is an unbiased estimate of the unknown parameter $\mathit{\theta }$.

In the LGRM, the LS estimate is equal to the maximum-likelihood (ML) estimate. The likelihood function $\mathcal{L}$($\mathit{\theta }$) and its logarithm $\ell \left(\mathit{\theta }\right)$ are given by:

$$\mathcal{L}\left(\mathit{\theta }\right)=C·exp\left[-\frac{1}{2}(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }){}^{\mathrm{T}}\mathit{G}(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta })\right],$$

(5)

$$\ell \left(\mathit{\theta }\right)=log\left(C\right)-\frac{1}{2}(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }){}^{\mathrm{T}}\mathit{G}(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }),C\in \mathbb{R}.$$

(6)

The gradient $\nabla \ell \left(\mathit{\theta }\right)$ is equal to:

$$\begin{array}{c}\nabla \ell \left(\mathit{\theta }\right)=\frac{\partial \ell \left(\mathit{\theta }\right)}{\partial \mathit{\theta }}=(\mathit{y}-\mathit{c}-\mathit{X}\mathit{\theta }){}^{\mathrm{T}}\mathit{G}\mathit{X},\end{array}$$

(7)

and its matrix of second derivatives, the Hesse matrix $\mathit{H}\left(\mathit{\theta }\right)$ is equal to:

$$\begin{array}{c}\mathit{H}\left(\mathit{\theta }\right)=\frac{{\partial }^2\ell \left(\mathit{\theta }\right)}{\partial {\mathit{\theta }}^2}=-\mathit{X}{}^{\mathrm{T}}\mathit{G}\mathit{X}.\end{array}$$

(8)

Setting the gradient to zero yields the ML estimate, which is identical to the LS estimate.

Under mild regularity conditions, the covariance matrix of an unbiased estimator cannot be smaller than the inverse of the Fisher information matrix; this is the famous Cramér-Rao inequality [15]. The (expected) Fisher information matrix $\mathcal{I}\left(\mathit{\theta }\right)$ of the parameters $\mathit{\theta }$ is defined as the negative expectation of the Hesse matrix of the log-likelihood function:

$$\begin{array}{c}\mathcal{I}\left(\mathit{\theta }\right)=-\mathsf{E}\left[\left(\mathit{\theta }\right)\right]=\mathit{X}{}^{\mathrm{T}}\mathit{G}\mathit{X}.\end{array}$$

(9)

Equation (4) shows that in the LGRM the covariance matrix of the LS/ML estimate is equal to the inverse of the Fisher information matrix. The LS/ML estimator is therefore efficient: no other unbiased estimator of $\mathit{\theta }$ can have a smaller covariance matrix.

2.3. Nonlinear Gaussian Regression Models

In most applications to track fitting, the track model that describes the dependence of the observed track positions on the track parameters is nonlinear. The general nonlinear regression Gaussian model has the following form:

$$\begin{array}{c}\mathit{y}=\mathit{f}\left(\mathit{\theta }\right)+\mathit{ϵ},\mathrm{with}\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathit{V}),\end{array}$$

(10)

where $\mathit{f}\left(\mathit{\theta }\right)$ is a nonlinear smooth function. The estimation of $\mathit{\theta }$ usually proceeds via Taylor expansion of the model function to first order:

$$\mathit{y}\approx \mathit{f}\left({\mathit{\theta }}_0\right)+\mathit{X}(\mathit{\theta }-{\mathit{\theta }}_0)+\mathit{ϵ}=\mathit{X}\mathit{\theta }+\mathit{c}+\mathit{ϵ},\mathrm{with}$$

(11)

$$\mathit{X}={\left.\frac{\partial \mathit{f}}{\partial \mathit{\theta }}\right|}_{{\mathit{\theta }}_0},\mathit{c}=\mathit{f}\left({\mathit{\theta }}_0\right)-\mathit{X}{\mathit{\theta }}_0,$$

(12)

where ${\mathit{\theta }}_0$ is a suitable expansion point and the Jacobian $\mathit{X}$ is evaluated at ${\mathit{\theta }}_0$. The estimate ${\mathit{\theta }}_1$ is obtained according to Equation (4). The model function is then re-expanded at ${\mathit{\theta }}_1$, and the procedure is iterated until convergence.

3. Linear Regression in Gaussian Mixture Models

If the observation errors do not follow a normal distribution, the Gaussian linear model can be generalized to a GMM, which is more flexible, but still retains some of the computational advantages of the Gaussian model. Applications of Gaussian mixtures can be traced back to Karl Pearson more than a century ago [16]. In view of the applications in Section 4, the following discussion is restricted to Gaussian mixtures with two components with identical means. The probability density function (PDF) $g\left(y\right)$ of such a mixture has the following form:

$$\begin{array}{c}g\left(y\right)=p·\phi (y;\mu ,{\sigma }_1)+(1-p)·\phi (y;\mu ,{\sigma }_0),0<p<1,{\sigma }_0<{\sigma }_1,\end{array}$$

(13)

where $\phi (y;\mu ,\sigma )$ denotes the normal PDF with mean $\mu$ and variance ${\sigma }^2$. The Gaussian mixture (GM) with this PDF is denoted by $\mathcal{GM}(p,\mu ,{\sigma }_1,{\sigma }_0)$. The component with index 0 is called the narrow component, the component with index 1 is called the wide component in the following.

3.1. Homoskedastic Mixture Models

We first consider the following linear GMM:

$$\begin{array}{c}\mathit{y}=\mathit{X}\mathit{\theta }+\mathit{c}+\mathit{ϵ},\mathrm{with}{ϵ}_i\sim \mathcal{GM}(p,0,{\sigma }_1,{\sigma }_0),i=1,\dots ,n,\end{array}$$

(14)

where $n=\mathsf{dim}\left(\mathit{y}\right)$. The variance of ${ϵ}_i$ is equal to:

$$\begin{array}{c}\mathsf{var}\left[{ϵ}_i\right]={\sigma }_{\mathrm{tot}}^2=p·{\sigma }_1^2+(1-p)·{\sigma }_0^2,i=1,\dots ,n,\end{array}$$

(15)

and the joint covariance matrix of $\mathit{ϵ}$ is equal to

$$\begin{array}{c}\mathit{V}=\mathsf{Var}\left[\mathit{ϵ}\right]={\sigma }_{\mathrm{tot}}^2·\mathit{I}.\end{array}$$

(16)

As the total variance is the same for all i, the GMM is called homoskedastic. By using just the total variances, it can be approximated by a homoskedastic LGRM, and the estimation of the regression parameters proceeds as in Equation (4), which can be simplified to:

$$\begin{array}{c}\check{\mathit{\theta }}={\left(\mathit{X}{}^{\mathrm{T}}\mathit{X}\right)}^{-1}\mathit{X}{}^{\mathrm{T}}(\mathit{y}-\mathit{c}),\check{\mathit{C}}={\sigma }_{\mathrm{tot}}^2·{\left(\mathit{X}{}^{\mathrm{T}}\mathit{X}\right)}^{-1}.\end{array}$$

(17)

If the association of the observations to the components is known, it can be encoded in a binary vector $\mathit{z}\in Z={\{0,1\}}^n$, where $z_i=j$ indicates the component with variance ${\sigma }_j^2$, for $i=1,\dots ,n$. In the case of independent draws from the mixture PDF in Equation (13), the $z_i$ are independent Bernoulli variables with $P\left(1\right)=p,P\left(0\right)=1-p$. The probability of $\mathit{z}=(z_1,\dots ,z_n)$ is therefore equal to:

$$\begin{array}{c}P\left(\mathit{z}\right)=\prod _{i=1}^n{p}^{z_i}{(1-p)}^{1-z_i}=p^k{(1-p)}^{n-k},\mathrm{with}k=\sum _{i=1}^n{z}_i.\end{array}$$

(18)

Conditional on the known associations $\mathit{z}$, the model can be interpreted as a heteroskedastic LGRM with a diagonal covariance matrix ${\mathit{V}}_{\mathit{z}}$:

$$\begin{array}{c}{\mathit{V}}_{\mathit{z}}=\mathsf{Var}\left[\mathit{ϵ}\right|\mathit{z}]=\mathsf{diag}({\sigma }_{z_1}^2,\dots ,{\sigma }_{z_n}^2),\mathrm{and}{\mathit{G}}_{\mathit{z}}={\mathit{V}}_{\mathit{z}}^{-1}.\end{array}$$

(19)

It follows from the binomial theorem that:

$$\begin{array}{c}\mathit{V}=\sum _{\mathit{z}\in Z}P\left(\mathit{z}\right)·{\mathit{V}}_{\mathit{z}}.\end{array}$$

(20)

The GMR is thus a mixture of heteroskedastic Gaussian linear regressions with mixture weights $P\left(\mathit{z}\right),\mathit{z}\in Z$. In each of the heteroskedastic LGRMs the regression parameters are estimated by:

$$\begin{array}{c}{\widehat{\mathit{\theta }}}_{\mathit{z}}={\left(\mathit{X}{}^{\mathrm{T}}{\mathit{G}}_{\mathit{z}}\mathit{X}\right)}^{-1}\mathit{X}{}^{\mathrm{T}}{\mathit{G}}_{\mathit{z}}(\mathit{y}-\mathit{c}),{\widehat{\mathit{C}}}_{\mathit{z}}={\left(\mathit{X}{}^{\mathrm{T}}{\mathit{G}}_{\mathit{z}}\mathit{X}\right)}^{-1}.\end{array}$$

(21)

Conditional on $\mathit{z}$, the estimator is normally distributed with mean $\mathit{\theta }$ and covariance matrix ${\widehat{\mathit{C}}}_{\mathit{z}}$, and therefore unbiased and efficient. The unconditional distribution of $\widehat{\mathit{\theta }}$ is obtained by summing over all $\mathit{z}$ with weights $P\left(\mathit{z}\right)$. Its density therefore reads:

$$\begin{array}{c}f(\widehat{\mathit{\theta }})=\sum _{\mathit{z}\in Z}P\left(\mathit{z}\right)·\phi (\widehat{\mathit{\theta }};\mathit{\theta },{\widehat{\mathit{C}}}_{\mathit{z}}),\end{array}$$

(22)

where $\phi (\widehat{\mathit{\theta }};\mathit{\theta },{\widehat{\mathit{C}}}_{\mathit{z}})$ is the multivariate normal density of $\widehat{\mathit{\theta }}$ with mean $\mathit{\theta }$ and covariance matrix ${\widehat{\mathit{C}}}_{\mathit{z}}$. $\widehat{\mathit{\theta }}$ is therefore unbiased, and its covariance matrix is obtained by summing over the $2^n$ possible values of $\mathit{z}$ with weights $P\left(\mathit{z}\right)$:

$$\begin{array}{c}\widehat{\mathit{C}}=\sum _{\mathit{z}\in Z}P\left(\mathit{z}\right)·{\widehat{\mathit{C}}}_{\mathit{z}}.\end{array}$$

(23)

The unconditional estimator has minimum variance among all estimators that are conditionally unbiased for all $\mathit{z}\in Z$, as shown by the following theorem.

Theorem 1.

Let $\tilde{\mathit{\theta }}$ be an estimator of $\mathit{\theta }$ with $\mathsf{E}[\tilde{\mathit{\theta }}|\mathit{z}]=\mathit{\theta }$ for all $\mathit{z}\in Z$. Then its covariance matrix $\tilde{\mathit{C}}$ is not smaller than $\widehat{\mathit{C}}$ in the Loewner partial order:

$$\begin{array}{c}\widehat{\mathit{C}}\le \tilde{\mathit{C}}.\end{array}$$

(24)

Proof. 

As $\mathsf{E}[\tilde{\mathit{\theta }}|\mathit{z}]=\mathit{\theta }$, $\tilde{\mathit{C}}$ can be written in the following form:

$$\begin{array}{c}\tilde{\mathit{C}}=\sum _{\mathit{z}\in Z}P\left(\mathit{z}\right)·{\tilde{\mathit{C}}}_{\mathit{z}},\mathrm{with}{\tilde{\mathit{C}}}_{\mathit{z}}=\mathsf{Var}[\tilde{\mathit{\theta }}|\mathit{z}].\end{array}$$

(25)

As ${\widehat{\mathit{\theta }}}_{\mathit{z}}$ is efficient conditional on $\mathit{z}$ and $\tilde{\mathit{\theta }}|\mathit{z}$ is conditionally unbiased, it follows that

$$\begin{array}{c}{\widehat{\mathit{C}}}_{\mathit{z}}=\mathsf{Var}[\widehat{\mathit{\theta }}|\mathit{z}]\le \mathsf{Var}[\tilde{\mathit{\theta }}|\mathit{z}]={\tilde{\mathit{C}}}_{\mathit{z}},\mathit{z}\in Z.\end{array}$$

(26)

Summing over all $\mathit{z}$ with weights $P\left(\mathit{z}\right)$ proves the theorem. □

Collorary. 

$\widehat{\mathit{C}}\le \check{\mathit{C}}$, with equality if $p\in \{0,1\}$.

Proof. 

As $\mathsf{E}[\check{\mathit{\theta }}|\mathit{z}]=\mathit{\theta }$ for all $\mathit{z}$, the estimator $\check{\mathit{\theta }}$ fulfills the premise of the theorem. Its conditional covariance matrix is obtained by linear error propagation:

$$\begin{array}{c}{\check{\mathit{C}}}_{\mathit{z}}=\mathsf{Var}[\check{\mathit{\theta }}|\mathit{z}]={\left(\mathit{X}{}^{\mathrm{T}}\mathit{X}\right)}^{-1}\mathit{X}{}^{\mathrm{T}}{\mathit{V}}_{\mathit{z}}\mathit{X}{\left(\mathit{X}{}^{\mathrm{T}}\mathit{X}\right)}^{-1}.\end{array}$$

(27)

As $\widehat{\mathit{\theta }}|\mathit{z}$ is efficient, ${\widehat{\mathit{C}}}_{\mathit{z}}\le {\check{\mathit{C}}}_{\mathit{z}}$ and $\widehat{\mathit{C}}\le \check{\mathit{C}}$. If $p\in \{0,1\}$, ${\mathit{V}}_{\mathit{z}}$ is homoskedastic and ${\widehat{\mathit{C}}}_{\mathit{z}}={\check{\mathit{C}}}_{\mathit{z}}={\sigma }_{\mathrm{tot}}^2·{\left(\mathit{X}{}^{\mathrm{T}}\mathit{X}\right)}^{-1}$. □

The effective improvement in precision of the heteroskedastic estimator $\widehat{\mathit{\theta }}$ w.r.t.regarding to the homoskedastic estimator $\check{\mathit{\theta }}$ is difficult to quantify by simple formulas, but easy to assess by simulation studies. The results of three such studies in the context of track fitting are presented in Section 4.1, Section 4.2 and Section 4.3, respectively.

3.2. Heteroskedastic Mixture Models

The GMM in Equation (14) can be further generalized to the heteroskedastic case, where the mixture parameters depend on i:

$$\begin{array}{c}\mathit{y}=\mathit{X}\mathit{\theta }+\mathit{c}+\mathit{ϵ},\mathrm{with}{ϵ}_i\sim \mathcal{GM}(p_i,0,{\sigma }_{i,1},{\sigma }_{i,0}),i=1,\dots ,n.\end{array}$$

(28)

The joint covariance matrix of $\mathit{ϵ}$ is now diagonal with

$$\mathit{V}=\mathsf{Var}\left[\mathit{ϵ}\right]=\mathsf{diag}({\sigma }_{1,\mathrm{tot}}^2,\dots ,{\sigma }_{n,\mathrm{tot}}^2),\mathrm{with}$$

(29)

$${\sigma }_{i,\mathrm{tot}}^2=p_i·{\sigma }_{i,1}^2+(1-p_i)·{\sigma }_{i,0}^2,i=1,\dots ,n.$$

(30)

The approximating LGRM that uses just the total variances is now heteroskedastic with covariance matrix $\mathit{V}$. The corresponding estimator $\check{\mathit{\theta }}$ with its covariance matrix $\check{\mathit{C}}$ is computed as in Equation (4).

Conditional on the associations coded in $\mathit{z}$, the covariance matrix ${\mathit{V}}_{\mathit{z}}$ now reads:

$$\begin{array}{c}{\mathit{V}}_{\mathit{z}}=\mathsf{Var}\left[\mathit{ϵ}\right|\mathit{z}]=\mathsf{diag}({\sigma }_{1,z_1}^2,\dots ,{\sigma }_{n,z_n}^2).\end{array}$$

(31)

As in Section 3.1, the unconditional covariance matrix $\widehat{\mathit{C}}$ is the weighted sum of the conditional covariance matrices ${\widehat{\mathit{C}}}_{\mathit{z}}$:

$$\begin{array}{c}\widehat{\mathit{C}}=\sum _{\mathit{z}\in Z}P\left(\mathit{z}\right)·{\widehat{\mathit{C}}}_{\mathit{z}},\mathrm{with}P\left(\mathit{z}\right)=\prod _{i=1}^n{p}_i^{z_i}{(1-p_i)}^{1-z_i}.\end{array}$$

(32)

It can be proved as before that $\widehat{\mathit{C}}\le \check{\mathit{C}}$.

In the application to track fitting, heteroskedastic mixture models can be useful in non-uniform tracking detectors, in which the mixture model obtained by the calibration depends on the layer or on the sensor or on the possible cluster types.

A nonlinear GMM can be approximated by a linear GMM in the same way as a nonlinear Gaussian model can be approximated by a LGRM, see Section 2.3.

4. Track Fitting with Gaussian Mixture Models

4.1. Simulation Study with Straight Tracks

The setting of the first simulation study is a toy detector that consists of n equally spaced sensor planes in the interval from $0\mathrm{cm}$ to $115\mathrm{cm}$. The track model is a straight line in the $(x,y)$-plane. No material effects are simulated. The track position is recorded in each sensor plane, with a position error generated from a Gaussian mixture, see Equation (13). The total variance in each layer is equal to ${\sigma }_{\mathrm{tot}}^2={\left(0.005\mathrm{cm}\right)}^2$, corresponding to a standard deviation of ${\sigma }_{\mathrm{tot}}=0.005\mathrm{cm}$. As the main objective of the simulation study is the comparison of the homoskedastic estimator of Equation (17) with the heteroskedastic estimator of Equation (21), the actual value of ${\sigma }_{\mathrm{tot}}^2$ is irrelevant.

${10}^5$ tracks were generated for each of 20 combinations of p and $r={\sigma }_0/{\sigma }_1$: $p=0.6,0.7,0.8,0.9;$$r=0.1,0.2,0.3,0.4,0.5$. The corresponding values of ${\sigma }_1$ and ${\sigma }_0$ are collected in

Table 1. Values of ${\sigma }_1$ and ${\sigma }_0$ (in µm) for $p=0.6,0.7,0.8,0.9$ and $r=0.1,0.2,0.3,0.4,0.5$.

	p = 0.6		p = 0.7		p = 0.8		p = 0.9
	${\sigma }_1$	${\sigma }_0$	${\sigma }_1$	${\sigma }_0$	${\sigma }_1$	${\sigma }_0$	${\sigma }_1$	${\sigma }_0$
r = 0.1	64.3	6.4	59.6	6.0	55.8	5.6	52.7	5.3
r = 0.2	63.7	12.7	59.3	11.9	55.6	11.1	52.6	10.5
r = 0.3	62.7	18.8	58.6	17.6	55.3	16.6	52.4	15.7
r = 0.4	61.4	24.5	57.8	23.1	54.8	21.9	52.2	20.9
r = 0.5	59.7	29.9	56.8	28.4	54.2	27.1	52.0	26.0

. The number n of planes varied from 3 to 20. For each pair $(p,r)$ and each value of n the exact covariance matrix was computed according to Equation (23).

The track parameters were chosen as the intercept d and the slope k. All tracks were fitted with the homoskedastic and with the heteroskedastic error model. Figure 1 shows the estimated standard deviation (rms) of both fits, along with the exact values, for the intercept d and the slope k, with $p=0.8,r=0.2$ in the range $3\le n\le 20$. The standard deviations of the heteroskedastic fit are not only smaller, but also fall at a faster rate than the ones of the homoskedastic fit, i.e., faster than $1/\sqrt{n}$. The agreement between the estimated and the exact standard deviations is excellent.

Figure 2 shows the relative standard deviation (heteroskedastic rms over homoskedastic rms) for the slope k and all 20 combinations of p and r. The corresponding plot for the intercept d looks very similar. A clear trend can be observed: the ratio gets smaller with decreasing p and with decreasing r. The largest improvement is therefore seen for $p=0.6,r=0.1$, where the probability of very precise measurements ($r=0.1$) is the largest ($1-p=0.4$).

For large values of n the ratio curves approach a constant level, because asymptotically both rms values fall at the same rate, namely $1/\sqrt{n}$, see Figure 3. In the experimentally relevant region $15\le n\le 20$ the ratio is already close to its lower limit in most, but not all, cases. In these cases, the superior performance of the heteroskedastic fit is fully exploited, and adding more measurements yields little further improvement. The exception is the case $p=0.9$, where the limit is reached much later, at $n=200$ in the extreme case.

4.2. Simulation Study with Circular Tracks

The setting of the second simulation study is a toy detector that consists of n equally spaced concentric cylinders with radii in the interval $[R_{\min },R_{\max }]$, with $R_{\min }=10,R_{\max }=100\mathrm{cm}$. The track model is a helix through the origin. Only the parameters in the bending plane ($(x,y)$-plane) are estimated. The track parameters are: $\mathrm{\Phi }$, the azimuth of the track position at $R=R_{\min }$; $\varphi$, the azimuth of the track direction at $R=R_{\min }$; and $\kappa$, the curvature of the circle.

${10}^6$ tracks were generated for each of the same 20 combinations of p and $r={\sigma }_0/{\sigma }_1$ as in Section 4.1. The track position in the bending plane was recorded in each cylinder, with a position error generated from the Gaussian mixtures in Table 1. Each sample of ${10}^6$ tracks consists of 10 subsamples with circles of 10 different radii, namely $\rho /\mathrm{cm}=100,200,\dots ,1000$. The number n of planes varied from 3 to 20.

As the track model is nonlinear, estimation was performed according to the procedure in Section 2.3. All tracks were fitted with the homoskedastic and with the heteroskedastic error model, and for each pair $(p,r)$ and each value of n the exact covariance matrix was computed according to Equation (23). Figure 4 shows the estimated standard deviation (rms) of both fits, along with the exact values, for $\mathrm{\Phi }$, $\varphi$ and $\kappa$, with $p=0.8,r=0.2$ in the range $3\le n\le 20$. The standard deviations of the heteroskedastic fit are not only smaller, but also fall at a faster rate than the ones of the homoskedastic fit, i.e., faster than $1/\sqrt{n}$. The agreement between the estimated and the exact standard deviations is again excellent.

Figure 5, Figure 6 and Figure 7 show the relative standard deviation for the three track parameters and all 20 combinations of p and r. Because of the large correlation between $\varphi$ and $\kappa$, Figure 6 and Figure 7 are virtually identical. The same trend as before can be observed: the ratio gets smaller with decreasing p and with decreasing r. The largest improvement is again seen for $p=0.6,r=0.1$, where the probability of very precise measurements ($r=0.1$) is the largest ($1-p=0.4$). Compared to the results of the straight-line fits, the gain in precision is somewhat smaller.

4.3. Simulation Study with Circular Tracks and Multiple Scattering

The third study investigates the effect of multiple Coulomb scattering (MCS) on the results presented in the previous subsection. To this end, every layer of the toy detector is assumed to have a thickness of 1% of a radiation length. The simulation and the reconstruction of the tracks is modified accordingly, see for example [17].

${10}^5$ tracks were generated for each of four values of the transverse momentum, i.e., $p_{\mathrm{T}}=10,20,50,100\mathrm{GeV}/c$ and for each of the 20 combinations of p and $r={\sigma }_0/{\sigma }_1$ as in Section 4.2. The dip angle of the tracks was set to $\lambda =\pi /4$. The track position in the bending plane was recorded in each cylinder, with a position error generated from the Gaussian mixtures in Table 1 plus a random deviation due to MCS. The number n of planes varied from 3 to 20.

It is to be expected that the benefit provided by the heteroskedastic regression is diminished by the introduction of MCS. This is confirmed by the results in Figure 8, Figure 9, Figure 10 and Figure 11, which show the relative standard deviation of the curvature $\kappa$ for the four values of $p_{\mathrm{T}}$ and the 20 combinations of p and r. Although the gain in precision of the heteroskedastic estimator is rather poor for the lowest value of $p_{\mathrm{T}}$, see Figure 8, it steadily improves with increasing $p_{\mathrm{T}}$ and approaches the optimal level for the highest value of $p_{\mathrm{T}}$, see Figure 7 and Figure 11.

5. Summary and Conclusions

We have evaluated the gain in precision that can be achieved by GMR in track fitting using a mixture model of position errors rather than a simple Gaussian model with an average resolution. The results, both with straight and with circular tracks, show the expected behavior: the gain rises if the mixture weight (proportion) of the narrow component becomes larger, and it rises when the ratio of the width of the narrow component and the width of the wide component becomes smaller. Unsurprisingly, the gain is negatively affected by MCS, depending on the track momentum and on the amount of material.

It is also instructive to look at the results in the light of the findings reported in [12,13,14], in particular the claim that the heteroskedastic regression, which is in fact the GMR described here, achieves a linear growth of the resolution with the number n of detecting layers, instead of the usual $\sqrt{n}$ behavior [13]. A glance at Figure 3 shows that in fact the heteroskedastic standard deviations initially fall faster than the homoskedastic ones, i.e., faster than $1/\sqrt{n}$. Eventually, however, the relative standard deviations level out at a constant ratio, showing that asymptotically both standard deviations obey the expected $1/\sqrt{n}$ mode of falling with n. The observation reported in [13] is clearly a small-sample effect, which—in contrast to most other cases—is very significant in the mixture with $p=0.8,r=0.1$, the one assumed in [13]. Fortunately, it is in the range of n that is relevant for realistic tracking detectors, i.e., somewhere between 5 and 20 for a silicon tracker. It seems therefore absolutely worthwhile to optimize the local reconstruction and error calibration of the position measurements as far as possible.