Robust Fitting of Spectra to Splines with Variable Knots

 

R. L. Coldwell

 

Department of Physics, University of Florida, Gainesville, FL 32611 and
Constellation Technology Corporation, 7887 Brian Dairy Road, Largo, FL 33777

 

 

Spectra consist of continuum features that vary over many channels, and peaks that vary over few channels. In a fit to the continuum the peaks appear as outliers. Robust methods in which the weights associated with the peaks are reduced allow the continuum to be fitted almost independently of the peaks. This requires a very smooth background. Cubic splines, which are continuous with continuous first and second derivatives, are a good choice for this task when the locations of the discontinuities in the third derivatives, the knots, are included in the parameters optimized in the fitting process. An extension to the Marquardt method allows a Newton-Raphson method to be used in this optimization.

 

INTRODUCTION*

* This work supported in part by the Department of Defense Nuclear Treaty Programs Office through the United States Army Space and Missile Defense Command agent for the Defense Advanced Research Projects Agency and performed at the Pinellas Science, Technology, and Research Center under Grant DASG-609-610-007.

This is an html version of a paper submitted for possible inclusion in the proceedings of the 15th International Conference
the Applications of Accelerators in Research and Industry, Denton Texas Nov 4-7, 1998
A version with the vu-graphs actually used is in Stanfit3.html 

Spectra consist of continuum and peaks. It seems intuitive that the continuum should be smooth and as close to the bottom of the peaks as possible. This paper shows that this can be accomplished by fitting the continuum as the exponential of a cubic spline with adjustable knots. Most of the ability to fit the continuum comes simply from the stiffness of the splines. The addition of outlier regression, however, enables this stiffness to fit the continuum with no reference to the peaks.

The methods described here are embodied in the code RobWin being developed at Constellation Technology with windows programming by George P. Lasche, and analysis programming by the author.

The first few sections give a sketch of the mathematics of the non-linear regression methods that allow the continuum to be fitted with respect to the knot locations.

Curve fitting is tested only when the error estimates on the data points are very small. The data shown here was counted for 128 hours with distant sources and for 16 hours with close sources. In the first, an example, which easily separates into continuum and peaks, a nuclear spectrum taken with an intrinsic germanium detector, is shown. Then in the second, which tests the method presented here more severely, an 152Eu source was placed relatively close to a thin HgI2 Constellation Technology spectral-grade detector.

 

NON-LINEAR MINIMIZATION

 

 

Define

(1)

In equation 1 represents the set of constants in the fit and e i is the error estimate for each data point fi. The object is to minimize c 2 with respect to . The such that the derivatives of c 2 with respect to its components are all zero, starting from a sufficiently close to has components given by equation 2.

(2)

To find a minimum with the Newton-Raphson method of equation 2, it is necessary that be less than the radius of convergence of a quadratic expansion of equation 1. When this is not true the quantity minimized should be +c 2, where is a multidimensional version of the Marquardt parameter. For sufficiently large values of , this added term dominates and the minimum is found at , while in the limit, the Newton Raphson method is recovered. Owing to the ease with which the problem of linear dependence of the constants is solved with an extremely small value of some work with this gives the false impression in references (1) and (2) that the actual value of is unimportant. The Robmin method used in this work solves iteratively for the decrease in c 2 that can be achieved in each step. It uses the difference between the vector of first derivatives and the array of second derivatives for the predicted changes in to find the ratio of the l K's to each other. Finally it uses a one dimensional Newton's method to find the final for which the predicted c 2 is as much lower than the present one as can be achieved. The method is working properly far from the minimum when every third prediction results in a lower c 2. In the event that convergence appears to be established, but is very slow, Aitkin's extrapolation (3) is used in the early steps. At the extrema in c 2, |l |® 0. The best description of the current Robmin is in reference (4), while descriptions of earlier versions are in (5) and (6). In practice, the above simply gives such that the derivatives of c 2 are zero. If the weighted least squares fit is carried out as given with data that varies by many orders of magnitude, there are numerous extrema in addition to the desired global minimum. The solution is to begin with large error estimates in equation 1, slowly decreasing them to the maximum of those given by the actual statistics and smaller and smaller percentage error estimates.

 

 

SPLINES

 

The spline is the function representing the data with the smoothest possible second derivatives. (7) (8). A spline can be defined as

(3) Where

(4)

With c1=1 and x 1 and all others 0, this is

FIGURE 1 A single knot spline and its derivatives

Consider a peak like bump which is 1 at x=0, extending from -1 to 1 and is zero elsewhere, which we want to fit to a spline. This b-spline is

(5)

FIGURE 2 The bspline of equation 5

Five knots from -1 to +1, spanning the range of twice the fwhm of the peak are required to produce this. It is interesting to attempt to produce this result with only three knots.

FIGURE 3. Attempt to produce a bspline with 3 knots

The practical value of this is that a cubic spline will be able to reproduce a peak only if it has approximately 5 knots within a few half widths of the peak. In addition, the natural tendency of the splines defined here is to become large for small values of x. This is also the natural tendency for spectral data. The only coding needed to take advantage of this reluctance of a spline to produce a peak is the requirement that new knots be introduced at a distance from the old knots.

A positive-definite continuum is best fitted to the exponential of this spline. The hard part of fitting splines to data is the fitting of the knots. Allowing these to vary, however, significantly decreases the number of constants required to fit the continuum and more importantly eliminates the usual curve fit problem of oscillations in the fitting function.

The final detail needed to keep the background below the peaks is the outlier regression, which increases the error estimates for data above the continuum and decreases it for those below. This is described in detail in reference (6) pp. 50-54.

 

 

