%-----------------------------------------------------------------------
%; Copyright (C) 1996
%; Associated Universities, Inc. Washington DC, USA.
%;
%; This program is free software; you can redistribute it and/or
%; modify it under the terms of the GNU General Public License as
%; published by the Free Software Foundation; either version 2 of
%; the License, or (at your option) any later version.
%;
%; This program is distributed in the hope that it will be useful,
%; but WITHOUT ANY WARRANTY; without even the implied warranty of
%; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%; GNU General Public License for more details.
%;
%; You should have received a copy of the GNU General Public
%; License along with this program; if not, write to the Free
%; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
%; MA 02139, USA.
%;
%; Correspondence concerning AIPS should be addressed as follows:
%; Internet email: aipsmail@nrao.edu.
%; Postal address: AIPS Project Office
%; National Radio Astronomy Observatory
%; 520 Edgemont Road
%; Charlottesville, VA 22903-2475 USA
%-----------------------------------------------------------------------
\documentstyle{article}
\setlength{\topmargin}{-8mm}
\setlength{\textwidth}{160mm}
\setlength{\textheight}{220mm}
\setlength{\oddsidemargin}{0mm}
\setlength{\evensidemargin}{45mm}
\setlength{\parindent}{3em}
\setlength{\unitlength}{10mm}
\input{epsf.sty}
%\def\plotfiddle#1#2#3#4#5#6#7{ % \begin{picture}(10,8)
% \centering \leavevmode
% \vbox to#2{\rule {0pt}{#2}}
% \special{psfile=#1 voffset=#7 hoffset=#6 vscale=#5 hscale=#4 angle=#3}}
% % \end{picture}
% \plotfiddle {taper.ps}{60mm}{0}{1}{1}{0}{0}
\def\plotone#1{\begin{picture}(12,7)(-2,-0.8)
\centering \leavevmode
\epsfysize=6\unitlength \epsfbox{#1}
\end{picture}}
\def\plottwo#1#2{\centering \leavevmode
\epsfxsize=5.5\unitlength \epsfbox{#1} \hfil
\epsfxsize=5.5\unitlength \epsfbox{#2}}
\newcommand{\nl}{\newline}
\title{\bf ERRORS IN TWO DIMENSIONAL GAUSSIAN FITS}
\vspace{2 mm}
\author{ L.~Kogan
\vspace{2 mm}\\
\small National Radio Astronomy Observatory, Socorro, New Mexico,
USA\\}
%\date{September~ 5,~1993}
\vspace{2mm}
\begin{document}
\maketitle
\vspace{5mm}
\begin{abstract}
Two-dimensional elliptical Gaussian fits are used in astronomy for
accurate measurements of source parameters such as central position,
peak flux density and angular size. Condon \cite{jjcon} has provided
an analysis of error propagation for two-dimensional elliptical
Gaussian fits in coordinates along the major and minor axes. This memo
concerns the implementation of error estimates within the fitting task
JMFIT. The problem of error estimates can be considered more generally
in the RA-DEC coordinate system for an ellipse of arbitrary position
angle. This analysis is presented and a difference at the different
ratio of the beam and the fitted gausian size is discussed. Three
sets of formulae for the errors have been obtained: 1) the case of the
point source (the beam size $\simeq$ the gaussian size); 2) the case
of a resolved source (the beam size $\ll$ the gaussian size); and 3)
The intermediate case.
\end{abstract}
\section{Introduction}
The fitting of two-dimensional elliptical Gaussians is common in
astronomy and the appropriate error analysis for the six fitted
parameters is very important. Condon \cite{jjcon} has provided an
analysis of error propagation for two-dimensional elliptical Gaussian
fits in coordinates along the major and minor axes. Such a special
coordinate system does not allow the determination of all six
parameter errors for an ellipse of arbitrary position angle. Using the
linear transformation between two cartesian coordinate systems, the
following expression can be used however:
\begin{eqnarray}
D\,X_1 = D\,X \cdot \sin^2 \varphi \;+\; D\,Y \cdot \cos^2 \varphi
\nonumber \\
D\,Y_1 = D\,X \cdot \cos^2 \varphi \;+\; D\,Y \cdot \sin^2 \varphi
\end{eqnarray}
\begin{tabbing}
where~ \= \(D\,X_1,\; D\,Y_1\) are variances of the central position at the RA-DEC coordinate system;\\
\> \(D\,X,\; D\,Y\) are variances of the central position at the
coordinates along the major \\
\> and minor axes;\\
\> \( \varphi\) is position angle of the major axis
\end{tabbing}
Such a formula is implemented in the AIPS task JMFIT. It merits
further investigation however, because, a dependence on the variance
of position angle may be expected. This problem can be considered more
generally by performing an analysis of error propagation for
two-dimensional Gaussian fits in an RA-DEC coordinate system.\\
Many formulae in this memo were taken from the Condon \cite{jjcon} paper.
\section{Error~ Propagation~ in an ~RA-DEC ~Coordinate ~System:
Independent Samples.}
Condon \cite{jjcon} obtained the following matrix using coordinates
along the major and minor axes:
\begin{equation}
M = \frac{\pi A}{2}
\left(
\begin{array}{cccccc}
\frac{A\sigma_y}{\sigma_x} & 0 & 0 & 0 & 0 & 0 \\
0 & \frac{A\sigma_x}{\sigma_y} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{A\sigma_x\sigma_y}{2} & 0 & 0 & 0 \\
0 & 0 & 0 & \frac{3A\sigma_y}{2\sigma_x} & \frac{A}{2} & \sigma_y \\
0 & 0 & 0 & \frac{A}{2} & \frac{3A\sigma_x}{2\sigma_y} & \sigma_x \\
0 & 0 & 0 & \sigma_y & \sigma_x & \frac{2\sigma_x \sigma_y}{A} \\
\end{array}
\right)
\begin{array}{c}
X \\ Y \\ \beta \\ W_x \\ W_y \\ A
\end{array}
\label{eq:m}
\end{equation}
\begin{tabbing}
where~
\=the right-hand column denotes the sequence of the fitted parameters;\\
\> X,Y are coordinates of the ellipse center;\\
\> $\beta$ is a parameter determining the position angle of the major axis;\\
\> $W_x, W_y$ are the width of the gaussian at the major and minor axis direction;\\
\> A is the maximum value of the gaussian;
\end{tabbing}
The widths of the Gaussian and location of its center are measured
throughout this memo in pixels. The parameter sequence has been
rearranged for the matrix to have diagonal type to simplify evaluation
of the inverse matrix. In an RA-DEC coordinate system, matrix
(\ref{eq:m}) takes the form
\begin{equation}
M = \frac{\pi A}{2}
\left(
\begin{array}{cccccc}
CS & S2F & 0 & 0 & 0 & 0 \\
S2F & SC & 0 & 0 & 0 & 0 \\
0 & 0 & DF^{-1} & 0 & 0 & 0 \\
0 & 0 & 0 & \frac{3A}{2\,\alpha} \frac{\theta_m}{\theta_M}
& \frac{A}{2\alpha} & \frac{\theta_m}{\alpha}\\
0 & 0 & 0 & \frac{A}{2\alpha} & \frac{3A}{2\,\alpha} \frac{\theta_M}{\theta_m}
& \frac{\theta_M}{\alpha}\\
0 & 0 & 0 & \frac{\theta_m}{\alpha} & \frac{\theta_M}{\alpha}
& \frac{2 \theta_M \theta_m}{\alpha A}
\end{array}
\right)
\begin{array}{c}
X \\ Y \\ \varphi \\ \theta_M\\ \theta_m\\ A
\end{array}
\label{eq:m1}
\end{equation}
\begin{tabbing}
where~
\=the right column denotes the sequence of the fitted parameters;\\
\> X,Y are coordinates of the ellipse center;\\
\> $\varphi$ is the position angle of the major axis;\\
\> $\theta_M, \theta_m$ are the length of the major and minor ellipse axes at the half amplitude point;\\
\> A is the maximum value of the Gaussian; \\
\> $\alpha = 8\,\ln 2$\\
\> $CS = A\, \left(\frac{\theta_m}{\theta_M} \sin^2 \varphi +
\frac{\theta_M}{\theta_m} \cos^2 \varphi \right)$\\
\> $SC = A\, \left(\frac{\theta_m}{\theta_M} \cos^2 \varphi +
\frac{\theta_M}{\theta_m} \sin^2 \varphi \right)$\\
\> $S2F = \frac{A}{2} \frac{\theta_M^2 - \theta_m^2}
{\theta_M \theta_m} \sin 2\varphi$\\
\> $DF^{-1} = \frac{A}{2 \alpha} \frac{\left(\theta_M^2 - \theta_m^2\right)^2}
{\theta_M \theta_m}$\\
\end{tabbing}
Matrices (\ref{eq:m}) and (\ref{eq:m1}) differ by two additional
nonzero terms which demonstrate the correlation of the center
coordinate errors, if the major and minor axes do not coincide with
the coordinate axes. The inverse matrix $M^{-1}$ (the error matrix)
can be evaluated in this case as a combination of three inverse
matrices: a 2x2 top-left matrix, a 1x1 middle matrix and a 3x3
bottom-right matrix. The resulting error matrix takes the form:
\pagebreak
\begin{equation}
M^{-1} = \frac{2}{\pi A}
\left(
\begin{array}{cccccc}
\frac{SC}{D} & -\frac{S2F}{D} & 0 & 0 & 0 & 0 \\
-\frac{S2F}{D} & \frac{CS}{D} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1}{DF^{-1}} & 0 & 0 & 0 \\
0 & 0 & 0 & \frac{\alpha}{A}\frac{\theta_M}{\theta_m}
& 0 & -\frac{\alpha}{2\theta_m}\\
0 & 0 & 0 & 0 & \frac{\alpha}{A}\frac{\theta_m}{\theta_M}
& -\frac{\alpha}{\theta_M}\\
0 & 0 & 0 & -\frac{\alpha}{2\theta_m}& -\frac{\alpha}{2\theta_M}
& \frac{A \alpha}{\theta_M \theta_m}
\label{eq:m-1}
\end{array}
\right)
\begin{array}{c}
X \\ Y \\ \varphi \\ \theta_M\\ \theta_m\\ A
\end{array}
\end{equation}
\begin{tabbing}
where~\= \\
\> $D = CS \cdot SC \,-\, S2F^2$ is the determinant of the top-left 2x2 matrix.
\end{tabbing}
In spite of the fact that the initial expression for $D$ includes a
rather complicated dependence on position angle $\varphi$, after
simplification this dependence disappears and the determinant is given
by the very simple formula: $D = A^2$ !!!
If samples of amplitude measurements are independent and have the same
rms $\sigma$, the variance of any differentiable function ${\bf F}$ of
the six fitted Gaussian parameters is given by the following function
of $\sigma$ and elements of matrix $M^{-1}$:
\begin{equation}
VAR(F) = \sigma ^2 \sum_{i=1}^{6}
\sum_{j=1}^{6} M_{ij}^{-1} \frac {\partial F}
{\partial p_i} \frac {\partial F}{\partial p_j}
\label{eq:sum}
\end{equation}
In particular, the variance of the fitted parameters is determined
exclusively by the diagonal elements of the matrix :
\begin{equation}
VAR(p_i) = \sigma ^2 M_{ii}^{-1}
\end{equation}
\begin{eqnarray}
\frac{2}{\rho^2} & = &
\frac{2\alpha}{\pi\theta_M \theta_m}\frac{\sigma^2}{A^2} = \frac{VAR (A)}{A^2} =
\frac{VAR (\theta_m)}{\theta_m^2} =\frac{VAR (\theta_M)}{\theta_M^2} =
\frac{VAR (X)}{ \theta_M^2 \sin^2 \varphi + \theta_m^2 \cos^2 \varphi}\alpha =
\nonumber \\
& = &
\frac{VAR (Y)}{ \theta_m^2 \sin^2 \varphi + \theta_M^2 \cos^2 \varphi}\alpha =
\frac{VAR (\varphi)}{2} \left( \frac{\theta_M^2 - \theta_m^2}{\theta_m \theta_M}
\right)^2
\label{eq:ro}
\end{eqnarray}
The number $\rho$ specified the overall signal-to-noise ratio of the
Gaussian was implemented by Condon \cite{jjcon} and is determined by
the following equation in the case of independent samples:
\begin{equation}
\rho^2 = \frac{\pi}{\alpha} \frac{A^2}{\sigma^2} \theta_M \theta_m
\end{equation}
The integrated amplitude is defined by:
\begin{equation}
I = \frac{2\pi}{\alpha} \theta_M \theta_m A
\end{equation}
Its variance inferred from equation~ (\ref{eq:sum}) equals:
\begin{eqnarray}
VAR(I) &=& \sigma^2 [ M_{66}^{-1}\left(\frac{\partial I}{\partial A}\right)^2
+ M_{55}^{-1}\left(\frac{\partial I}{\partial \theta_m}\right)^2
+ M_{44}^{-1}\left(\frac{\partial I}{\partial \theta_M}\right)^2 +\nonumber
\\
&+& 2 M_{56}^{-1}\left(\frac{\partial I}{\partial A}\right)
\left(\frac{\partial I}{\partial \theta_m}\right) +
2 M_{46}^{-1}\left(\frac{\partial I}{\partial A}\right)
\left(\frac{\partial I}{\partial \theta_M}\right)
]
\label{eq:int}
\end{eqnarray}
The fourth and fifth items of the sum in equation~ (\ref{eq:int}) are
negative because errors of maximum value and width of the ellipse are
anti-correlated (as seen by examining the corresponding coefficients
of the matrix equation~(\ref{eq:m-1}). Forming the sum in equation~
(\ref{eq:int}) shows that the absolute value of the second and third
items are equal to the absolute value of the fourth and fifth items
respectively. So they cancel and the variance of the integrated
amplitude is given by the first term only:
\begin{equation}
\frac{VAR(I)}{I^2} = \frac{VAR(A)}{A^2}
\label{eq:inti}
\end{equation}
This result is valid only for independent samples. The general case
is more complicated and is partially analyzed in the following
section.
\section{Correlated Noise.}
Generally speaking, there is no image having independent pixel noise,
because any image is the result of convolution of an actual image with
the synthesized beam. Mathematically it is extremely difficult to
provide a full analysis of error propagation for dependent pixel noise
if the beam and fitted component have arbitrary size and orientation.
Condon \cite{jjcon} has analyzed the case of circular beam. Let's
consider the two limiting cases: 1) beam size is much bigger original
source size i.e. image coincides with beam; 2) beam size is much
smaller than original source size. For these two limiting cases
Condon's analysis can be expanded to a general form of the beam by
substituting the square of the beam diameter $\theta_N^2$ by the
product of the major and minor axes of the beam $\theta_{bM
}\theta_{bm}$. Now let's look at the two limiting cases separately.\\
\\
\underline{\em Expanded source. Beam size is much smaller than
original source size.}\\
Following Condon \cite{jjcon} a new expression for $\rho$ has to be
used instead of equation (\ref {eq:ro}):
\begin{equation}
\frac{2}{\rho^2} = 8 \frac{\sigma^2}{A^2} \frac{\theta_{bM}\theta_{bm}}
{\theta_{M}\theta_{m}}
\label{eq:ro1}
\end{equation}
Equation (\ref {eq:ro}) can be used to determined errors of the six
fitted parameters but using $\rho$ taken from equation (\ref{eq:ro1}).
The error in the integrated amplitude (\ref{eq:inti}) is transformed
to:
\begin{equation}
\frac{VAR(I)}{I^2} = \frac{VAR(A)}{A^2} + \frac{\theta_{bM}\theta_{bm}}
{\theta_{M}\theta_{m}} \left[\frac{VAR(\theta_M)}{\theta_M^2} + \frac{VAR(\theta_m)}{\theta_m^2} \right]
\label{eq:int2}
\end{equation}
\\
\underline{\em Point source. Image coincides with beam.}\\
Parameter $\rho$ is given in this case by the following expression:
\begin{equation}
\frac{2}{\rho^2} = \frac{\sigma^2}{A^2}
\label{eq:ro2}
\end{equation}
The six fitted parameters can be obtained from (\ref {eq:ro}) using
equation (\ref{eq:ro2}) for $\rho$. In particular $VAR(A) = \sigma^2$.
The error in integrated amplitude now takes the form:
\begin{equation}
\frac{VAR(I)}{I^2} = \frac{VAR(A)}{A^2} + \frac{VAR(\theta_M)}{\theta_M^2} + \frac{VAR(\theta_m)}{\theta_m^2}
\end{equation}
\\
\underline{\em Intermediate case. The source is partially resolved}\\
This case was analyzed by Condon \cite{jjcon} for a circular beam. For
a general beam this problem is very difficult. A reasonable approach
to the problem can be interpolation of expression for
$\frac{2}{\rho^2}$ using equation (\ref {eq:ro}) and (\ref{eq:int2}):
\begin{equation}
\frac{2}{\rho^2}=\frac{VAR(A)}{A^2} = \left\{
\begin{array}{ll}
8 \frac{\sigma^2}{A^2} \frac{\theta_{bM}\theta_{bm}}{\theta_{M}\theta_{m}}
& \mbox{if~~ $0.0 \leq \frac{\theta_{bM}\theta_{bm}}{\theta_{M}\theta_{m}} < 0.1$}\\
\frac{\sigma^2}{A^2}\left[0.8 \;+\; \frac{1}{4}\left(\frac{\theta_{bM}\theta_{bm}}{\theta_{M}\theta_{m}} - 0.1\right)\right]
& \mbox{if~~ $0.1 \leq \frac{\theta_{bM}\theta_{bm}}{\theta_{M}\theta_{m}} < 0.9$}\\
\frac{\sigma^2}{A^2} & \mbox{if~~ $0.9\leq \frac{\theta_{bM}\theta_{bm}}{\theta_{M}\theta_{m}} < 1.0$}
\end{array}
\right.
\label{eq:three}
\end{equation}
\section{Conclusion. What Has Been Implemented Within JMFIT?}
Several modifications in the errors calculation within JMFIT were
implemented as a result of this memo:\\
1. The old version of JMFIT had one set of formulae independent of
the ratio of the beam size to the image size. The new version
(15OCT96) analyses this ratio and uses three sets of the formulae as
given in equation (\ref{eq:three}) above.\\
2. The old version of JMFIT provided two sets of errors based on the
rms obtained from the residual map and the rms from the data itself,
which caused confusion. Only one set of errors based on the rms
obtained from the data itself has been retained. A new adverb has been
added to allow the user to force a pre-determined rms of samples to be
used in the error calculations.\\
3. The variance of the maximum value $A$ can be halved by fixing the
values of width and position angle of the gaussian ellipse. It follows
from the anti-correlation of the errors of the ellipse size and the
maximum value \cite{jjcon}. This reduction of the error of the maximum
value has been implemented in JMFIT when adverb GWIDTH is used to fix
the major and minor axes and position angle leaving only three free
parameters: maximum value and $X,Y$ coordinates of the ellipse center.\\
4. A typographical error in the calculation of the position angle
error has been corrected.
\begin{thebibliography}{99}
\bibitem{jjcon} J.J. Condon, Preprint of National Radio Astronomy
Observatory May, 1996
\end{thebibliography}
\end{document}