Changing the grouping scheme of a data set within Sherpa

Sherpa Threads (CIAO 4.4 Sherpa v1)

Overview

Synopsis:

In order to use Gaussian statistics to fit a model to a data set, it is often necessary to "group" the data - i.e., combine channels until you have enough counts - before use. It is possible to set and change the grouping of a file after it has been read into Sherpa by using the group commands: set_grouping, group, group_counts, group_snr, group_adapt, group_adapt_snr, group_bins, and group_width.

This thread shows how you can use these functions when fitting PHA data.

Last Update: 15 Dec 2011 - reviewed for CIAO 4.4: added a warning about filtering/grouping source and background data sets

Getting Started
- Loading the data
Grouping the data
Fitting the data
Scripting It
Notes
Summary
History
Images
- Figure 1: Data grouped by 15 and 30 counts per group
- Figure 2: Fit to both data sets

Getting Started

Loading the data

This thread uses the same data set as used in the Introduction to Fitting PHA Spectra thread. We load the data set twice - using load_pha - so that we can easily see the effect of changing the grouping scheme (the screen output has been omitted for clarity):

sherpa> load_pha(1, "3c273.pi")
...
sherpa> load_pha(2, "3c273.pi")
...

We note that the data in 3c273.pi was created so that each group contained at least 15 counts - by using the NUM_CTS grouptype option of dmgroup - as can be seen by using the dmhistory tool.

sherpa> !dmhistory 3c273.pi dmgroup | tr ' ' "\012"
dmgroup
infile="3c273.pi"
outfile="./3c273.tmp"
grouptype="NUM_CTS"
grouptypeval="15"
binspec=""
xcolumn="channel"
ycolumn="counts"
tabspec=""
tabcolumn=""
stopspec=""
stopcolumn=""
errcolumn=""
clobber="no"
verbose="0" 
maxlength="0"

The command was preceded by "!" to tell Sherpa to execute it as a shell command, and "| tr ' ' "\012"" was used to add a new-line character between parameter values (to make the screen output easier to read).

Grouping the data

We use the group_counts function to group the second data set by 30 counts per group, which is double that of the first data set. The "WARNING" message can be ignored as no filters have been applied to the data. The two data sets are then plotted using logarithmic scaling on both axes:

sherpa> group_counts(2, 30)
WARNING: grouping flags have changed, noticing all bins

sherpa> set_xlog()
sherpa> set_ylog()

sherpa> plot("data", 1, "data", 2)

The resulting plot (Figure 1) shows how the data looks before and after re-grouping.

Figure 1: Data grouped by 15 and 30 counts per group

[Plot of data set 1 grouped by 15 counts per bin, and data set 2 grouped by 30 counts per bin]

[Print media version: Plot of data set 1 grouped by 15 counts per bin, and data set 2 grouped by 30 counts per bin]

Figure 1: Data grouped by 15 and 30 counts per group

The top plot shows the data as read in to Sherpa - which is binned by 15 counts per group - and the bottom plot shows the data set after group_counts has been called to re-group the data to 30 counts per group.

In this thread, we use the group_counts command as an example. The other available grouping commands are:

group_bins() - this function divides the channels into the specified number ("num") of bins.
group_width() - this function divides the channels such that there are "num" bins in each group.
group_snr() - this function allows the user to adaptively group PHA spectral data by signal-to-noise ratio, i.e., group bins of data until each group exceeds at least the minimum specified signal-to-noise ratio.
group_adapt() - this function allows the user to adaptively group PHA spectral data by counts, i.e., group bins of data until the number of counts in each group exceeds the minimum number of counts specified in the 'min' argument, keeping bright features ungrouped while grouping low signal-to-noise regions.
group_adapt_snr() - this function allows the user to adaptively group PHA spectral data by signal-to-noise ratio, i.e., group bins of data until each group exceeds at least the minimum specified signal-to-noise ratio.

Warning about filtering and grouping source and background data:

If a source data set and its associated background data are to both be grouped with a Sherpa analysis session (not considered in this thread), then due to a bug in CIAO 4.4, you must ensure that the background data is grouped before the source, as shown in the example below; otherwise, any data filter in place for the background will be removed after the grouping command is issued.

sherpa> notice(0.5,7.)             # source and background filtered to
                                   # the 0.5-7. keV energy range
 

sherpa> group_counts(20, bkg_id=1) # background data must be grouped before
                                   # source data to retain background filter                            

sherpa> group_counts(20)           # both source and background filtered and gro

Fitting the data

Once the data has been re-grouped you can use it just like any other data set. Here we repeat the fit made in the Introduction to fitting PHA spectra to see what difference the alternate grouping scheme makes.

sherpa> notice_id(2, 0.1, 6.0)
sherpa> subtract(2)
sherpa> set_source(2, xsphabs.abs1 * powlaw1d.p1)
sherpa> abs1.nh = 0.07
sherpa> freeze(abs1)
sherpa> guess(2, p1)
sherpa> fit(2)
 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