EXAMPLE WITH NARROW PEAKS

 

 

The following is an almost ideal example. It is the longest count time spectra from a test devised by George P. Lasche, of Constellation Technology. Kenneth W. Neufeld and Jan A. Nobel took the spectrum at Constellation Technology's Pinellas Plant just north of St. Petersburg Florida. It is composed of 16,382 channels over a nominal energy range from 0 to 3 MeV, and was counted for 128 hours. It was taken with a set of 10 certified laboratory calibration radioactive sources, each of about 10 microcuries, mounted together at a distance 20.48 meters (67.2 feet) with a 100%-efficient high-purity germanium detector.

FIGURE 4. White line is a 60-constant approximation to the continuum, fitted as the exponential of a spline.

The results of simply fitting the spectra in terms of the exponential of a spline are shown in figure 4. It appears fair, but the residuals contain negative regions, and a close up of the low energy region indicates that the numerous peaks in that region have confused the spline, which fits too high. The full featured fit, whose low energy region is shown in figure 5 and a high energy region in figure 6, adds outlier suppression after half the requested number of knots have been inserted and is biased to put the first group of knots at low channel numbers. This is a fit to the spectrum with no peaks, only a continuum. After this fit, the difference between the background and the peaks is ready to be fitted as either peaks or as a nuclide response function.


FIGURE 5. Low energy region with outlier regression.

FIGURE 6. Continuum near the 2614 KeV 208Tl line.

Note that the residuals in this fit are almost entirely positive owing to the utilization in the outlier regression equations of the fact that the continuum is lower than the peaks.

 

 

EXAMPLE WITH WIDE PEAKS

 

 

The source was 462 ± 0.5% KBq of 152Eu held at 16.8 ± 0.5 cm from crystal. The detector, an HgI2 crystal about 1 in2 x 2mm, was set for 0 to 3000 ± 4 KeV. Amanda M. Gerrish of Constellation Technology made the settings and took the 16-hour spectrum in March 1998.

Somewhat fewer constants need to be used on the first fit to avoid fitting peaks as part of the continuum. The result of fitting with 48 constants is shown in figure 8.

FIGURE 7. 48 constants fitted from channel 63 to 8060

It is not immediately obvious that this fit to the background allows the underlying 152Eu to properly be fitted. That this is so is best illustrated by completing the fit.

 

The Nuclide Response Function

The code RobWin allows the user to select the isotopes to be fitted to a nuclide reference function, which is then generated from the Brookhaven Nuclear Data Library (9). The conversion from a set of lines and intensities to a response function requires an energy calibration, a full width at half-maximum fit as a function of energy, and detector efficiency as a function of energy. The code was set to generate these as fits to the data. In addition templates representing the response of the detector to single isolated peaks are needed. The code was started with Gaussian templates, then a standard template was constructed from the 1408 KeV line by fitting the difference between the data and the background, the results of which are shown in figure 8.

FIGURE 8. High-energy template as fitted

The 121 KeV peak was also made into a standard and a small set of unspecified peaks was allowed. The fitting proceeded by making small changes, and allowing the rest of the fit to adapt to them. At one point the continuum and the efficiency function conspired to make the continuum go low. The continuum fit was then restarted with the unspecified peaks present to help at low energies, but without the nuclides. Then the nuclide fit was re-introduced, but the continuum not allowed to vary. The final fit was at a local minimum and all constants were allowed to vary. The results near the origin and in an intermediate region are shown in figures 9 and 10.

FIGURE 9. Final fit to HgI2 spectrum at low energies.


FIGURE 10. Fit in the intermediate energy region

 

The efficiency of the detector is in the nuclide response function as an exponential of a linear term and a 2-knot spline. The final rather innocuous looking efficiency is shown in figure 11. If it had been known in the beginning, as would be the case in future fits to this detector, the fit would have been much easier.


FIGURE 11. Efficiency of HgI2 detector.

 

CONCLUSION

 

The final c 2 of 10,240 for 7998 channels indicates that the fit is good, but not perfect. The ability to fit the continuum as the exponential of a cubic spline, however, has been demonstrated.

 

 

REFERENCES

 

 

1. Bevington, P. R., Data Reduction and Error Analysis for the
Physical Sciences, New York, McGraw Hill (1969)

2. Press, W. H., Teukolsky, S. A., Vetterling, W. T. , Flannery,
B. P., Numerical Recipes : The Art of Scientific Computing,

Cambridge University Press, Boston 1986 p. 523

3. Hildebrand, F. B., Introduction to Numerical Analysis,
McGraw-Hill (1956) pp. 443-445.

4. Alexander, S. A. and Coldwell, R. L. "Atomic Calculations
Using Variational Monte Carlo", in Recent Advances in
Quantum Monte Carlo Methods, ed W. A. Lester Jr., World
Scientific (1977) pp. 56-58

5. Coldwell, R.L., in Radiative Properties of Hot-Dense Matter,
Singapore, World Scientific (1983), pp. 315-349

6. Coldwell, R. L. and Bamford, G. J., The Theory and
Operation of Spectral Analysis Using Robfit, AIP, New York
(1991) pp. 50-54.

7. Hollady, J. C.,"A Smoothest Curve Approximation", Math.
Tables and Aids to Computation, 11, 233-243 (1957)
8. Ahlberg, J. H., Nilson E. N., and Walsh, J.E., The Theory of
Splines and Their Applications. New York, Academic Press
(1967), p. 3.

9. PCNUDAT Nuclear Data master database used by permission of NNDC at B.N.L.