
%Paper: hep-lat/9305014
%From: seibert@surya11.cern.ch (David Seibert)
%Date: Mon, 17 May 93 11:46:00 +0200


\documentstyle[12pt]{article}
\pagestyle{empty}
\raggedright \raggedbottom
\setlength{\textwidth}{15.8cm}
\setlength{\oddsidemargin}{0pt}
\setlength{\evensidemargin}{0pt}
\setlength{\textheight}{22.7cm}
\setlength{\topmargin}{0.2in}
\setlength{\headheight}{0pt}
\setlength{\headsep}{0pt}
\setlength{\footheight}{0pt}
\setlength{\footskip}{20pt}
\renewcommand{\section}[1]{\refstepcounter{section}
\vspace{24pt}\noindent{\bf\arabic{section}.\quad #1}
\vspace*{12pt}}
\newcommand{\ulsect}[1]{\vspace{18pt}\noindent{\bf #1}
\vspace*{12pt}}
%\renewcommand{\baselinestretch}{1.5}

\begin{document}
\begin{flushright} CERN-TH.6892/93\\
\end{flushright}
\vspace*{10mm}
\begin{center}
{\bf Undesirable effects of covariance matrix techniques for error
analysis}\\[10mm]
David Seibert$^*$\\[5mm] Theory Division, CERN, CH-1211 Geneva 23,
Switzerland\\[10mm]
{\bf Abstract}\\
\end{center}
\hspace*{12pt}A version of the covariance matrix technique for treating
correlated data points can amplify systematic errors that are present
in the data.  In addition, this technique systematically underestimates
errors.  Both the amplification and the underestimate of the error are
uncontrolled, and arbitrarily inaccurate results can be obtained.  I
discuss alternative procedures for handling correlated data.\\
\vfill
\begin{center}
{\em Submitted to Phys.\ Rev.\ D}
\end{center}
\vfill
CERN-TH.6892/93\\
May 1993\\
\vspace*{10mm}
\footnoterule
\vspace*{3pt}
$^*$On leave until October 12, 1993 from:
Physics Department, Kent State University, Kent, OH 44242 USA.
Internet: seibert@surya11.cern.ch.\\

\newpage\setcounter{page}{1} \pagestyle{plain}
     \setlength{\parindent}{12pt}

%\section{Introduction}

Recently there has been some interest in the analysis of correlated
data, and people seeking more sophisticated analysis techniques have
often used the covariance matrix technique [1].  DeGrand [2], DeTar
and Kogut [3] and Gottlieb {\it et al}. [4] use the technique to
analyze lattice gauge theory results, while Abreu {\it et al}.\ [5]
and Wosiek [6] use it to analyze scaled factorial moment data.  I
show here that this analysis technique can produce systematic errors
that are not obtained with simpler, more robust techniques.

The covariance matrix technique is simple in principle -- transform
the data to an uncorrelated basis, fit the data in this basis, then
transform back to the laboratory frame.  However, some of the results
obtained by this procedure are very odd.  In particular, Gottlieb
{\it et al}. [4], Toussaint [7] and Wosiek [6] find that this procedure
can produce best-fit lines that fall below all data points, and even
below all error bars!

In this paper, I first discuss the proposed treatment of correlated
data, and show that in a {\em gedanken} experiment without systematic
errors this treatment produces exactly the desired results, but
underestimates the error.  In a very similar {\em gedanken} experiment
with arbitrarily small systematic errors, this procedure amplifies the
errors in the data; therefore, this treatment of data is not robust.  I
use these simple {\em gedanken} experiments instead of the scaled
factorial moment data for purposes of presentation, as the effect
observed in the two data sets is the same.  I then discuss
alternative procedures for handling correlated data.

%\section{Error estimation}

The experimental procedure is very simple.  Consider $N$ trial
measurements of $I$ data points, $y_i$.  Calculate the covariance matrix
from these data,
\begin{equation}
C_{ij} = \frac {1} {N-1} \left[ \sum_{n=1}^N (y_{i,n} -y_i) (y_{j,n} -y_j)
\right]~, \label{eCij}
\end{equation}
where
\begin{equation}
y_i ~=~ \frac {1} {N} \sum_{n=1}^N y_{i,n}~.
\end{equation}
Fit the data with the curve $f_i ( \alpha )$, where $\alpha$ is the set
of free parameters, by minimizing
\begin{equation}
\chi^2 = \sum_{i=1}^I \sum_{j=1}^I \, (y_i -f_i) \,
(C^{-1})_{ij} (y_j -f_j)~. \label{echi2}
\end{equation}

I illustrate this procedure with a {\em gedanken} experiment to measure
the mean voltage of a generator that produces random voltages $v$ with
probability distribution $p(v)$.  The generator charges a capacitor, and
I then measure the voltage on the capacitor twice, calling the two
measurements $y_1$ and $y_2$.  Each measurement has some (uncorrelated)
``measurement noise'' in addition to the fluctuations due to the random
voltage generator; I assume for now that this noise is different for the
two measurements.

After $N$ trials, the experimentally determined covariance matrix is
\begin{equation}
C = \frac {1} {N-1} \left( \begin{array} {cc}
                           \sigma^2 + e_1^2   & \sigma^2 \\
                           \sigma^2           & \sigma^2 + e_2^2
                \end{array} \right)~, \label{eCija}
\end{equation}
where
\begin{equation}
\sigma^2 = \int dv \, p(v) \, v^2 - \left[ \int dv \, p(v) \, v
\right]^2 \label{esp}
\end{equation}
is the contribution to the covariance matrix from the distribution of
random voltages, and $e_{1(2)}^2$ is the contribution of noise from
the set of first (second) measurements.  Fitting to the function $f_1 =
f_2 = V$, $\chi^2$ is minimized when
\begin{equation}
V = \frac {e_1^{-2} \langle y_1 \rangle + e_2^{-2} \langle
y_2 \rangle} {e_1^{-2} + e_2^{-2}}~, \label{eV1}
\end{equation}
where $\langle y_{1(2)} \rangle$ is the average value for the first
(second) measurement.  It is clear from Eq.~(\ref{eV1}) that $V$ is the
average of $\langle y_1 \rangle$ and $\langle y_2 \rangle$, properly
weighted for measurement error, and so the analysis procedure is very
successful at fitting the curve for this {\em gedanken} experiment.

%\section{Systematic effects}

The determination of the experimental error is not so simple.  The
minimum value of $\chi^2$ is
\begin{equation}
\chi^2_{min} ~=~ \frac {(N-1) \left( \langle y_1 \rangle - \langle y_2
\rangle \right)^2} {e_1^2 + e_2^2}~,
\end{equation}
and I could say that the experimental errors for the two points are
$e_1 / \sqrt{N-1}$ and $e_2 / \sqrt{N-1}$, respectively.  This gives
an error for $V$ of $(1/2)\sqrt{(e_1^2+e_2^2)/(N-1)}$.
However, this error reflects only the accuracy with which the voltages
were measured, ignoring fluctuations of the voltage itself.

If I want to determine the average voltage of the fluctuating source,
there is also error due to the fact that the measured voltages are
only a sample of the full voltage distribution.  Including this error,
the experimental error for $V$ is
$(1/2)\sqrt{(4\sigma^2+e_1^2+e_2^2)/(N-1)}$.  Thus, the experimental
error for $V$ is larger than the measurement error, and also
larger than the ``naive'' estimate from the fluctuations of the
individual points, $(1/2)\sqrt{(2\sigma^2+e_1^2+e_2^2)/(N-1)}$.
Unfortunately, it is not possible to calculate exactly  the different
errors from the data, so they must be estimated instead.

The simplest case is when the measurement error is small.  Then, I can
fit with high accuracy to the function $y_1=y_2=V$.  Fitting
both data points separately, the errors in $V$ are
$\sqrt{(\sigma^2+e_1^2)/(N-1)}$ and $\sqrt{(\sigma^2+e_2^2)/(N-1)}$;
the larger of these is greater than the true experimental error in
$V$.  The difference between this estimate and the true
error is small as long as all measurement errors are small compared to
$\sigma$.  Thus, the error can be constrained experimentally, as it is
bounded from above by measurement-to-measurement fluctuations of single
points.

Now, I modify the {\em gedanken} experiment slightly, by assuming that
the capacitor discharges somewhat between the two measurements.  I
therefore assume that the first measurement is unchanged, but that the
voltage is reduced by a factor $\gamma$ at the time of the second
measurement.  I could alternatively assume that the scales of the
voltmeters are slightly different, but I wish to have all systematic
effects occur before measurement rather than during it.

In this case, the experimentally determined covariance matrix is
\begin{equation}
C = \frac {1} {N-1} \left( \begin{array} {cc}
                  \sigma^2 + e_1^2   & \gamma \sigma^2 \\
                  \gamma \sigma^2    & \gamma^2 \sigma^2 + e_2^2
                \end{array} \right)~. \label{eCijb}
\end{equation}
After an infinite number of measurements $\langle y_1 \rangle = \overline{v}$
and $\langle y_2 \rangle = \gamma \overline{v}$, where $\overline{v}$ is the
mean value of the random voltage.  Fitting again to the function $f_1 = f_2 =
V$, $\chi^2$ is minimized when
\begin{equation}
V = \frac {\gamma e_1^2 + e_2^2} {e_1^2 + e_2^2 +
(\gamma-1)^2 \sigma^2} \, \overline{v}~. \label{eVb}
\end{equation}

The procedure gives {\em systematic} errors for this second experiment,
as $V$ is always less than $\overline{v}$.  This is not totally
unexpected, because the discharge of the capacitor between the
measurements gives a systematically lower value of $y_2$.  If
$\gamma=1$, as in the first experiment, all systematic errors vanish.
For any other value of $\gamma$, however, $V$ can have any value
between zero and $\overline{v}$.  By contrast, a naive least-squares fit
always yields a value for $V$ between $\gamma \overline{v}$ and
$\overline{v}$.  Thus, the covariance matrix technique can produce large
systematic errors from arbitrarily small intrinsic systematic errors.

Now I modify the {\it gedanken} experiment again,  by measuring three
voltages.  This increases the similarity to the scaled factorial moment
data [5, 6], which are expected to be fit by a straight line.  I measure
the voltages at $x=-1$, 0 and 1.  For each trial, the voltages are
$\gamma^{-1}v$, $v$ and $\gamma v$ respectively, and measurement error is
$e$ at all three points.  Thus, for $\gamma \neq 1$, the points do not
fall on a straight line, leading to systematic error.

The covariance matrix after $N$ trials is
\begin{equation}
C = \frac {1} {N-1} \left( \begin{array} {ccc}
    \gamma^{-2} \sigma^2 + e^2 & \gamma^{-1} \sigma^2 & \sigma^2 \\
    \gamma^{-1} \sigma^2 & \sigma^2 + e^2 & \gamma \sigma^2 \\
    \sigma^2 & \gamma \sigma^2 & \gamma^2 \sigma^2 + e^2
                \end{array} \right)~. \label{eCijc}
\end{equation}
Fitting to a straight line $y=v_0+\beta x$, I obtain (again after an
infinite number of trials)
\begin{eqnarray}
v_0 = \frac {2 (\gamma+1+\gamma^{-1}) (3\sigma^2+e^2) e^2 \overline{v}}
{\left[ (\gamma-2+\gamma^{-1})^2\sigma^2+6e^2 \right]
\left[ (\gamma^2+1+\gamma^{-2})\sigma^2+e^2 \right]}, \label{ev0} \\
\beta = \frac {(\gamma-\gamma^{-1})
\left[ (\gamma^2-4\gamma-3-4\gamma^{-1}+\gamma^{-2})\sigma^2+3e^2 \right]
e^2 \overline{v}} {\left[ (\gamma-2+\gamma^{-1})^2\sigma^2+6e^2 \right]
\left[ (\gamma^2+1+\gamma^{-2})\sigma^2+e^2 \right]}. \label{ebeta}
\end{eqnarray}
{}From eqs.~(\ref{ev0}) and (\ref{ebeta}), it is clear that $v_0, \beta
\rightarrow 0$ as $e/\sigma \rightarrow 0$, independent of the value
of $\gamma$.  Thus, the systematic errors on these two quantities are
again uncontrolled.  In fact, for $\gamma \simeq 1$, even the sign
of $\beta$ can be wrong when $\sigma^2 > e^2/3$.

%\section{Conclusions}

The amplification of the systematic errors disappears for fitting
functions that go through all of the data points.  This is trivial for
the {\em gedanken} experiments, as there are only two or three data
points and thus a two- or three-parameter fit will remove the systematic
effects.  However, it is clearly impossible in general, as all useful
fitting functions have fewer degrees of freedom than the number of data
points.

There are three possible alternative procedures for experimenters.  The
first method is to ignore the correlations and to perform a ``naive''
least-squares fit, assuming that all of the errors are equal.  The
second method is to estimate the covariance matrix theoretically, {\em
assuming that the data points are given by the curve $f_i (\alpha)$},
and then to minimize $\chi^2$ as in the above procedure.  This is
alternative is robust, but may involve a lot of computation.

The third method is to ignore the correlations as in the first method,
but to attempt to estimate the experimental errors for each of the data
points.  This will result in a somewhat better fit than the first
method, but obviously requires more work for the experimenter.  It is
usually less arduous than the second method, however, and will probably
produce equally satisfactory results.

In all cases, the experimental error in any quantity should not be
determined from $\chi^2$ used in the best-fit, either using the ``naive''
$\chi^2$ or the more sophisticated version constructed from the
covariance matrix.  The best approach is probably to give each parameter
of the best-fit curve the largest variance that reproduces
measurement-to-measurement fluctuations of individual data.  This will
in general give an overestimate of the experimental error, but it may
not be possible to improve on this procedure.

\bigskip %\ulsect{Acknowledgement}

This material is based upon work supported by the North Atlantic Treaty
Organization under a Grant awarded in 1991.

%\vfill \eject

\ulsect{References}

\begin{list}{\arabic{enumi}.\hfill}{\setlength{\topsep}{0pt}
\setlength{\partopsep}{0pt} \setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt} \setlength{\leftmargin}{\labelwidth}
\setlength{\rightmargin}{0pt} \setlength{\listparindent}{0pt}
\setlength{\itemindent}{0pt} \setlength{\labelsep}{0pt}
\usecounter{enumi}}

\item W.T. Eadie {\it et al}., {\it Statistical Methods in Experimental
Physics} (North-Holland, Amsterdam, 1971), pp.\ 62--66.

\item T.A. DeGrand, Phys.\ Rev.\ D {\bf 36}, 176 (1987).

\item C. DeTar and J.B. Kogut, Phys.\ Rev.\ D {\bf 36}, 2828 (1987).

\item S. Gottlieb {\it et al}., Phys.\ Rev.\ D {\bf 38}, 2245 (1988).

\item Abreu, P. {\it et al}., Phys.\ Lett.\ B {\bf 247}, 137 (1990).

\item B. Wosiek, Acta Physica Polonica {\bf B21}, 1021 (1990).

\item D. Toussaint, in {\it From Actions to Answers -- Proceedings of
the 1989 Theoretical Advanced Study Institute in Elementary Particle
Physics} (World Scientific, Singapore, 1990; T. DeGrand and D.
Toussaint, eds.).

\end{list}

\vfill \eject

\end{document}

%\vfill \eject

\ulsect{Figure caption}

\begin{list}{\arabic{enumi}.\hfill}{\setlength{\topsep}{0pt}
\setlength{\partopsep}{0pt} \setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt} \setlength{\leftmargin}{\labelwidth}
\setlength{\rightmargin}{0pt} \setlength{\listparindent}{0pt}
\setlength{\itemindent}{0pt} \setlength{\labelsep}{0pt}
\usecounter{enumi}}

\item

\end{list}

