Simultaneously Fitting Two Data Sets
Sherpa Threads (CIAO 4.9 Sherpa v1)
The Sherpa syntax easily allows the modeling of multiple, independent data sets. In this example, we have a source observed on two different occasions. If we assume that the source is constant, then the only change will be in the instrument response. We account for this by defining each spectrum independently and fitting with the same model. This allows us to increase the significance of our fit, without making assumptions inherent in the process of averaging spectra and responses.
This thread may also be used to fit source and background data which are independently grouped. Note that the "standard method" of doing this is described in the Simultaneously Fitting Source and Background Spectra thread.
While the sample data files used in this thread are available as sherpa.tar.gz, note also that they are generated with the Extract ACIS Spectra and Response Files for Pointlike Sources thread.
Last Update: 8 Nov 2016 - reviewed for CIAO 4.9, no content change.
- Statistical Issues: Background Subtraction
- Reading FITS Data
- Defining Instrument Responses
- Defining Source Models
- Fitting the Data
- Link Parameters and Fit Again
- Scripting It
A typical dataset may contain multiple spectra, one of which contains contributions from the source and the background, and one or more which contain background counts alone. (The background itself may contain contributions from the cosmic X-ray background, the particle background, and so on, but we ignore this complication.)
The proper way to treat background data is to model them. However, many X-ray astronomers subtract background data from the raw data.
Why should one not subtract background?
- It reduces the amount of statistical information in the analysis—the final fit parameter values will be a less accurate estimate of the true values.
- The background subtracted data are not Poisson-distributed; one cannot fit them with the Poisson likelihood or the Cash statistic, even in the low-counts limit. For example, subtracting a background can give negative counts; this is definitely not Poissonian!
- Fluctuations, particularly in the vicinity of localized features, can adversely affect analysis.
Simultaneous fitting is identical to the basic method of fitting data, except for the specification of the data set ID. Most Sherpa commands have the option of specifying the ID of the data set on which to invoke the command. By default, this number is assumed to be 1; the ID is assigned when the data are loaded into Sherpa.
The spectral data from the two observations of the source are stored in the FITS files pi2278.fits and pi2286.fits. Data Set pi2278.fits is input with the load_pha() command and given ID "1":
sherpa> load_pha(1, "pi2278.fits") statistical errors were found in file 'pi2278.fits' but not used; to use them, re-read with use_errors=True
The second data set is loaded similarly, but is given ID "2":
sherpa> load_pha(2, "pi2286.fits") statistical errors were found in file 'pi2286.fits' but not used; to use them, re-read with use_errors=True
Sherpa will automatically read in any instrument response files that are defined in the header of the spectrum (keywords RESPFILE and ANCRFILE).
sherpa> load_arf(1, "arf2278.fits") sherpa> load_rmf(1, "rmf2278.fits") sherpa> load_arf(2, "arf2286.fits") sherpa> load_rmf(2, "rmf2286.fits")
We want to model each of the spectra with an absorbed power-law, so a power-law component (powlaw1d, called pl1) and an absorbing component (xsphabs, called abs1) are defined. The same model expression is assigned to each data set:
We can use the Sherpa guess command to guess the initial parameter values and ranges for the power-law model component. To have Sherpa automatically query for the initial parameter values when a model is established, set paramprompt(True) (it is False by default).
sherpa> guess(pl1) sherpa> show_model() Model: 1 apply_rmf(apply_arf((11619.081430486 * (xsphabs.abs1 * powlaw1d.pl1)))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- abs1.nH thawed 1 0 100000 10^22 atoms / cm^2 pl1.gamma thawed 1 -10 10 pl1.ref frozen 1 -3.40282e+38 3.40282e+38 pl1.ampl thawed 6.4442e-06 6.4442e-09 0.0064442 Model: 2 apply_rmf(apply_arf((11619.081410727 * (xsphabs.abs1 * powlaw1d.pl1)))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- abs1.nH thawed 1 0 100000 10^22 atoms / cm^2 pl1.gamma thawed 1 -10 10 pl1.ref frozen 1 -3.40282e+38 3.40282e+38 pl1.ampl thawed 6.4442e-06 6.4442e-09 0.0064442
The guess command makes an initial guess at parameter values to ensure convergence, but it is always a good idea to check that the initial range of values (soft limits) is sensible for the data being fit.
Data above 5.5 keV is ignored in both data sets to cut out energies with no data. While this is not strictly necessary for this data, it does produce a nicer plot. Ignoring the higher energies may be especially useful for ACIS-S data, where the high-energy background is stronger.
The comma after the value indicates an open range, i.e. everything above "5.5" will be ignored:
In order to fit both data sets simultaneously, fit() is called with both IDs (it can also be called with no arguments to fit all data sets defined in the session):
sherpa> fit(1,2) Datasets = 1, 2 Method = levmar Statistic = chi2gehrels Initial fit statistic = 28.7326 Final fit statistic = 5.72177 at function evaluation 47 Data points = 14 Degrees of freedom = 11 Probability [Q-value] = 0.891269 Reduced statistic = 0.520161 Change in statistic = 23.0109 abs1.nH 1.12165 pl1.gamma 2.23491 pl1.ampl 4.27375e-05
The calc_stat_info command (see also get_stat_info—similar to the CIAO 3.4 GOODNESS command—may be used to access goodness-of-fit statistics for each individual data set included in the simultaneous fit, plus the simultaneous fit, without having to re-run the fit.
sherpa> calc_stat_info() Dataset = 1 Statistic = chi2gehrels Fit statistic value = 3.62858 Data points = 7 Degrees of freedom = 4 Probability [Q-value] = 0.458599 Reduced statistic = 0.907145 Dataset = 2 Statistic = chi2gehrels Fit statistic value = 2.09319 Data points = 7 Degrees of freedom = 4 Probability [Q-value] = 0.718624 Reduced statistic = 0.523297 Datasets = [1, 2] Statistic = chi2gehrels Fit statistic value = 5.72177 Data points = 14 Degrees of freedom = 11 Probability [Q-value] = 0.891269 Reduced statistic = 0.520161
sherpa> plot_fit() sherpa> plot_fit(2, overplot=True) # or simply 'plot("fit", 1, "fit", 2)' sherpa> print_window("two_fits.ps")
If the second spectrum is similar, but not identical to the first spectrum, a different model component could be used in one of the source definitions. Additionally, parameter linking can be used to fix parameter values to each other.
For example, if the intensity of the source is expected to change while the power-law slope does not, we could define a new power-law component and link the slopes together. The result is another free parameter in the fit: the normalization of the second power-law component.
Now fit the data again:
sherpa> fit(1,2) Datasets = 1, 2 Method = levmar Statistic = chi2gehrels Initial fit statistic = 1.68851e+10 Final fit statistic = 5.71885 at function evaluation 31 Data points = 14 Degrees of freedom = 10 Probability [Q-value] = 0.838305 Reduced statistic = 0.571885 Change in statistic = 1.68851e+10 abs1.nH 1.11887 pl1.gamma 2.23165 pl1.ampl 4.22737e-05 pl2.ampl 4.28572e-05
Only the amplitude of the new component (pl2.ampl) is fit. In this case, the amplitudes of the two power-laws are nearly identical. The final model parameter values may be printed:
sherpa> show_model() Model: 1 apply_rmf(apply_arf((11619.081430486 * (xsphabs.abs1 * powlaw1d.pl1)))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- abs1.nH thawed 1.11887 0 100000 10^22 atoms / cm^2 pl1.gamma thawed 2.23165 -10 10 pl1.ref frozen 1 -3.40282e+38 3.40282e+38 pl1.ampl thawed 4.22737e-05 6.4442e-09 0.0064442 Model: 2 apply_rmf(apply_arf((11619.081410727 * (xsphabs.abs1 * powlaw1d.pl2)))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- abs1.nH thawed 1.11887 0 100000 10^22 atoms / cm^2 pl2.gamma linked 2.23165 expr: pl1.gamma pl2.ref frozen 1 -3.40282e+38 3.40282e+38 pl2.ampl thawed 4.28572e-05 0 3.40282e+38
This thread is complete, so we can exit the Sherpa session:
The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing exec(open("fit.py").read()) on the Sherpa command line.
The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:
sherpa> script(filename="sherpa.log", clobber=False)
(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, et cetera.)
The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.9 syntax to you.
|14 Nov 2007||rewritten for CIAO 4.0 Beta 3|
|09 Dec 2008||figures moved inline with text|
|11 Dec 2008||updated for Sherpa 4.1|
|16 Feb 2009||example of guess functionality added|
|29 Apr 2009||new script command is available with CIAO 4.1.2|
|17 Dec 2009||updated for CIAO 4.2|
|19 Mar 2010||photoelectric absorption model xswabs replaced with xsphabs|
|13 Jul 2010||updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.|
|15 Dec 2010||updated for CIAO 4.3: calc_stat_info is available for accessing goodness-of-fit statistics for each individual data set included in a simultaneous fit|
|15 Dec 2011||added an additional option for simultaneously plotting fits of two different data sets, using plot_fit(); reviewed for CIAO 4.4 (no changes)|
|03 Dec 2013||reviewed for CIAO 4.6: no changes|
|30 Jan 2015||updated for CIAO 4.7: no content change|
|02 Dec 2015||updated for CIAO 4.8: cleaned up description on using guess.|
|08 Nov 2016||reviewed for CIAO 4.9, no content change.|