Loving ISIS - Confessions of a Former XSPEC User
Fitting
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:
constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bkn2pow(1)+gaussian(1)) 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.9625612The 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.
