Step-by-Step Guide to Estimating Errors and Confidence Levels

Getting Started

To obtain the sample data files used in this thread, download the sherpa.tar.gz file as described in the "Sherpa Getting Started" thread.

Finding the best fit

First, we check the current Sherpa settings. By default, Sherpa uses the Levenberg-Marquardt method to optimize the fit of a model to data with chi-squared statistics, using the Gehrels variance function. In this example, we choose to fit with the default optimization method, and change the fit statistic from default to the Chi2Xspecvar statistic. (The available optimization methods and statistics may be listed with list_methods and list_stats, and can be changed with set_method and set_stat. Detailed information on each option may be found in the ahelp documentation, as well as on the Sherpa Optimization Methods and Statistics pages.)

sherpa> clean()
sherpa> show_method()
sherpa> show_stat()
sherpa> set_stat("chi2xspecvar")

sherpa> load_data("source_grouped_pi.fits")
sherpa> notice(0.5,8.0)
sherpa> plot_data()
sherpa> log_scale(XY_AXIS)

sherpa> set_source(xswabs.abs1*powlaw1d.p1)
sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2xspecvar
Initial fit statistic = 6.99384e+09
Final fit statistic   = 124.032 at function evaluation 33
Data points           = 133
Degrees of freedom    = 130
Probability [Q-value] = 0.630974
Reduced statistic     = 0.95409
Change in statistic   = 6.99384e+09
   abs1.nH        2.36636     
   p1.gamma       1.49037     
   p1.ampl        0.00190285  

The data we will be fitting is read into the session with the load_data command, which automatically loads the instrument response associated with a source data set when the ARF and RMF filenames are recorded in the header of the source data file, as shown above; this also applies to background data associated with the source (not considered in this example). The plot in Figure 1 results from the plot_data command above, showing the source data that is to be fit, between 0.5-8.0 keV.

Figure 1: Source spectrum

Figure 1: Source spectrum

The source model is set to an absorbed power-law with the set_source command , and fit with the fit command, producing the results shown above. The fit and delta-chi residuals may be plotted with the plot or plot_fit_delchi command.

sherpa> plot("fit",1,"delchi",1)
sherpa> current_plot("plot1")
sherpa> log_scale(XY_AXIS)
sherpa> current_plot("plot2")
sherpa> log_scale(X_AXIS)
sherpa> limits(X_AXIS,1,10)

The resulting plot is shown in Figure 2.

Figure 2: Best fit model with residuals

Figure 2: Best fit model with residuals

The show_fit and get_fit_results commands are also available for checking the quality of the fit:

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: Chi2XspecVar
Chi Squared with data variance (XSPEC style)

Fit:Dataset               = 1
Method                = levmar
Statistic             = chi2xspecvar
Initial fit statistic = 6.99384e+09
Final fit statistic   = 124.032 at function evaluation 33
Data points           = 133
Degrees of freedom    = 130
Probability [Q-value] = 0.630974
Reduced statistic     = 0.95409
Change in statistic   = 6.99384e+09
   abs1.nH        2.36636     
   p1.gamma       1.49037     
   p1.ampl        0.00190285  

Confidence limits for individual parameters

After finding the best-fit model parameter values, we calculate the confidence limits (parameter bounds) for these parameters using either the confidence, projection, or covariance methods. The confidence limits are defined by a required confidence level in our analysis, typically the 68.3% or 90% level, corresponding to 1 sigma and 1.6 sigma for the normal distribution. The table below displays the relationship between sigma and confidence level for one significant parameter.

Confidence intervals for a normal distribution:

We will use the confidence method to estimate 1 sigma errors on the gamma parameter of the power-law model component:

sherpa> print get_conf()
name         = confidence
numcores     = 4
verbose      = False
openinterval = False
max_rstat    = 3
maxiters     = 200
soft_limits  = False
eps          = 0.01
fast         = False
maxfits      = 5
remin        = 0.01
tol          = 0.2
sigma        = 1
parallel     = True

sherpa> conf(p1.gamma)
p1.gamma lower bound:   -0.0840758
p1.gamma upper bound:   0.085951
Dataset               = 1
Confidence Method     = confidence
Fitting Method        = levmar
Statistic             = chi2xspecvar
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.gamma          1.49037   -0.0840758     0.085951

