Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 3 Nov 2016

URL: http://cxc.harvard.edu/sherpa/threads/confidence_manual/

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

Sherpa Threads (CIAO 4.9 Sherpa v1)


Overview

Synopsis:

This thread uses the "native" Sherpa interface to estimate errors or confidence levels for parameters in a fit (to data of any dimensionality).

Related Links:

Last Update: 3 Nov 2016 - reviewed for CIAO 4.9; updated outputs and fixed typos.


Contents


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 χ2 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.

[TIP]
Tip

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()
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

sherpa> show_stat()

sherpa> set_stat("chi2xspecvar")
Statistic: Chi2XspecVar
Chi Squared with data variance (XSPEC style).

    The variance in each bin is estimated from the data value in that
    bin. See also `Chi2DataVar`, `Chi2Gehrels`, and `Chi2ModVar`.

    The calculation of the variance is the same as `Chi2DataVar`
    except that if the number of counts in a bin is less than 1
    then the variance for that bin is set to 1.

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 29
Data points           = 133
Degrees of freedom    = 130
Probability [Q-value] = 0.630974
Reduced statistic     = 0.954091
Change in statistic   = 6.99384e+09
   abs1.nH        2.36631     
   p1.gamma       1.49034     
   p1.ampl        0.00190277  

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.

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 Δχ2 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.

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).

    The variance in each bin is estimated from the data value in that
    bin. See also `Chi2DataVar`, `Chi2Gehrels`, and `Chi2ModVar`.

    The calculation of the variance is the same as `Chi2DataVar`
    except that if the number of counts in a bin is less than 1
    then the variance for that bin is set to 1.

    

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

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σ and 1.6σ for the normal distribution. The table below displays the relationship between σ and confidence level for one significant parameter.

Confidence intervals for a normal distribution:

Confidence σ Δχ2 Δlog L
68.3% 1.0 1.00 0.50
90.0% 1.6 2.71 1.36
95.5% 2.0 4.00 2.00
99.0% 2.6 6.63 3.32
99.7% 3.0 9.00 4.50

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

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

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

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.192951
p1.gamma lower bound:	-0.13318
abs1.nH upper bound:	0.207426
p1.gamma upper bound:	0.13824
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.49034     -0.13318      0.13824
   abs1.nH           2.36631    -0.192951     0.207426

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.192951
p1.ampl lower bound:	-0.00034192
p1.gamma lower bound:	-0.13318
abs1.nH upper bound:	0.207426
p1.gamma upper bound:	0.13824
p1.ampl upper bound:	0.000430493
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.36631    -0.192951     0.207426
   p1.gamma          1.49034     -0.13318      0.13824
   p1.ampl        0.00190277  -0.00034192  0.000430493

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 in Sherpa as of 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.

Common WARNING messages returned by confidence methods

WARNING: hard minimum hit for parameter <parameter name>

When the confidence, projection, and covariance methods are used to estimate confidence intervals for thawed model parameters after a fit, sometimes a hard upper or lower limit will be reached for one or more parameters. This produces the message "WARNING: hard minimum hit for parameter <parameter name>", along with a row of dashes in the appropriate place in the function output.

[NOTE]
Note

The covariance method can also return a null value for an upper/lower limit when the parameter space at the minimum is non-quadratic for a given parameter. The covariance matrix calculations assume that the parameters follow the normal distribution. If the parameter space is non-smooth, then the covariance calculations fail and Sherpa returns "-----".

Example confidence output:

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

This occurs when the parameter bound found by one of the confidence methods lies outside the hard limit boundary for a model parameter—this could result from an issue with the signal-to-noise of the data, the applicability of the model to the data, systematic errors in the data, among others things.

A parameter hard limit represents either a hard physical limit (e.g., temperature is not allowed to go below zero), a mathematical limit (e.g., prevent a number from going to zero or below, when the log of that number will be taken), or the limit of what a float or double can hold (the fit should not be driven above or below the maximum or minimum values a variable can hold). For this reason, model parameter hard limits should not be changed by the user.

WARNING: The confidence level lies within <interval>

Another warning message which may be returned by confidence is that a model parameter lies within the stated range:

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

This occurs where confidence cannot locate the root (minimum value of the fit statistic function) even though the root is bracketed within an interval (perhaps due to poor resolution of the data or a discontinuity). In such cases, when the openinterval option of confidence is set to False (default), the confidence function will not be able to find the root within the set tolerance and the function will return the average of the open interval which brackets the root. If the option openinterval is set to True, then confidence will print the minimal open interval which brackets the root (not to be confused with the lower and upper bound of the confidence interval). The most accurate thing to do is to return an open interval where the root is localized/bracketed rather than the average of the open interval (since the average of the interval is not a root within the specified tolerance).

