Parametric Gradient of the Chi-Squared Statistic

In Forward Fitting (or MEM or PIXONS) one often needs the gradient of the χ^{2} statistic with respect to the model parameters. The χ^{2} statistic for photon count rates observed (C_{obs}) and modeled (C_{mod}) is:

χ^{2} = \frac{1}{N} ∑_{j} \frac{[C_{mod}^{(j)} − C_{obs}^{(j)}]^{2}}{σ_{j}^{2}} ,  (1)

where σ_{j}^{2} = C_{mod} in the case of pure Poisson statistics. An alternative form of the χ^{2} statistic is for visibilities observed (V_{obs}) and modeled (V_{mod}) is:

χ^{2} = \frac{1}{2N} ∑_{j} \frac{|V_{mod}^{(j)} − V_{obs}^{(j)}|^{2}}{σ_{j}^{2}}  (2)

where the sum is over the index j of time bins, and the variance σ_{j}^{2} can either be independent of the model parameters (as in radio astronomy) or can be proportional to the mean count rate (\frac{1}{2}|V_{obs}|).

1 Visibilities

For the case of visibilities, the gradient of χ^{2} is:

\frac{∂χ^{2}}{∂a_{k}} = \frac{1}{2N} \frac{∂}{∂a_{k}} ∑_{j} \frac{|V_{mod}^{(j)} − V_{obs}^{(j)}|^{2}}{σ_{j}^{2}}  (3)

It might be thought that evaluation of this expression is only possible numerically, but it is straightforward in many cases to write an analytic expression for (3). This will speed up the search through parameter space, and eliminates the insidious round-off errors commonly found in numerical differencing schemes.

To illustrate this, we consider the particular case:

1. σ_{j}^{2} is not a function of the parameters a_{k},

2. V_{mod} derives from a sum of Gaussian sources.

Both of these conditions can be relaxed, and similar analytic expressions can be found where σ_{j}^{2} derives from Poisson statistics, and for more general basis functions than Gaussians. For the Orthogonal Gaussian model, the visibility is given by a sum over terms like:

V_{mod}^{(j)} = a_{0}e^{−a_{1}u_{j}^{2}−a_{2}v_{j}^{2}}e^{ia_{3}u_{j}+ia_{4}v_{j}},  (4)

for sources of width (\sqrt{a_{1}}, \sqrt{a_{2}})/π, location (a_{3}, a_{4})/2π, and flux a_{0}. For multiple Gaussians, the a_{k} parameters would take on an additional subscript n, and V_{mod}^{(j)} would be a sum over n, but for simplicity we ignore that subscript in the following. The derivatives of V_{mod}^{(j)} with respect to a_{k} are simple:

\frac{∂}{∂a_{k}} V_{mod}^{(j)} = f_{k} V_{mod}^{(j)}  (5)

where the vector f = [1/a_{0}, −u_{j}^{2}, −v_{j}^{2}, iu_{j}, iv_{j}]. Note that f has 3 real elements and 2 imaginary elements. This leads to simplifications in the derivatives below. Recalling our assumption that σ_{j}^{2} is independent of a_{k}, the gradient of the numerator in (3) is:

\frac{∂}{∂a_{k}} |V_{mod}^{(j)} − V_{obs}^{(j)}|^{2} = (V_{mod}^{(j)∗} − V_{obs}^{(j)∗})f_{k}V_{mod}^{(j)} + (V_{mod}^{(j)} − V_{obs}^{(j)})f_{k}^{∗}V_{mod}^{(j)∗} (6)

The above expression is real because the 3rd and 4th terms on the RHS are the complex conjugates of the 1st and 2nd terms. If f_{k} is real, the RHS becomes:

RHS = f_{k}[2|V_{mod}^{(j)}|^{2} − (V_{obs}^{(j)∗}V_{mod}^{(j)} + V_{obs}^{(j)}V_{mod}^{(j)∗})]  (7)

and if f_{k} is pure imaginary,

RHS = f_{k} [V_{obs}^{(j)}V_{mod}^{(j)∗} − V_{obs}^{(j)∗}V_{mod}^{(j)}]  (8)

So the equation for the gradient simplifies considerably. We can get a still simpler form by using the amplitude and phase of V_{obs} and by using the location parameter a_{3} u_{j} + a_{4} v_{j} as a phase:

V_{obs} = |V_{obs}| e^{iψ}  (9)
ϕ = a_{3} u_{j} + a_{4} v_{j}  (10)

Then we find the gradient to be:

\frac{∂χ^{2}}{∂a_{k}} = \frac{a_{0}}{N} ×  (11)
∑j (f_{k}/σ_{j}^{2})e^{−a_{1}u_{j}^{2}−a_{2}v_{j}^{2}}[a_{0} − |V_{obs}|cos(ψ_{j} − ϕ_{j})] k = 0, 1, 2 (12)
∑j (if_{k}/σ_{j}^{2})e^{−a_{1}u_{j}^{2}−a_{2}v_{j}^{2}}[−|V_{obs}|sin(ψ_{j} − ϕ_{j})] k = 3, 4  (13)

These can be combined in the form:

\frac{∂χ^{2}}{∂a_{k}} = \frac{a_{0}}{N} Re\{∑_{j}(f_{k}/σ_{j}^{2})e^{−a_{1}u_{j}^{2}−a_{2}v_{j}^{2}}[a_{0} − V_{obs}e^{−i(a_{3}u_{j}+a_{4}v_{j})}]\}  (14)