To access the 90% confidence limits on a parameter, the sigma field of the conf_opt variable should be changed to 1.6.

sherpa> set_conf_opt("sigma",1.6)

sherpa> conf(p1.gamma, abs1.nH)
abs1.nH lower bound:	-0.192986
p1.gamma lower bound:	-0.133226
p1.gamma upper bound:	0.138241
abs1.nH upper bound:	0.207404
Dataset               = 1
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2xspecvar
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.gamma          1.49036    -0.133226     0.138241
   abs1.nH           2.36635    -0.192986     0.207404

We have also used conf to calculate the uncertainty on the optimized hydrogen column density parameter of the absorption model, abs1.nH.

To estimate errors on all the thawed parameters, conf should be called with no parameter names. Since sigma is still set to 1.6 for the confidence method, the following calculates the 90% confidence limits for all the thawed parameters:

sherpa> conf()
abs1.nH lower bound:	-0.192986
p1.ampl lower bound:	-0.000341981
p1.gamma lower bound:	-0.133226
abs1.nH upper bound:	0.207404
p1.gamma upper bound:	0.138241
p1.ampl upper bound:	0.000430029
Dataset               = 1
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2xspecvar
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nH           2.36635    -0.192986     0.207404
   p1.gamma          1.49036    -0.133226     0.138241
   p1.ampl        0.00190283 -0.000341981  0.000430029

The covariance command behaves similarly to conf, although the fields in the state object are different. While it is quicker than the confidence method, it is less accurate; i.e., it is always symmetric because it uses the diagonal elements of the covariance matrix and ignores correlations between the parameters.

Note that the computationally intensive confidence function has been parallelized as of Sherpa in CIAO 4.2, to make use of multi-core systems (i.e., laptops or desktops with 2 or 4 cores). Users with multi-core machines should see significant improvements in efficiency compared to previous releases of Sherpa.

sherpa> conf()
   ...

WARNING: hard minimum hit for parameter bpow1.gamma2
WARNING: hard maximum hit for parameter bpow1.gamma2
WARNING: hard minimum hit for parameter bpow1.eb
WARNING: hard maximum hit for parameter bpow1.eb
   ...

Dataset               = 1
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = neldermead
Statistic             = cstat
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.gamma1       1.54147   -0.0292891   0.0292709
   bpow1.gamma2       8.10056        -----       -----
   bpow1.eb           9.49083       -----        -----
   bpow1.ampl         0.022806 -0.000378395  0.000383854

sherpa>  conf(g15.Sigma)
    ...
g15.Sigma -: WARNING: The confidence level lies within (8.706380e-05,
9.252185e-05)
    ...
Datasets              = 1, 2
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2datavar
confidence 1.64-sigma (89.8995%) bounds:
  Param            Best-Fit  Lower Bound  Upper Bound
  -----            --------  -----------  -----------
  g15.Sigma     0.000997626 -0.000907834  0.000597058

sherpa> set_conf_opt('verbose',1)
sherpa> conf(g15.Sigma)
#
# f[  2.931742e+00   2.957942e-02   1.471941e+00   9.976265e-04
1.851972e-04
  1.840837e+00   3.667986e-03   7.820442e-05   1.864564e+00   2.562143e-03
  1.415699e-04   2.004898e+00   6.288115e-04   1.259512e-04] =
8.569929e+02
# sigma = 1.640000e+00
# target_stat = 8.596825e+02
# tol = 1.000000e-02
# smin = [-2.    0.    1.45  0.    0.    1.82  0.    0.    1.85  0.
0.    1.96
 0.    0.  ]
# smax = [  9.000000e+00   1.000000e+24   1.500000e+00   1.000000e-02
1.000000e-03
  1.850000e+00   1.000000e-02   1.000000e-03   1.900000e+00   1.000000e-02
  1.000000e-03   2.040000e+00   1.000000e-03   1.000000e-03]
# hmin = [ -3.402823e+38   0.000000e+00   0.000000e+00   0.000000e+00
0.000000e+00
  0.000000e+00   0.000000e+00   0.000000e+00   0.000000e+00   0.000000e+00
  0.000000e+00   0.000000e+00   0.000000e+00   0.000000e+00]
