Splines

Method

The constant CNMPKR is passed to the code CNMFIT.FOR. It determines the ratio of the background width to the foreground width. This is used in the SUBROUTINE DATAC in the bliep.for part of the code. The default size for this ratio is 4.

Data is compressed and peaks removed in the bliep.for part of the code. This began with the work in ERR.DOC which showed that the correct statistical weighting for points partly in one region and partly in another is 1/err² in one region and 1-1/err² in the second region. The integral of a Gaussian is . It is fitted to a pade approximate in Gaussians\nlfit2\Welcome to nlfit2.htm.

This function is 0.5 at z = 0 and is such that . Data is fitted by minimizing

In frequency space the product of the two AiGausses is shown in doc1/A trial d.htm#To_the_frequency_domain to be.

This weighting greatly suppresses frequency components greater than w. This is reinforced by a second “robust” fitting to the same data with the e(t) for these points increased to account for the fact that they are far enough from the fit to be outliers. The inverse of ½ the second derivative matrix is the error matrix for this. Thus this procedure for points spaced more than w apart yields a set of statistically independent points and their standard deviations.

The second part fitting the background involves fitting a spline to these points. The best set of results and hints for how to use bkgfit is in doc1\Aug23-04.doc. A tricky part mentioned in the paper below involves increasing the error on accurate sets of points. The method used here is described in doc1\Alpha for bkg.doc.

References

The paper described in

http://www.phys.ufl.edu/~coldwell/Const/denton/Denton3.html [. “Robust Fitting of Spectra to Splines with Variable Knots”, R.L. Coldwell, CP475, Application of Accelerators in Research and Industry, edited by J.L. Duggan and I.L Morgan AIP Press, New York (1999)] and the more detailed Stanfit3.html link under hint at the best method for fitting the continuum.

07/18/2004 10:19a Gaussians\Gaussians.doc Detailed simulation of sets of Gaussians separated by size.

The separation of data into foreground and background is described in #datac. The code fits data at constant CNMPKR in width tapered by AiGauss. That is the weights with a Gaussian fwhm are tapered both sides of a nominal end point such that inside the region the weights are 1-AiGauss and outside they are AiGauss. This makes the errors statistically correct for points that enter into two intervals. The minimium distance between intervals is set such that no point appears in more than two intervals.

The code CNMFIT.FOR calls bliep.for – which contains datac. Five times as many data points are generated as there are terms in a fitting spline. A linear fit to the data points requires ~ 1 second for enough terms ~ 36 to well represent the function. A better fit for enough terms with variable knots ~20 requires on the order of a minute.

Spline fitting

08/01/2004 10:05 AM 124,416 doc1\Cubic Spline Interpolation.doc is Press’s presentation of the codes SPLINE.FOR and SPLINT.FOR. These use x, y, and y” to interpolate a set of data points. The y” is solved for by requiring continuity of the first derivatives.

08/11/2004 10:26 AM 152,064 doc1\ypp.doc is Bob’s Lagrange derivation of the Press result.

08/11/2004 10:36 AM 101,376 doc1\The two forms.doc Contains code that converts from the x,y,y” form to the form

The partials with respect to c_j, d_k, and x_k are easy to evaluate in this form making it appropriate for fitting the data generated by datac.

08/01/2004 10:05 AM 124,416 doc1\Cubic Spline Interpolation.htm Flirts with midpoint knots, picture showing effects of end point derivatives.

08/11/2004 10:07 AM 139,776 doc1\Fitting Sections.doc Considers the boundary conditions for fitting data in short sections.

08/11/2004 12:52 PM 62,976 doc1\Fitting Sections 2.doc A bit more detail. In fact the fitting method actually used was to vary a knot, the six above it and the six below it. The knots were moved toward lower channels 2 knots at a time. The idea is that the top two knots and the bottom 2 represent the boundary conditions, while the rest are local. In fact, the non-local effect of each knot on regions much lower is probably the primary reason for the method working. The fact that this limits the number of constants to 26 makes the matrix invertible which possibly is more important than any other consideration. This section will be important for LIGO data.

