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 3C273, 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 3C273 using the Sherpa 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")
sherpa> notice(0.3,7.0)
sherpa> subtract()

sherpa> set_stat("chi2datavar")
sherpa> set_method("simplex")
sherpa> set_model(xsphabs.gal*(xspowerlaw.p1+xsbremss.kt))
sherpa> gal.nh=0.07
sherpa> set_par(kt.kt,0.5)
sherpa> freeze(gal.nh)
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2datavar
Initial fit statistic = 1.31027e+11
Final fit statistic   = 36.7109 at function evaluation 1062
Data points           = 44
Degrees of freedom    = 40
Probability [Q-value] = 0.619124
Reduced statistic     = 0.917772
Change in statistic   = 1.31027e+11
   p1.PhoIndex    2.10157
   p1.norm        0.000215161
   kt.kT          0.0530235
   kt.norm        0.147625

Having performed an initial fit to our data set - choosing χ² 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 therange 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 subcomponent(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> sample_flux(component, 0.3, 7 , num=10, correlated=True)
WARNING: hard minimum hit for parameter kt.norm
original model flux = 8.72859e-13, + 1.88539e-14, - 3.48363e-14
model component flux = 1.25365e-12, + 5.35084e-14, - 1.84121e-13
           [array([  8.72858729e-13,   8.91712669e-13,   8.38022471e-13]),
 array([  1.25365405e-12,   1.30716242e-12,   1.06953262e-12]),
 array([[  8.54885244e-13,   2.10849665e+00,   2.17146898e-04,
          4.28110030e-02,   2.90078826e-01],
       [  8.72858729e-13,   2.09006734e+00,   2.18353875e-04,
          4.33819275e-02,   3.14472625e-01],
       [  9.07300410e-13,   2.05072001e+00,   2.14032920e-04,
          4.78936823e-02,   3.15353473e-01],
       [  8.38022471e-13,   2.09316882e+00,   2.16832225e-04,
          3.10165642e-02,   3.94096754e-01],
       [  8.91712669e-13,   2.10297168e+00,   2.14691410e-04,
          5.32346550e-02,   1.67086120e-01],
       [  8.24219037e-13,   2.16342389e+00,   2.21704616e-04,
          6.73794277e-02,   0.00000000e+00],
       [  8.90767760e-13,   2.14669854e+00,   2.23486640e-04,
          4.88731833e-02,   2.33116485e-01],
       [  8.71740014e-13,   2.07213101e+00,   2.08726255e-04,
          4.92086607e-02,   2.40499458e-01],
       [  8.64695643e-13,   2.12024844e+00,   2.26299783e-04,
          3.41467467e-02,   3.90466761e-01],
       [  9.32534970e-13,   2.10032802e+00,   2.18909538e-04,
          5.43316237e-02,   1.96510611e-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-sigma, 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 output arrays in a list 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.110017
kt.norm +: WARNING: The confidence level lies within (1.693518e+00, 1.697021e+00)
kt.norm upper bound:    1.54764
original model flux = 3.91493e-13, + 8.40617e-14, - 4.17343e-14
model component flux = 4.59132e-15, + 1.3732e-13, - 4.58955e-15

sherpa> print sample1
.
.
.
[  3.68900622e-13,   2.16021926e+00,   2.08249480e-04,
          5.17067662e-02,   0.00000000e+00],
       [  3.39145762e-13,   2.10079019e+00,   1.89412404e-04,
          6.58446554e-02,   2.01480322e-02]])]

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 'print sample1' command returns the content of the 'sample1' list we defined, which we can use to access the original and model component flux values, and any of the 1000 arrays of model parameters and associated fluxes from the simulations we ran. For example, to access the model flux for our full model expression, the flux from our specified model components, and the parameter and flux results from the first simulation, we would do:

sherpa> sample1[0]
array([  3.91493137e-13,   4.75554822e-13,   3.49758824e-13])

sherpa> sample1[1]
array([  4.59132135e-15,   1.41911588e-13,   1.76836919e-18])

sherpa> sample1[2]
.
.
.
   [  3.68900622e-13,   2.16021926e+00,   2.08249480e-04,
          5.17067662e-02,   0.00000000e+00],
       [  3.39145762e-13,   2.10079019e+00,   1.89412404e-04,
          6.58446554e-02,   2.01480322e-02]])

The third command above returns the full set of model parameter values for the full model, from the first of one thousand simulations. To obtain the values for one or a subset of the model parameters from a given simulation - e.g., the photon index values from the power-law model of the first simulation run - we would specify the corresponding item in the array 'sample1' list:

sherpa> sample1[2][:,1]

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, while plot_cdf shows the cumulative distribution function with lower, median, and upper quantiles.

sherpa> plot_pdf(sample1[2][:,1])
sheroa> plot_cdf(sample1[2][:,1])

The plot produced by the first command above displays the histogram of power-law photon index values divided by the number of simulated values.

Figure 1: Power-law photon index probability density

Figure 1: Power-law photon index probability density

The second plot shows the cummulative distribution of the photon index with the median and 68% quantiles marked.

Figure 2: Power-law photon index cumulative distribution

Figure 2: Power-law photon index cumulative distribution

Supplying the Scales for the Simulations

sample_flux assumes the multivariate 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.10157   -0.0644017    0.0644017
   p1.norm       0.000215161 -1.22168e-05  1.22168e-05
   kt.kT           0.0530235   -0.0151479    0.0151479
   kt.norm          0.147625        -----     0.209316

sherpa> matrix = get_covar_results().extra_output
sherpa> print matrix
[[  4.14757463e-03   5.41450870e-07  -2.65570070e-04   3.12973286e-03]
 [  5.41450870e-07   1.49249492e-10  -5.37740376e-08   6.39457271e-07]
 [ -2.65570070e-04  -5.37740376e-08   2.29458304e-04  -3.10396400e-03]
 [  3.12973286e-03   6.39457271e-07  -3.10396400e-03   4.38132623e-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(kt+p1,0.5,7,num=1000,correlated=True,scales=get_covar_results().extra_output)
original model flux = 7.59443e-13, + 3.15572e-14, - 3.20532e-14
model component flux = 8.63795e-13, + 3.39857e-14, - 3.67103e-14

sherpa> print sample2[2][:,1]
[ 2.13854308  2.1723029   2.01691641  2.16713482  2.14343459  2.06153498
  2.13413177  2.09300187  2.14492589  2.23027366  2.04289602  2.1392189
  2.15751496  2.06901516  2.15352481  2.18385192  2.07055837  2.17068898
  2.05333147  2.05551043  1.99095613  2.18172993  2.08646942  2.07696355
  2.11088105  2.0353201   2.05290183  2.08322651  2.10553819  2.14458421
  2.12975472  2.10372147  2.09307727  2.16081589  2.13037643  2.10208767

...
  
  2.1475424   2.11513573  2.17544644  2.02685666  2.16201298  2.07505881
  2.0779257   2.15500343  2.11544695  2.13293005  2.05692154  2.05882872
  2.02840333  2.02961024  2.15485285  2.28450665  2.16241881  2.10912928
  2.14594701  2.13628782  2.03780289  1.99432489  2.13941771  2.06684829
  2.16819322  2.11703822  2.1442987   2.10240431]

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 'execfile(flux_script.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.