Calculating Model Flux and Flux Uncertainty

Getting Started

Please follow the Obtaining Data Used in Sherpa thread. In this example, we use the spectral data products for 3C 273, located in the sample_flux sub-directory.

Fitting a Model to the Spectrum

To begin, we load the spectrum and instrument response files for quasar 3C 273 using the load_pha command, and then fit it with a source model composed of the sum of an absorbed, non-thermal power-law model, and a thermal bremsstrahlung model component. We subtract off the background before fitting, and include only the 0.3–7.0 keV data range in our analysis.

sherpa> load_pha("3c273.pi")
WARNING: systematic errors were not found in file '3c273.pi'
statistical errors were found in file '3c273.pi' 
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
WARNING: systematic errors were not found in file '3c273_bg.pi'
statistical errors were found in file '3c273_bg.pi' 
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi

sherpa> notice(0.3, 7.0)
dataset 1: 0.00146:14.9504 -> 0.2482:9.8696 Energy (keV)

sherpa> subtract()
sherpa> set_stat("chi2datavar")
sherpa> set_method("simplex")
sherpa> set_source(xsphabs.gal * (xspowerlaw.p1 + xsbremss.kt))
sherpa> gal.nh = 0.07
sherpa> kt.kt = 0.5
sherpa> freeze(gal)
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2datavar
Initial fit statistic = 1.31212e+11
Final fit statistic   = 36.4503 at function evaluation 1435
Data points           = 44
Degrees of freedom    = 40
Probability [Q-value] = 0.630838
Reduced statistic     = 0.911256
Change in statistic   = 1.31212e+11
   p1.PhoIndex    2.09439     
   p1.norm        0.000213361 
   kt.kT          0.0526971   
   kt.norm        0.141333    

For reference, the best-fit looks like Figure 1:

sherpa> plot_fit(xlog=True, ylog=True)

Figure 1: Best-Fit Model

Figure 1: Best-Fit Model

Having performed an initial fit to our data set—choosing \(\chi^{2}\) statistics with variance calculated from the data, and the robust Neldermead-Simplex fit optimization—we are ready to run the sample_flux function to obtain an estimate of the energy flux for our source 3C 273, with the corresponding model parameters and associated flux uncertainty.

Sampling the Model Flux

The sample_flux function takes a number of arguments, the key ones being the model component(s) for which to calculate the associated flux and flux uncertainty, the number of flux simulations to run, and the range of data to consider in the analysis. The returned values include the flux from the full model expression used in the fitting (labeled "original model flux"), and the emission specific to the entered model sub-component(s) (labeled "model component flux").

In this thread example we consider the unabsorbed flux due to the sum of the power-law and bremsstrahlung model components (p1+kt), in the 0.3–7 keV energy range, using only ten simulations (in standard analysis, a larger number, e.g. 100 or more, should be used for a better sampling of the parameter distributions). We also set the Boolean 'correlated' parameter to 'True' to include correlations between the parameters in evaluating the uncertainty.

sherpa> component = p1 + kt
sherpa> print(component)
(xspowerlaw.p1 + xsbremss.kt)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   p1.PhoIndex  thawed      2.09439           -3           10           
   p1.norm      thawed  0.000213361            0        1e+24           
   kt.kT        thawed    0.0526971       0.0001          200        keV
   kt.norm      thawed     0.141333            0        1e+24           

sherpa> sample_flux(component, 0.3, 7 , num=10, correlated=True)
WARNING: hard minimum hit for parameter kt.norm
original model flux = 8.74346e-13, + 5.26023e-14, - 2.69053e-14
model component flux = 1.20476e-12, + 1.06352e-13, - 7.62144e-14

[array([8.74345524e-13, 9.26947862e-13, 8.47440234e-13]),
 array([1.20476151e-12, 1.31111323e-12, 1.12854710e-12]),
 array([[9.30448273e-13, 2.11670536e+00, 2.25943504e-04, 5.32346046e-02,
         1.50612411e-01, 0.00000000e+00, 3.81991383e+01],
        [7.97334299e-13, 2.04243058e+00, 1.99564006e-04, 6.84301063e-02,
         0.00000000e+00, 1.00000000e+00, 5.14706557e+01],
        [7.72562890e-13, 2.02968980e+00, 1.88165962e-04, 3.92915321e-02,
         3.53227871e-01, 0.00000000e+00, 4.67410821e+01],
        [8.62077222e-13, 2.05298429e+00, 2.04528880e-04, 4.75095136e-02,
         2.58904266e-01, 0.00000000e+00, 3.73707110e+01],
        [9.52655209e-13, 2.04474053e+00, 2.17037889e-04, 6.10347174e-02,
         8.75530216e-02, 0.00000000e+00, 3.93423060e+01],
        [8.42097587e-13, 2.13959034e+00, 2.15628651e-04, 4.29280632e-02,
         2.76485506e-01, 0.00000000e+00, 3.84064136e+01],
        [9.17946805e-13, 2.03196941e+00, 2.22190102e-04, 4.38093027e-02,
         2.37653807e-01, 0.00000000e+00, 4.34527083e+01],
        [8.61178469e-13, 2.09183946e+00, 2.17948627e-04, 3.95636494e-02,
         3.06846156e-01, 0.00000000e+00, 4.05559455e+01],
        [8.74345524e-13, 2.17413439e+00, 2.33372359e-04, 3.42397260e-02,
         3.72128811e-01, 0.00000000e+00, 4.54029629e+01],
        [8.16672348e-13, 2.06041404e+00, 2.06500568e-04, 6.93343760e-02,
         0.00000000e+00, 1.00000000e+00, 4.97757245e+01],
        [9.02845303e-13, 2.09696727e+00, 2.20840314e-04, 4.84862508e-02,
         2.10116128e-01, 0.00000000e+00, 3.71136706e+01]])]

