Loving ISIS - Confessions of a Former XSPEC User



We've loaded the data, added systematic errors to the PCA data, grouped the PCA and HEXTE data, extended the model grids to accommodate the reflect model, defined the model, and set its parameters. Now let's actually fit the data.

ISIS provides multiple minimization methods: mpfit and lmdif (both the same levenburg-marquardt method; the former has stricter tolerances and hence may take longer to converge than the latter), marquardt (another l-m method), minim (a simplex method), and subplex. The ISIS command set_fit_method(); can be used to toggle among them, as well as change the control parameters used when they are run (e.g., maximum number of function evaluations). I have created aliases: mpfit, lmdif, marq, minim, and subplex, to toggle among them (albeit without changing any of their internal parameters). In practical terms, lmdif and marquardt are faster than minim and subplex. I have found the two most useful ones to be lmdif and subplex. The former is fast and reliable, but the latter, albeit slow, can find that odd, interesting corner of parameter space that other methods overlook. (For this reason, subplex is also very good at crashing poorly written model function code. This will happen from time to time with some especially complex local model codes. Unfortunately, if this happens inside of Fortran code, this will crash your whole ISIS session.) lmdif is the default fitting method, and it generally works very well.

Unlike XSPEC, ISIS will not automatically evaluate a model upon its definition. To do so, use the ISIS command eval_counts();. The ISIS function renorm_counts(); provides an analog to the XSPEC command renorm, while the ISIS function fit_counts(); (which I have aliased to fit();) is an analog to the XSPEC command fit. The above ISIS functions will return 0 upon successful completion (-1 if they fail). These return values can be very useful when invoking a fit command internally to a script. Both fit_counts and eval_counts can be called with an optional variable argument, e.g.,

     isis> () = fit_counts(&info);
which will return a structure variable, info, that has useful information, such as the value of the fit statistic (e.g., chi squared). From the analysis script, here is the ISIS fit to our data, using a reflected, broken power law model (plus relativistic line!):


Error Bars

ISIS provides several functions, e.g., conf(); and vconf();, that are the functional equivalent of the XSPEC err command. As many of us have experienced, we often stumble across new, lower chi squared solutions while doing error searches. To that end, ISIS provides:

     (pmin,pmax) = conf_loop(params[],level,tolerance;save,prefix="parameter_file.");
which will loop through a set of parameters, calculating error bars and refitting the model if a new lower chi^2 is found. The new parameters and the final error bars are saved to files if the save and prefix qualifiers are set. This function is "transparently" parallelized. It will automatically perform the error bar search using multiple cores on any multi-core machine. The number of cores, and the niceness level of the spawned processes, are controllable, e.g.,.
     (pmin,pmax) = conf_loop(params[],level,tolerance;save,prefix="file.",nice=19,num_slaves=2);

Here is an example output for the doubly broken power law plus gaussian line model from the analysis script:

 idx  param            tie-to  freeze         value         min         max
  1  constant(1).factor    0     0         1.178409  0.01841264    12108.66
  2  phabs(1).nH           0     1              0.6           0      100000  10^22
  3  highecut(1).cutoffE   0     0         41.41828    35.72186    49.18991  keV
  4  highecut(1).foldE     0     0          144.095    125.5314     166.222  keV
  5  bkn2pow(1).norm       0     0          183.901  0.09618433       10000
  6  bkn2pow(1).PhoIndx1   0     0        0.8327529   0.7298967   0.9163523
  7  bkn2pow(1).BreakE1    0     0     0.0004672722  1.097883e-06    3.152033  keV
  8  bkn2pow(1).PhoIndx2   0     0          1.68408    1.677428    1.691277
  9  bkn2pow(1).BreakE2    0     0         11.05349    10.48958    11.56851  keV
 10  bkn2pow(1).PhoIndx3   0     0         1.477259    1.457069     1.49763
 11  gaussian(1).norm      0     0      0.001813016  0.001595308  0.002065615
 12  gaussian(1).LineE     0     1              6.4           0     1000000  keV
 13  gaussian(1).Sigma     0     0         0.660237   0.5403519   0.7872875  keV
 14  constant(2).factor    0     1                1           0       1e+10
 15  constant(3).factor    0     0        0.9498044   0.9373023   0.9625612
The min/max values give the 90% confidence level ranges of the parameters.

Confidence Contours

As in XSPEC, ISIS lets you do confidence contours on pairs of parameters, using the ISIS commands:

     x = conf_grid(idx,lo,hi,nstep);    % Create x-axis grid
     y = conf_grid(idy,lo,hi,nstep);    % Create y-axis grid
     cntr = conf_map_counts(x,y);       % Make the error contours
     plot_conf(cntr);                   % Plot the error contours
     save_conf(cntr,"file");            % Save the error contours
ISIS has a couple of especially nice features here. First, you can save the contours to a FITS file (and then later reload them). Second, you can define a 'failure hook' to let the contour mapping program know what to do if it fails to resolve a fit at any given point. I have defined conf_map(); to toggle between the lmdif and subplex fit methods, should a fit fail. (It switches to whichever one was not used, and then back again after that fit point.) Results for error contours of the low and high X-ray energy power laws in the doubly broken power law plus gaussian model are shown here:


Next up: Fluxes and Equivalent Widths.

This page was last updated Oct 7, 2013 by Michael Nowak. To comment on it or the material presented here, send email to mnowak@space.mit.edu.
Valid HTML 4.01! Made with JED. Viewable With Any Browser.