08/11/2004 10:26 AM 152,064 doc1/ypp.htm – details the Lagrange interpolation of y and y”, produces the same results as doc1/Cubic Spline Interpolation.htm based on Press

08/11/2004 10:36 AM 101,376 doc1/The two forms.htm Detailed derivation of the code SUBROUTINE YPPTOC(X,Y,YPP,C0,C1,C2,C3,D,N) for converting x,y,y” to sum on c and d terms. This code is found in CNMFIT.FOR

08/14/2004 05:03 PM 24,576 doc1/Alpha for bkg.htm The notion that the error should be increased introduced in the general reference above is implemented here by solving for alpha such that 2N=(f-f_A)**2/((alpha*f)**2+e²)

08/21/2004 06:00 PM 18,432 doc1/Adding knots.htm A penalty is added to keep knots from drifting outside the fitting region.

08/22/2004 01:26 PM 194,560 doc1/Aug21-04.doc Documents some early results and glitches. Also has an unkdis display with approximate knot locations.

08/23/2004 05:52 AM 491,520 doc1/August22-04.doc Results for large numbers of knots.

08/23/2004 07:03 AM 336,896 doc1/Aug23-04.doc

Bkgfit

The code is

bkgfit.for à bkgfit.wpj Driver and overhead

C#INCLUDE bliep.for 07/28/2004 05:19 PM 62,464 doc1/Datac.doc

SUBROUTINE DATAC(F,WX,NDAT,FWHMS,IX,FI,WXF,NN,CNMPKR)

doc1/Datac.doc – This is the workhorse code. The fact that AIGAU(0)=0.5 and that AIGAU(-x)+AIGAU(x)=1

Is used to fit the data with a taper equal to a multiple, CNMPKR, of the fwhm. Peaks cannot be fully fitted with this form. The “robust” part of the code increases the error for the positive parts of the fit not fitted. A second pass then produces a background estimate for the data fitted. The minimum separation distance of the peaks fitted is made such that no point is fitted more than twice. Those fitted twice have their weights reduced to account for this fitting.

Note that NN may be returned as less than it was called with. It is limited to the number of data points divided by CNMPKR * fwhm.

C#INCLUDE CNMFIT.for à called by bkgfit – calls datac and outputs data if requested.

C#INCLUDE NLFITSUB.for à On watching this, I changed the distance to Aitkin’s from 5 to 25. The res had to be watched to avoid overflow in converting from a large double precision value to a real.

C#INCLUDE rbrack.for à xml file overhead

C#INCLUDE ROBMIN à sminv used by datac in fitting data sections, used directly by nlfitsub

C#INCLUDE RNUMB à xml file overhead

C#INCLUDE WSYSTEM àwcom\openwsys.for

C#INCLUDE BINDEX à i.o. overhead

C#INCLUDE XMLCOPY à xml file overhead

C#INCLUDE PKPOLY à used by bkgfit to subtract peaks from data

C#INCLUDE SPLINE à from Press

C#INCLUDE SPLINT. à from Press

The display part of the code is described in ..\..\..\unkdis\Welcome.doc.

Also relevant originally these were in vmspline2

06/09/2004 12:25 PM 41,472 doc1/Transforms.htm

07/13/2004 05:21 PM 26,112 doc1/Convolution Prime.htm A derivation of the convolution theorem

07/18/2004 06:01 AM 897,024 doc1/Peak fitting.htm

07/19/2004 09:43 AM 1,261,568 doc1/A trial d.htm – very useful. Shows results of narrow peaks superimposed on a background d of wide peaks.

07/21/2004 11:07 AM 62,464 doc1/Curve fitting.doc – the original method proposed – differs somewhat from the above.

The AiGauss function is defined and fitted in

Gaussians\nlfit2\Welcome to nlfit2.htm

Gaussians\Gauss integral fit.doc Simple curve showing that the Pade approximate fits AiGauss

Gaussians\AiGauss\Welcome to AiGauss.doc