The screen output shows the 0.3–7 keV energy flux with upper and lower bounds in units of [ergs/s/cm²], for the full model fit to the data, followed by the flux for the specified model component, p1+kt. The flux is obtained as a median of all the flux values in the simulated parameter sets. The default upper and lower bounds are set at 68%, or 1σ, confidence, and may be changed via the 'confidence' parameter.

By default, sample_flux assumes the units of a standard X-ray data set (the 'Xrays' parameter is set to 'True'), returning the energy flux in [ergs/s/cm²]. If 'Xrays' is set to 'False' instead, then the units are determined by the type of input data. For example, for an optical spectrum entered in the units of [W/m²/Hz], defined in frequency (Hz), sample_flux will integrate the spectrum and give the flux in [W/m²].

The sample_flux function also returns a list of arrays. The first two arrays contain the flux with upper and lower bounds for the original model and the specified model component(s), respectively. The last array contains all the simulated parameters with a corresponding flux for the full model expression fit to the data.

To consider another example of the usage of the sample_flux function, we estimate the flux of an individual model component from our full model expression, the thermal bremsstrahlung model component 'kt', within the 0.5–2 keV range at 90% confidence. This time, we choose to use the default 'False' setting for the 'correlated' parameter, increase the number of simulations to 1000, and store the returned values in a variable called sample1.

sherpa> sample1 = sample_flux(kt, 0.5, 2, num=1000, confidence=90)
WARNING: hard minimum hit for parameter kt.norm
WARNING: Covariance failed for 'kt.norm', trying Confidence...
kt.norm lower bound:    -0.10618
kt.norm upper bound:    5.89866
original model flux = 3.906e-13, + 5.81815e-14, - 4.19351e-14
model component flux = 3.32545e-15, + 9.06334e-14, - 3.32291e-15

sherpa> len(sample1)
        3
sherpa> sample1[0]
        array([3.89814872e-13, 4.60003736e-13, 3.53452651e-13])
sherpa> sample1[1]
        array([4.21469347e-15, 1.13921936e-13, 2.47708320e-18]
sherpa> sample1[2].shape
        (1001, 7)

The screen output from sample_flux shows the flux values and 90% bounds from both the full model and the single model component 'kt', along with messages from the covariance/confidence calculation of model parameter uncertainties.

The sample1 variable can be used to access the original and model component flux values—sample1[0] and sample1[1] respectively—and the flux, parameter values, and statistic for each simulation (sample1[2]). The following returns the fluxes and photon index values:

sherpa> fluxes = sample1[2][:, 0]
sherpa> phoindxs = sample1[2][:, 1]
sherpa> plot_scatter)(phoindxs, fluxes, xlabel=r'$\Gamma$', ylabel='flux')

The result of the scatter plot is Figure 2.

Figure 2: Scatter Plot: Photon-Index versus Flux

Figure 2: Scatter Plot: Photon-Index versus Flux

sherpa> print(get_source())
(xsphabs.gal * (xspowerlaw.p1 + xsbremss.kt))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   gal.nH       frozen         0.07            0        1e+06 10^22 atoms / cm^2
   p1.PhoIndex  thawed      2.09439           -3           10           
   p1.norm      thawed  0.000213361            0        1e+24           
   kt.kT        thawed    0.0526971       0.0001          200        keV
   kt.norm      thawed     0.141333            0        1e+24           

sherpa> sample1[2][0]

array([3.82274270e-13, 2.03273525e+00, 2.09775085e-04, 5.37676331e-02,
       2.66725095e-01, 0.00000000e+00, 5.08529834e+01])

The first element of this array is the flux, and then the thawed parameter values in the same order as listed by the get_source command, with the last element being the statistic value.

The plot_trace function lets us look at trace plots (i.e. a value per iteration). For example, to look at how the fluxes and photon-indexes vary you could say (Figure 3):

sherpa> plt.clf()
sherpa> plt.subplot(2, 1, 1)
sherpa> plot_trace(fluxes, name='flux', clearwindow=False)
sherpa> plt.title('')
sherpa> plt.subplot(2, 1, 2)
sherpa> plot_trace(phoindxs, name=r'$\Gamma$', clearwindow=False)
sherpa> plt.title('')
sherpa> plt.subplots_adjust(hspace=0.05)