# hmax = [  3.402823e+38   3.402823e+38   3.402823e+38   3.402823e+38
3.402823e+38
  3.402823e+38   3.402823e+38   3.402823e+38   3.402823e+38   3.402823e+38
  3.402823e+38   3.402823e+38   3.402823e+38   3.402823e+38]
#
# Note: for the intermediate steps, the notation:
        par.name -/+: f( x ) = stat
# ==> `stat` is the statistic when parameter `par.name` is frozen at `x`
# while searching for the `lower/upper` confidence level, repectively.
#
g15.Sigma -: f( 4.812517e-04 ) = -2.565633e+00
g15.Sigma -: f( 0.000000e+00 ) = 3.505966e+01
g15.Sigma -: f( 2.406259e-04 ) = -2.565434e+00
g15.Sigma -: f( 2.113694e-04 ) = -2.565393e+00
g15.Sigma -: f( 1.653199e-04 ) = -2.565381e+00
g15.Sigma -: f( 1.343754e-04 ) = -2.564794e+00
g15.Sigma -: f( 1.079722e-04 ) = -2.565100e+00
g15.Sigma -: f( 8.706380e-05 ) = 3.778055e+01
g15.Sigma -: f( 1.068516e-04 ) = -2.563338e+00
g15.Sigma -: f( 1.022595e-04 ) = -2.562500e+00
g15.Sigma -: f( 9.961908e-05 ) = -2.563027e+00
g15.Sigma -: f( 9.721570e-05 ) = -2.564524e+00
g15.Sigma -: f( 9.532733e-05 ) = -2.562276e+00
g15.Sigma -: f( 9.377811e-05 ) = -2.562519e+00
g15.Sigma -: f( 9.252185e-05 ) = -2.563134e+00
g15.Sigma -: f( 9.149979e-05 ) = -2.563034e+00
g15.Sigma -: f( 9.066944e-05 ) = -2.563196e+00
g15.Sigma -: WARNING: The confidence level lies within (8.706380e-05,
9.252185e-05)
g15.Sigma lower bound:  -0.000907834
g15.Sigma +: f( 1.514001e-03 ) = -8.109195e-01
g15.Sigma +: f( 2.546751e-03 ) = 2.553184e+01
g15.Sigma +: f( 2.030376e-03 ) = 7.845440e+00
g15.Sigma +: f( 2.030376e-03 ) = 7.845440e+00
g15.Sigma +: f( 2.030376e-03 ) = 7.845440e+00
g15.Sigma +: f( 1.772189e-03 ) = 2.467092e+00
g15.Sigma +: f( 1.772189e-03 ) = 2.467092e+00
g15.Sigma +: f( 1.772189e-03 ) = 2.467092e+00
g15.Sigma +: f( 1.643095e-03 ) = 5.832027e-01
g15.Sigma +: f( 1.643095e-03 ) = 5.832027e-01
g15.Sigma +: f( 1.643095e-03 ) = 5.832027e-01
g15.Sigma +: f( 1.578548e-03 ) = -1.721177e-01
g15.Sigma +: f( 1.578548e-03 ) = -1.721177e-01
g15.Sigma +: f( 1.578548e-03 ) = -1.721177e-01
g15.Sigma +: f( 1.610822e-03 ) = 1.906909e-01
g15.Sigma +: f( 1.610822e-03 ) = 1.906909e-01
g15.Sigma +: f( 1.610822e-03 ) = 1.906909e-01
g15.Sigma +: f( 1.594685e-03 ) = 5.746982e-03
g15.Sigma upper bound:  0.000597058
Datasets              = 1, 2
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2datavar
confidence 1.64-sigma (89.8995%) bounds:
  Param            Best-Fit  Lower Bound  Upper Bound
  -----            --------  -----------  -----------
  g15.Sigma     0.000997626 -0.000907834  0.000597058
method = levmar, stat = chi2datavar, in 17.921696 secs

How does the fit surface vary for a parameter (interval-projection)?

In order to visually inspect the model parameter space we can "project" the statistics into the 1-D or 2-D plane with int_proj and reg_proj, respectively. This allows for checking the shape of the parameter space around the best fit parameter values, as well as evaluate correlations between parameters.

