Thursday, January 24, 2013

XSPEC v12.8 - modifications to fit command

The new XSPEC release includes some improvements to the fit command. The most visible is an additional column in the output during the fit. |beta|/N is the norm of the vector of derivatives of the statistic with respect to the parameters divided by the number of parameters. At the best fit this should be zero so provides another measure of how well the fit is converging. |beta|/N can also be used as the criterion to stop the fit instead of the statistic delta although this is still considered experimental. This feature is activated by an additional argument on the fit command. For instance, fit 50 1.0e-3 20.0 will run the fit until either 50 iterations are completed, the statistic difference between successive iterations is less than 1.0e-3 or |beta|/N is less than 20.0. Setting either the statistic difference or |beta|/N stopping criteria to a negative number means they will not be used. These stopping criteria can also be set using the method command.

Another internal change is the option to use delayed gratification. This is turned on by giving the argument delay to the fit or method leven commands and turned off by giving the argument nodelay. Delayed gratification is a modification to the choice of the Levenberg-Marquardt damping factor and has been shown to speed convergence in many cases. The default is nodelay.

One other change included in the release was for the first iteration to only modify normalization parameters. We have subsequently determined that this introduces instabilities in some cases so have removed it in patch 12.8.0c. This patch also fixes an infinite loop case which has shown up in a few cases when running the error command.

Error in simftest.tcl

Dacheng Lin has pointed out that the Tcl script simftest.tcl which can be run from the XSPEC command line using the command simftest does not work correctly. The call to the lrt.tcl script within simftest switched the order of the names of the model command files. This has been fixed in the 12.8.0d patch available from the usual site. This patch also adds options to simftest and lrt to specify a filename into which the simulation results are written. Note that documentation on simftest and lrt can be obtained within XSPEC by typing the command name.

Thursday, January 17, 2013

XSPEC v12.8 - new statistics

XSPEC uses statistics for two purposes: 1) determining the best-fit parameter and its uncertainties; and 2) checking whether the model is a good fit to the data. In the past these have both been performed using the same statistic. There are two problems with this. Firstly, there are many goodness-of-fit statistical tests which do not use a statistic appropriate for parameter estimation. Secondly, using the same statistic for both parameter estimation and goodness of fit affects the properties of the statistic when used for the latter purpose.

XSPEC v12.8 now allows separate statistics to be defined for parameter estimation and goodness of fit. There are no changes to the statistic command for parameter estimation however to use a new goodness-of-fit statistic the command becomes statistic test whatever. At the end of a fit both statistic values are now written out. The goodness command now uses the test statistic. There is a new plot option to help visualizing the result of running the goodness command. The command plot goodness after running goodness generates a histogram plot of the test statistic values from the simulations as well as a line indicating the value for the actual data.

The new test statistics added include the familiar Kolmogorov-Smirnov along with the related Anderson-Darling and Cramer-von Mises. The runs statistic tests for runs of consecutive positive or negative residuals. Also added is the original Pearson chi-square statistic for Poisson data.

One new fit statistic has been added: the Whittle statistic for fitting power density spectra.

All these statistics are described in the manual in the extended appendix on statistics.

Wednesday, January 02, 2013

XSPEC v12.8 - Parallelization

XSPEC v12.8 is now available as part of HEAsoft 6.13.

This is the first of several posts on new features in v12.8. The parallel command tells XSPEC to use multiple cores on your machine. At the moment we have parallel processing in two places: in the fit command and in the error command. The standard Levenberg-Marquardt algorithm requires the calculation of the numerical derivative of the statistic with respect to each of the parameters. This is an obvious place to split the calculation among multiple cores and we do this if the command parallel leven N is given, where N is the number of cores to use. The more compute-intensive the model, the more time will be saved. In general however, this does not produce improvements linear in the number of cores.

The error command can be run on multiple parameters at the same time and using the parallel command here with parallel error N does produce a linear gain in N provided there are at least N parameter errors being calculated. The one current drawback with the parallel error calculation comes if the error finds a new minimum. At present all that will happen is that XSPEC will write out that a new error has been found but not update the best fit. To actually find the new best fit the error command must be repeated on individual parameters. We do intend to fix this.

The next command for which we intend to add parallel processing is steppar. We welcome other suggestions and feedback.