Figure 3: Trace Plot: Flux and Photon Index

Figure 3: Trace Plot: Flux and Photon Index

Visualizing the Results of the Flux Simulations

We can visualize the distribution of model parameter values resulting from the sample_flux simulations, for a particular model component, using the Sherpa plot_pdf or plot_cdf functions. The plot_pdf function displays a binned probability density function (Figure 4), while plot_cdf shows the cumulative distribution function with lower, median, and upper quantiles (Figure 5).

Figure 4: Power-Law Photon Index Probability Density

Figure 4: Power-Law Photon Index Probability Density

The probablity density for any of the parameters can be displayed using plot_pdf; in this case the p1.PhoIndex parameter:

sherpa> plot_pdf(sample1[2][:,1], xlabel="p1.PhoIndex", name="p1.PhoIndex")

Figure 5: Power-Law Photon Index Cumulative Distribution

Figure 5: Power-Law Photon Index Cumulative Distribution

The cumulative distribution of the photon index is shown by plot_cdf; the median value is marked by the orange dashed line and 68% quantiles with the blue dashed lines.

sherpa> plot_cdf(sample1[2][:,1], xlabel="p1.PhoIndex", name="p1.PhoIndex")

Supplying the Scales for the Simulations

sample_flux assumes a multi-variate Gaussian distribution for the flux simulations, with the Gaussian scales for each parameter given by the covariance matrix (unless covariance fails, in which case confidence is used to calculate the sizes). Correlations are taken into account when the 'correlated' parameter is set to 'True' and the non-diagonal elements from the covariance matrix set the scales.

There is an option to supply the scales to sample_flux for the simulations via the 'scales' parameter, instead of recalculating them with covariance or confidence. If 'correlated=True', the covariance function must be run before sample_flux in order to set the scales. If parameters are uncorrelated ('correlation=False'), a list of scales corresponding to each parameter is needed.

Considering an example using correlated parameters with 'correlated=True', we run the Sherpa covariance command and obtain the resulting covariance matrix to set the scales for the sample_flux simulations.

sherpa> covar()
WARNING: hard minimum hit for parameter kt.norm
Dataset               = 1
Confidence Method     = covariance
Iterative Fit Method  = None
Fitting Method        = neldermead
Statistic             = chi2datavar
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.PhoIndex       2.09439   -0.0630782    0.0630782
   p1.norm       0.000213361 -1.18751e-05  1.18751e-05
   kt.kT           0.0526971   -0.0137291    0.0137291
   kt.norm          0.141333        -----      0.18303

sherpa> matrix = get_covar_results().extra_output
sherpa> print(matrix)
[[ 3.97885901e-03  5.06249904e-07 -1.91128060e-04  2.04055192e-03]
 [ 5.06249904e-07  1.41018755e-10 -3.87223988e-08  4.19621646e-07]
 [-1.91128060e-04 -3.87223988e-08  1.88488592e-04 -2.44765348e-03]
 [ 2.04055192e-03  4.19621646e-07 -2.44765348e-03  3.35000987e-02]]

The scales for sample_flux may now be defined using the covariance results from get_covar_results, using the command directly, as shown below, or the 'matrix' variable defined above.

sherpa> sample2 = sample_flux(component, 0.5, 7, num=1000, correlated=True, scales=matrix)
original model flux = 7.57634e-13, + 3.21046e-14, - 2.92545e-14
model component flux = 8.55786e-13, + 3.55797e-14, - 3.35509e-14

sherpa> print(sample2[2][0:20, 1])
[2.20709007 2.0606933  2.14378853 2.04113953 2.10293264 2.05139932
 2.1680747  2.21411571 2.11991748 2.0808406  2.10909223 2.1087087
 1.948618   2.082921   2.08365138 2.09032091 2.05502353 2.07327993
 2.11088092 2.07389675]

Figure 6: Comparison of the Power-Law Distributions

Figure 6: Comparison of the Power-Law Distributions

Here we compare the power-law distribution from Figure 4 (blue line) to that obtained when using the covariance matrix (orange line):

sherpa> plot_pdf(sample1[2][:, 1])
sherpa> plot_pdf(sample2[2][:, 1], overplot=True)

Saving the Simulation Results

The arrays of model parameter values and associated fluxes from our two sample_flux runs may be written to file using the Sherpa save_arrays command, as shown below. Here, we write the power-law photon index and thermal bremsstrahlung temperature arrays from the first sample_flux run to the ASCII file sample1.dat, specifying data column names and overwrite the existing sample1.dat file.

sherpa> save_arrays("sample1.dat", [sample1[2][:,1], sample1[2][:,3]], ['Photon Index', 'Temperature'], clobber=True)

Scripting It

The file, fit.py , performs the primary commands used above. It can be executed by typing %fit -i 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, et cetera.)