Here we use the int_proj (interval-projecton) method to see how the fit statistic varies with the gamma parameter of the power-law component. Since we already know that the 90% errors of p1.gamma are approximately ±0.13, we choose to set the axis range manually:

sherpa> print get_int_proj()
x     = None
y     = None
min   = None
max   = None
nloop = 20
delv  = None
fac   = 1
log   = False

sherpa> int_proj(p1.gamma,min=1,max=2)

sherpa> print get_int_proj()
x     = [ 1.      1.0526  1.1053 ...,  1.8947  1.9474  2.    ]
y     = [ 161.8675  153.806   146.7844 ...,  144.311   149.5398  155.2686]
min   = 1
max   = 2
nloop = 20
delv  = None
fac   = 1
log   = False

The resulting plot is shown in Figure 3.

Figure 3: Plot of interval-projection results

Figure 3: Plot of interval-projection results

Confidence interval plot of the fit statistic versus the p1.gamma model parameter value.

The "confidence intervals" table above lists a range of common confidence levels and the corresponding change in chi-square values (i.e., the statistic value on the Y-axis in this plot). The parameters displayed by print get_int_proj() show how the fit statistic versus single model parameter plot is calculated. The parameters min and max are the minimum and maximum grid boundary values; if set to the default values of None, the grid boundaries are calculated automatically from the covariance. nloop is the bin size, which by default is used with the min and max grid boundaries in order to determine the step size, delv> (default is delv=None).

The int_unc (interval-uncertainty) command behaves similarly to int_proj , although the fields in the state object for the two methods are different.

How are two parameters correlated (region-projection)?

In this section we use the reg_proj (region-projection) method of Sherpa to see whether the p1.gamma and abs1.nh parameters are correlated.

From our earlier run we know that the 90% errors on the two parameters - when evaluated independently - are approximately 0.13 (gamma) and 0.2 (nH). However we decide to let the routine calculate plot limits automatically, and choose to display contours at the 1 and 1.6 sigma level (68.3% and 90% confidence levels).

sherpa> print get_reg_proj()
x0      = None
x1      = None
y       = None
min     = None
max     = None
nloop   = (10, 10)
fac     = 4
delv    = None
log     = (False, False)
sigma   = (1, 2, 3)
parval0 = None
parval1 = None
levels  = None

sherpa> reg_proj(p1.gamma,abs1.nH,sigma=[1,1.6])

sherpa> print get_reg_proj()
x0      = [ 1.1516  1.2269  1.3021 ...,  1.6786  1.7539  1.8292]
x1      = [ 1.8675  1.8675  1.8675 ...,  2.8652  2.8652  2.8652]
y       = [ 143.7669  145.3427  155.4813 ...,  141.6272  137.5717  138.9838]
min     = [ 1.1516  1.8675]
max     = [ 1.8292  2.8652]
nloop   = (10, 10)
fac     = 4
delv    = None
log     = [False False]
sigma   = [1, 1.6000000000000001]
parval0 = 1.49036863322
parval1 = 2.36635832376
levels  = [ 126.3275  128.4536]

The resulting plot is shown in Figure 4.

Figure 4: Plot of region-projection results

Figure 4: Plot of region-projection results

This plot shows the results of reg_proj() (confidence contour of fit statistic vs. two thawed parameter values), using the default values: 10 points were used along each axis and the range was calcualted automatically. The two contours are drawn at the 68.3% and 90% confidence levels. The plot needs to be re-evaluated using more points, and with a better choice of the axes.

The automatically-chosen limits have resulted in a poor-quality plot: there are not enough data points close to the best-fit location. The easiest way to improve on this is to change and re-run the function, increasing the number of points. We also elect to use a smaller parameter range along both axes to reduce the amount of wasted computation. In a complex case with a larger grid, it may be worthwhile to manually set the limits before running reg_proj, since it may take longer to create a plot.

sherpa> reg_proj(p1.gamma,abs1.nH,min=[1.2,2],max=[1.9,2.8],nloop=[21,21],sigma=[1,1.6]) 

