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 four fit minimization methods: lmdif (a
levenburg-marquardt method), 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:
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
also have found lmdif and minim to be the most robust,
although subplex has occasionally been useful for very noisy
data and/or models. 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,
I have written two functions:
(pmin,pmax) = error_loop(a,nmax,file);
error_save(a,pmin,pmax,file.err);
constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bkn2pow(1)+gaussian(1))
ITER PARI P_VALUE P_MIN P_MAX TIME
0 3 4.345e+01 3.880e+01 5.068e+01 Thu Aug 25 20:17:27 2005
0 4 1.309e+02 1.151e+02 1.464e+02 Thu Aug 25 20:17:53 2005
0 5 2.285e+02 5.051e+01 1.272e+03 Thu Aug 25 20:18:54 2005
0 6 8.297e-01 7.834e-01 9.171e-01 Thu Aug 25 20:20:32 2005
chi^2 = 7.89e+01, DoF = 87, time = Thu Aug 25 20:21:06 2005
1 7 3.755e-04 1.405e-04 1.089e-03 Thu Aug 25 20:22:00 2005
1 8 1.684e+00 1.677e+00 1.691e+00 Thu Aug 25 20:23:00 2005
1 9 1.105e+01 1.048e+01 1.157e+01 Thu Aug 25 20:23:31 2005
1 10 1.478e+00 1.458e+00 1.498e+00 Thu Aug 25 20:24:44 2005
1 11 1.817e-03 1.597e-03 2.068e-03 Thu Aug 25 20:25:16 2005
1 13 6.616e-01 5.406e-01 7.876e-01 Thu Aug 25 20:25:38 2005
1 15 9.487e-01 9.371e-01 9.603e-01 Thu Aug 25 20:26:18 2005
1 3 4.350e+01 3.880e+01 5.069e+01 Thu Aug 25 20:26:40 2005
1 4 1.308e+02 1.151e+02 1.464e+02 Thu Aug 25 20:27:07 2005
1 5 2.261e+02 5.053e+01 1.271e+03 Thu Aug 25 20:28:13 2005
1 6 8.303e-01 7.840e-01 9.171e-01 Thu Aug 25 20:29:55 2005
(The times are output, since for some very computationally expensive
fit models, it's nice to keep track of the speed of progress.) The
final error bar output file for the above is:
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 1 1 0 1e+10 2 phabs(1).nH 0 1 0.6 0 100000 10^22 3 highecut(1).cutoffE 0 0 43.49888 38.79909 50.68562 keV 4 highecut(1).foldE 0 0 130.795 115.0543 146.4468 keV 5 bkn2pow(1).norm 0 0 226.1095 50.52711 1271.178 6 bkn2pow(1).PhoIndx1 0 0 0.8302748 0.7839841 0.9170555 7 bkn2pow(1).BreakE1 0 0 0.0003754774 0.0001404727 0.001089473 keV 8 bkn2pow(1).PhoIndx2 0 0 1.684176 1.677435 1.691278 9 bkn2pow(1).BreakE2 0 0 11.04503 10.48013 11.57458 keV 10 bkn2pow(1).PhoIndx3 0 0 1.478041 1.457901 1.498382 11 gaussian(1).norm 0 0 0.001817382 0.001596692 0.002067671 12 gaussian(1).LineE 0 1 6.4 0 1000000 keV 13 gaussian(1).Sigma 0 0 0.6616106 0.5406061 0.7876207 keV 14 constant(2).factor 0 1 1 0 1e+10 15 constant(3).factor 0 0 0.9487277 0.937108 0.9603426
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 minim 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:
Now that we've done all our analysis, let's talk more about plotting our results.