The output from confidence may be checked by setting 'set_conf_opt('verbose',1)', and then re-running confidence for the relevant parameter(s).

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.1579, 1.2105, 1.2632, 1.3158, 1.3684, 1.4211,
  1.4737, 1.5263, 1.5789, 1.6316, 1.6842, 1.7368, 1.7895, 1.8421, 1.8947,
  1.9474, 2.    ]
y     = [ 161.8675, 153.806 , 146.7844, 140.7715, 135.7356, 131.6443, 128.4651,
  126.1651, 124.7112, 124.0705, 124.2098, 125.0964, 126.6977, 128.9813,
  131.9156, 135.469 , 139.6109, 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.

[Plot produced with int_proj]
[Print media version: Plot produced with int_proj]

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 χ2 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 σ 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.1499, 1.2255, 1.3012, 1.3769, 1.4525, 1.5282, 1.6038, 1.6795, 1.7551,
  1.8308, 1.1499, 1.2255, 1.3012, 1.3769, 1.4525, 1.5282, 1.6038, 1.6795,
  1.7551, 1.8308, 1.1499, 1.2255, 1.3012, 1.3769, 1.4525, 1.5282, 1.6038,
  1.6795, 1.7551, 1.8308, 1.1499, 1.2255, 1.3012, 1.3769, 1.4525, 1.5282,
  1.6038, 1.6795, 1.7551, 1.8308, 1.1499, 1.2255, 1.3012, 1.3769, 1.4525,
  1.5282, 1.6038, 1.6795, 1.7551, 1.8308, 1.1499, 1.2255, 1.3012, 1.3769,
  1.4525, 1.5282, 1.6038, 1.6795, 1.7551, 1.8308, 1.1499, 1.2255, 1.3012,
  1.3769, 1.4525, 1.5282, 1.6038, 1.6795, 1.7551, 1.8308, 1.1499, 1.2255,
  1.3012, 1.3769, 1.4525, 1.5282, 1.6038, 1.6795, 1.7551, 1.8308, 1.1499,
  1.2255, 1.3012, 1.3769, 1.4525, 1.5282, 1.6038, 1.6795, 1.7551, 1.8308,
  1.1499, 1.2255, 1.3012, 1.3769, 1.4525, 1.5282, 1.6038, 1.6795, 1.7551,
  1.8308]
x1      = [ 1.865 , 1.865 , 1.865 , 1.865 , 1.865 , 1.865 , 1.865 , 1.865 , 1.865 ,
  1.865 , 1.9764, 1.9764, 1.9764, 1.9764, 1.9764, 1.9764, 1.9764, 1.9764,
  1.9764, 1.9764, 2.0878, 2.0878, 2.0878, 2.0878, 2.0878, 2.0878, 2.0878,
  2.0878, 2.0878, 2.0878, 2.1992, 2.1992, 2.1992, 2.1992, 2.1992, 2.1992,
  2.1992, 2.1992, 2.1992, 2.1992, 2.3106, 2.3106, 2.3106, 2.3106, 2.3106,
  2.3106, 2.3106, 2.3106, 2.3106, 2.3106, 2.422 , 2.422 , 2.422 , 2.422 ,
  2.422 , 2.422 , 2.422 , 2.422 , 2.422 , 2.422 , 2.5334, 2.5334, 2.5334,
  2.5334, 2.5334, 2.5334, 2.5334, 2.5334, 2.5334, 2.5334, 2.6448, 2.6448,
  2.6448, 2.6448, 2.6448, 2.6448, 2.6448, 2.6448, 2.6448, 2.6448, 2.7562,
  2.7562, 2.7562, 2.7562, 2.7562, 2.7562, 2.7562, 2.7562, 2.7562, 2.7562,
  2.8676, 2.8676, 2.8676, 2.8676, 2.8676, 2.8676, 2.8676, 2.8676, 2.8676,
  2.8676]
y       = [ 143.9819, 145.5901, 155.8554, 174.8142, 202.3543, 238.2159, 281.9976,
  333.1706, 391.0964, 455.0492, 142.0187, 135.673 , 137.3918, 147.3268,
  165.5005, 191.8024, 225.9892, 267.6916, 316.4255, 371.6086, 148.513 ,
  135.4185, 129.7709, 131.8043, 141.6432, 159.2959, 184.6502, 217.4751,
  257.4269, 304.0591, 161.5629, 142.7884, 130.848 , 126.0312, 128.537 ,
  138.4646, 155.8074, 180.4503, 212.1713, 250.6474, 179.6386, 156.1178,
  138.8417, 128.1345, 124.247 , 127.346 , 137.5068, 154.7082, 178.8313,
  209.6608, 201.5157, 174.0552, 152.2845, 136.5473, 127.1276, 124.2405,
  128.0231, 138.5285, 155.7219, 179.48  , 226.2187, 195.5095, 169.976 ,
  149.9691, 135.7927, 127.6942, 125.8558, 130.3875, 141.322 , 158.6118,
  252.9727, 219.6036, 190.9387, 167.3278, 149.0841, 136.4752, 129.7146,
  128.9546, 134.2805, 145.7063, 281.1643, 245.6346, 214.3802, 187.7439,
  166.0403, 149.5478, 138.5003, 133.0803, 133.413 , 139.5612, 310.3094,
  273.0414, 239.6605, 210.4988, 185.8667, 166.046 , 151.2826, 141.7797,
  137.6916, 139.1186]