Data Set               = 2
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 7.58076e+10
Final fit statistic   = 28.0675 at function evaluation 25
Data points           = 22
Degrees of freedom    = 20
Probability [Q-value] = 0.107811
Reduced statistic     = 1.40338
Change in statistic   = 7.58076e+10
   p1.gamma       2.1615      
   p1.ampl        0.000232546 

sherpa> show_fit()
Optimization Method: LevMar
name    = levmar
ftol    = 1.19209289551e-07
xtol    = 1.19209289551e-07
gtol    = 1.19209289551e-07
maxfev  = None
epsfcn  = 1.19209289551e-07
factor  = 100.0
verbose = 0

Statistic: Chi2Gehrels

Fit:Data Set               = 2
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 7.58076e+10
Final fit statistic   = 28.0675 at function evaluation 25
Data points           = 22
Degrees of freedom    = 20
Probability [Q-value] = 0.107811
Reduced statistic     = 1.40338
Change in statistic   = 7.58076e+10
   p1.gamma       2.1615      
   p1.ampl        0.000232546 

sherpa>  covar(2)
Data Set               = 2
Confidence Method     = covariance
Fitting Method        = levmar
Statistic             = chi2gehrels
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.gamma           2.1615   -0.0763938    0.0763938
   p1.ampl       0.000232546 -1.40839e-05  1.40839e-05

which can be compared to the original results:

sherpa> notice_id(1, 0.1, 6.0)
sherpa> subtract(1)
sherpa> set_source(1, abs1 * p1)
sherpa> guess(p1)
sherpa> fit(1)
Data Set               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 38.3145
Final fit statistic   = 37.9079 at function evaluation 18
Data points           = 44
Degrees of freedom    = 42
Probability [Q-value] = 0.651155
Reduced statistic     = 0.902569
Change in statistic   = 0.406558
   p1.gamma       2.15853     
   p1.ampl        0.000224842 
  

sherpa> covar(1)
Data Set               = 1
Confidence Method     = covariance
Fitting Method        = levmar
Statistic             = chi2gehrels
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.gamma          2.15853   -0.0827859    0.0827859
   p1.ampl       0.000224842 -1.48255e-05  1.48255e-05

The following shows both fits (the extra commands are used to make the two plots have the same X axis, and to add plot labels):

sherpa> plot_fit()
sherpa> plot_fit(2, overplot=True)

sherpa> current_plot("plot1")
sherpa> reposition_plot(0.15, 0.475, 0.9, 0.9)
sherpa> current_plot("plot2")
sherpa> reposition_plot(0.15, 0.1 ,0.9 ,0.475)

sherpa> current_plot("plot1")
sherpa> set_xaxis(["ticklabel.visible",False])

sherpa> current_plot("all")
sherpa> limits(X_AXIS, 0.1, 7.0)

sherpa> current_plot("plot1")
sherpa> set_plot_xlabel("")

sherpa> current_plot("plot2")
sherpa> set_plot_title("")

sherpa> current_plot("plot1")
sherpa> add_label(0.5, 0.0006, "15 counts per group")
sherpa> set_label(["color","green"])
sherpa> current_plot("plot2")
sherpa> add_label(0.5, 0.0006, "30 counts per group")
sherpa> set_label(["color","green"])

which creates this plot (Figure 2).

Figure 2: Fit to both data sets

[Plot of fits to data set 1 grouped by 15 counts per bin, and data set 2 grouped by 30 counts per bin]

[Print media version: Plot of fits to data set 1 grouped by 15 counts per bin, and data set 2 grouped by 30 counts per bin]

Figure 2: Fit to both data sets

This plot shows the fits to the data when grouped by 15 counts per group (top plot) and 30 counts per group (bottom plot).

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") 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, etcetera.)

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.3 syntax to you.

Notes

The group_counts and related functions can also be used on data that was not grouped before being read into Sherpa; this can be useful for two reasons:

You can use the functions to find the best grouping scheme for your data without having to re-run the dmgroup tool and re-load the data into Sherpa.
You can fit the un-grouped data with the Cash statistic and then use the functions to make it easier to compare the fit to the data in plots.

Summary

This thread shows how you can use the grouping functionality in Sherpa 4.4 to change the grouping scheme of a PHA file once it has been read into Sherpa. This allows you to see how sensitive the fit results are to the grouping scheme by changing the number of counts per group or using a different method for grouping the data.

History

14 Dec 2004	updated for CIAO 3.2: script version and path
17 Jun 2005	updated information in Get Started on loading the script
21 Dec 2005	reviewed for CIAO 3.3: no changes
01 Dec 2006	reviewed for CIAO 3.4: no changes
14 Dec 2008	updated for CIAO 4.1: `sherpa_utils.sl` script replaced by new Sherpa 4.1 grouping functionality
29 Apr 2009	new `script` command is available with CIAO 4.1.2
17 Dec 2009	updated for CIAO 4.2: new `group_bins` and `group_width` commands
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: new `set_xlog` and `set_ylog` commands are available for setting the axis scale of plots to logarithmic
15 Dec 2011	reviewed for CIAO 4.4: added a warning about filtering/grouping source and background data sets