It is important to use these equations starting with a_{1} and a_{2} very small, because there is a spurious zero to the gradient of χ^{2} at a_{1} = ∞ or a_{2} = ∞.

2 Count Rates

The model count rate profile can be written as a function of the phase at map center Φ:

C_{mod} = a_{0}τ(1 + c_{1} cosΦ_{jm}) + higher harmonics  (15)

The phase term Φ_{jm} is the sum of the phase at map center Φ_{j} and the phase offset of a point source at Δx, Δy relative to map center:

Φ_{jm} = Φ_{j} + 2π (u_{j} Δx + v_{j} Δy)  (16)

The variables (u_{j}, v_{j}) are the Fourier coordinates on the u, v circle for a given subcollimator. In practice, during an iteration of the reconstruction, Φ_{j} need be computed only once, while the offset terms change during the iteration. The fundamental amplitude c_{1} for an orthogonal Gaussian source is given by:

c_{1} = e^{−u_{j}^{2}a_{1}−v_{j}^{2}a_{2}} · A_{m},  (17)

where A_{m} is the modulated amplitude (modamp) obtained from the calibrated event list for a point source at map center. For the case where the σ_{j} are independent of a_{k}, the derivatives of the χ^{2} function with respect to the a_{k} are:

\frac{∂χ^{2}}{∂a_{0}} = \frac{2}{N} ∑_{j} \frac{C_{mod}^{(j)} − C_{obs}^{(j)}}{σ_{j}^{2}} C_{mod}/a_{0}  (18)

\frac{∂χ^{2}}{∂a_{1}} = \frac{2}{N} ∑_{j} \frac{C_{mod}^{(j)} − C_{obs}^{(j)}}{σ_{j}^{2}} a_{0}τ(−u_{j}^{2})cosΦ_{jm}  (19)

\frac{∂χ^{2}}{∂a_{2}} = \frac{2}{N} ∑_{j} \frac{C_{mod}^{(j)} − C_{obs}^{(j)}}{σ_{j}^{2}} a_{0}τ(−v_{j}^{2})cosΦ_{jm}  (20)

\frac{∂χ^{2}}{∂a_{3}} = \frac{2}{N} ∑_{j} \frac{C_{mod}^{(j)} − C_{obs}^{(j)}}{σ_{j}^{2}} (−a_{0}τc_{1}u_{j})sinΦ_{jm}  (21)

\frac{∂χ^{2}}{∂a_{4}} = \frac{2}{N} ∑_{j} \frac{C_{mod}^{(j)} − C_{obs}^{(j)}}{σ_{j}^{2}} (−a_{0}τc_{1}v_{j})sinΦ_{jm}  (22)

(23)

For the pure Poisson statistics case, σ_{j}^{2} = C_{mod}, which is a function of the parameters a_{k}, so χ^{2} must be written as:

χ^{2} = \frac{1}{N} ∑_{j}\{C_{mod}^{(j)} − 2C_{obs}^{(j)} + [C_{obs}^{(j)}]^{2}/C_{mod}^{(j)}\};  (24)

and the χ^{2} gradient assumes the form:

\frac{∂χ^{2}}{∂a_{0}} = \frac{2}{N} ∑_{j}\{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}\}C_{mod}/a_{0}  (25)

\frac{∂χ^{2}}{∂a_{1}} = \frac{2}{N} ∑_{j}\{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}\}a_{0}τ(−u_{j}^{2})cosΦ_{jm}  (26)

\frac{∂χ^{2}}{∂a_{2}} = \frac{2}{N} ∑_{j}\{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}\}a_{0}τ(−v_{j}^{2})cosΦ_{jm}  (27)

\frac{∂χ^{2}}{∂a_{3}} = \frac{2}{N} ∑_{j}\{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}\}(−a_{0}τc_{1}u_{j})sinΦ_{jm}  (28)

\frac{∂χ^{2}}{∂a_{4}} = \frac{2}{N} ∑_{j}\{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}\}(−a_{0}τc_{1}v_{j})sinΦ_{jm}  (29)

There is a singularity in these equations when C_{mod}^{(j)}] = 0, and any (u, v) points for which this is true must be treated separately.

3 Cash Statistic

In the case of very low count rates, one must use what might be called the Cash statistic (Cash, 1979):

χ_{Cash}^{2} = \frac{2}{N} ∑_{j=1}^{N}[n_{j}ln(\frac{n_{j}}{e_{j}}) − (n_{j} − e_{j})],  (30)

where the n_{j} are integers, representing the observed counts per bin, and the e_{j} are arbitrary real numbers representing the model count rates. (In terms of previous variables, n_{j} = C_{obs} and e_{j} = C_{mod}.) In the limit of large n_{j} and e_{j}, χ_{Cash}^{2} becomes the standard χ^{2} statistic, and when e_{j} → n_{j}, χ_{Cash}^{2} → 0. The gradient of the Cash statistic with respect to the model e_{k} is:

\frac{∂χ_{Cash}^{2}}{∂e_{k}} = \frac{2}{N} (1 − n_{k}/e_{k})  (31)

Just as in the Poisson version above (equation 25 - 29), there is a singularity where e_{k} = 0, and such bins must be treated separately. Note that for low count rates, the derivatives in equations 25 - 29 must have the coefficients {1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]^{2}} replaced by 2{1 − [C_{obs}^{(j)}/C_{mod}^{(j)}]}, and these are equal when n_{j} is large.