sherpa> print get_reg_proj()
x0      = [ 1.2    1.235  1.27  ...,  1.83   1.865  1.9  ]
x1      = [ 2.   2.   2.  ...,  2.8  2.8  2.8]
y       = [ 136.67    134.4571  133.9499 ...,  138.6876  142.5511  147.6356]
min     = [1.2, 2]
max     = [1.8999999999999999, 2.7999999999999998]
nloop   = [21, 21]
fac     = 4
delv    = None
log     = [False False]
sigma   = [1, 1.6000000000000001]
parval0 = 1.49036863322
parval1 = 2.36635832376
levels  = [ 126.3275  128.4536]

The resulting plot is shown in Figure 5.

Figure 5: Improved region-projection results

Figure 5: Improved region-projection results

The results are improved by using more points along each axis and restricting the ranges of the two parameters used for the contour plot. However the contours still do not appear smooth.

Although the results are much better the contours still do not appear smooth. We now try changing the nloop value, which affects the calculated step-size to see whether this will improve the appearance of the plot:

sherpa> reg_proj(p1.gamma,abs1.nH,min=[1.2,2],max=[1.9,2.8],nloop=[100,100],sigma=[1,1.6])

The resulting plot is shown in Figure 6, which is a smooth contour plot.

Figure 6: Improved region-projection results (nloop=[100,100])

Figure 6: Improved region-projection results (nloop=[100,100])

By increasing the nloop field to 100 we have been able to create a sensible-looking contour plot.

When we utilize print get_reg_proj(), we are given information on the most recent confidence contour plot produced by reg_proj. The min and max lists of grid boundaries are calculated automatically from the covariance when they are set to the default "None". However, they may be set manually as min=[x_min,y_min] and max=[x_max,y_max]. nloop is a list of bin sizes of the x- and y-axes, which by default are used with the min and max grid boundaries to determine the x- and y-axes step sizes, given as the list, delv. The parameter of most interest for reg_proj is sigma, which is a list of the number of sigma at which to plot contours. The levels parameter is subsequently determined after executing reg_proj, which are the confidence level z-values for each sigma.

NOTE for int_proj and reg_proj: In the current version of Sherpa, the log parameter should be left at its default value of False, as the tools do not properly scale plots with logarithmic spacing.

The reg_unc (region-uncertainty) command behaves similarly, although the fields in the state object for the two methods are different. The two commands differ in that reg_unc fixes all other thawed parameters to their best-fit values, rather than being allowed to float to new best-fit values as in reg_proj. This makes reg_unc contours less accurate, but quicker to create.

Improving Visualization

While a little tedious, it is possible to create a color-coded plot with the array values determined by reg_proj, in a new window. We open a new window using add_window and assign the calculated arrays to variables. By using print get_reg_proj().x0, print get_reg_proj().x1, print get_reg_proj().y , we see that the x0, x1, and y arrays correspond to the power-law index, hydrogen column density, and uncertainty, respectively. It should be noted that these three data arrays are hidden when get_reg_proj is printed out due to a bug that will be addressed in the future.

sherpa> add_window()
sherpa> rp=get_reg_proj()
sherpa> xp=rp.x0
sherpa> yp=rp.x1
sherpa> zp=rp.y

We now can plot the contours at the confidence levels that have been calculated, using the ChIPS command add_contour and the Sherpa command get_reg_proj().levels.

sherpa> smin=get_reg_proj().levels
sherpa> add_contour(xp,yp,zp,[smin[0]],["color","red"])
sherpa> add_contour(xp,yp,zp,[smin[1]],["color","green"])
sherpa> limits(X_AXIS,1.2,1.9)
sherpa> limits(Y_AXIS,2,2.8)

sherpa> add_line(1.25,2.75,1.32,2.75,["color","red"])
sherpa> add_label(1.325,2.75,"1\sigma",["valign",0.5,"color","red"])
sherpa> add_line(1.25,2.72,1.32,2.72,["color","green"])
sherpa> add_label(1.325,2.72,"1.6\sigma",["valign",0.5,"color","green"])

This produces the color-coded plot shown in Figure 7.

Figure 7: Color Coded Contour Plot

Figure 7: Color Coded Contour Plot

A new confidence level plot using the 1.0 and 1.6 sigma contour levels.

The ChIPS tool pick is used to determine the location on the window to place the legend elements.

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