min     = [ 1.14989144  1.86501016]
max     = [ 1.83078663  2.86761199]
nloop   = (10, 10)
fac     = 4
delv    = None
log     = [False False]
sigma   = [1, 1.6]
parval0 = 1.49033903432
parval1 = 2.3663110745
levels  = [ 126.32751446  128.45362719]

The resulting plot is shown in Figure 4.

[Plot produced with reg_proj() command]
[Print media version: Plot produced with reg_proj() command]

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.305, 1.34 , 1.375, 1.41 , 1.445, 1.48 , 1.515,
  1.55 , 1.585, 1.62 , 1.655, 1.69 , 1.725, 1.76 , 1.795, 1.83 , 1.865,
  1.9  , 1.2  , 1.235, 1.27 , 1.305, 1.34 , 1.375, 1.41 , 1.445, 1.48 ,
  ..........
  1.2  , 1.235, 1.27 , 1.305, 1.34 , 1.375, 1.41 , 1.445, 1.48 , 1.515,
  1.55 , 1.585, 1.62 , 1.655, 1.69 , 1.725, 1.76 , 1.795, 1.83 , 1.865,
  1.9  ]
x1      = [ 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  ,
  2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.  , 2.04, 2.04, 2.04,
  2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04, 2.04,
  ..........
  2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76,
  2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 ,
  2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 , 2.8 ]
y       = [ 136.67  , 134.4571, 133.9499, 135.1674, 138.1226, 142.8232, 149.2707,
  157.4606, 167.3825, 179.02  , 192.3508, 207.3466, 223.9736, 242.1926,
  261.9591, 283.2241, 305.9336, 330.0299, 355.4514, 382.1334, 410.0082,
  ..........
  268.0711, 252.0218, 236.8634, 222.6299, 209.3541, 197.0678, 185.8014,
  175.5837, 166.4421, 158.4018, 151.4862, 145.7165, 141.1116, 137.6882,
  135.4602, 134.4389, 134.633 , 136.0483, 138.6876, 142.5511, 147.6356]
min     = [1.2, 2]
max     = [1.9, 2.8]
nloop   = [21, 21]
fac     = 4
delv    = None
log     = [False False]
sigma   = [1, 1.6]
parval0 = 1.49033903432
parval1 = 2.3663110745
levels  = [ 126.32751446  128.45362719]

The resulting plot is shown in Figure 5.

[Plot produced with reg_proj() using more points.]
[Print media version: Plot produced with reg_proj() using more points.]

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.

[Plot produced with reg_proj() using a larger nloop value.]
[Print media version: Plot produced with reg_proj() using a larger nloop value.]

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=[xmin,ymin] and max=[xmax,ymax]. 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 σ; at which to plot contours. The levels parameter is subsequently determined after executing reg_proj, which are the confidence level z-values for each σ.

[NOTE]
Log-space for int_proj and reg_proj

In the current version of Sherpa, the log parameter should be left at its default value of False in int_proj and reg_proj, 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.

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.

[Plot produced with get_reg_proj() and add_contour().]
[Print media version: Plot produced with get_reg_proj() and add_contour().]

Figure 7: Color Coded Contour Plot

A new confidence level plot using the 1.0 and 1.6 σ 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 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.


History

14 Jan 2005 reviewed for CIAO 3.2: no changes
21 Dec 2005 reviewed for CIAO 3.3: no changes
01 Dec 2006 reviewed for CIAO 3.4: no changes
02 Dec 2008 reviewed for CIAO 4.1: updated syntax for CIAO4.1
29 Apr 2009 new script command is available with CIAO 4.1.2
21 Jan 2010 updated for CIAO 4.2: the conf command is available
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
15 Jul 2010 updated to include information about warning messages returned by the confidence method.
03 Sep 2010 figures moved inline with text
30 Jan 2012 reviewed for CIAO 4.4 (no changes)
13 Dec 2012 reviewed for CIAO 4.5 (no changes)
11 Dec 2013 reviewed for CIAO 4.6: updated formatting, content unchanged
18 Mar 2015 reviewed for CIAO 4.7; updated plots and fixed typos, no content change.
10 Dec 2015 reviewed for CIAO 4.8; no content change.
03 Nov 2016 reviewed for CIAO 4.9; updated outputs and fixed typos.


Last modified: 3 Nov 2016
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxchelp@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2017. All rights reserved.