Step-by-Step Guide to Estimating Errors and Confidence Levels
Sherpa Threads (CIAO 4.10 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:
- The ahelp files for the Sherpa routines: covariance, conf, int_unc, int_proj, reg_unc, and reg_proj.
Last Update: 1 Jun 2018 - reviewed for CIAO 4.10; no content change.
Contents
- Getting Started
- Confidence limits for individual parameters
- How does the fit surface vary for a parameter (interval-projection)?
- How are two parameters correlated (region-projection)?
- Improving Visualization
- Scripting It
- History
- Images
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.
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.
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 δχ 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
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 | \(\sigma\) | \(\Delta \chi^{2}\) | \(\Delta \log{\mathcal{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 parameter. 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.
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 \(\pm 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.
Figure 3: Plot of interval-projection results
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.
Figure 4: Plot of region-projection results
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.
Figure 5: Improved region-projection results
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])
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 σ; at which to plot contours. The levels parameter is subsequently determined after executing reg_proj, which are the confidence level z-values for each σ.
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,r"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,r"1.6\sigma",["valign",0.5,"color","green"])
This produces the color-coded plot shown in Figure 7.
Figure 7: Color Coded Contour Plot
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.10 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. |
01 Jun 2018 | reviewed for CIAO 4.10; no content change. |