Calculating Uncertainties by Simulating Flux Distributions

Introduction to the Flux Error Functions

Sherpa includes functions for creating and plotting flux probability distributions of Sherpa models—that is, the probability distribution of flux values that accounts for uncertainties in the model parameter values.

All thawed model parameters are sampled from the Gaussian distribution, where the mean is set as the best-fit parameter value and the variance is determined by the diagonal elements of the covariance matrix. The univariate Gaussian is assumed by default; treating each parameter independently. If there are correlations between parameters then the multivariate Gaussian is available, which utilizes the off-diagonal elements of the covariance matrix. The flux is calculated for each sampled set of the parameters and recorded. The histogram of the simulated flux values gives the flux probability distribution for the assumed model.

The Sherpa flux error functions include the following:

• sample_flux — return a sample of Sherpa model parameters and the corresponding flux and flux uncertainty
• sample_energy_flux — sample from the energy flux distribution
• sample_photon_flux — sample from the photon flux distribution
• get_energy_flux_hist — create histograms of the simulated energy flux values
• get_photon_flux_hist — create histograms of the simulated photon flux values
• plot_energy_flux — plot the energy flux probability distribution histogram
• plot_photon_flux — plot the photon flux probability distribution histogram

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 pha_intro sub-directory.

Fitting a Model to the Spectrum

We begin by loading the source spectrum, background, ARF and RMF with load_pha; filtering the data to model the 0.5–7.0 keV energy range; and subtracting the background.

sherpa> load_pha("3c273.pi")
sherpa> notice(0.5,7.0)
dataset 1: 0.00146:14.9504 -> 0.4672:9.8696 Energy (keV)

sherpa> subtract()

Choosing \(\chi^{2}\) statistics, with variance calculated from the data, we fit an absorbed power-law model to the data using a neutral hydrogen column density of 0.07×10²² atoms/cm².

sherpa> set_stat("chi2datavar")
sherpa> set_source(xsphabs.abs1*xspowerlaw.p1)
sherpa> abs1.nh = 0.07
sherpa> freeze(abs1)
sherpa> print(get_source())
(xsphabs.abs1 * xspowerlaw.pl)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nH      frozen         0.07            0        1e+06 10^22 atoms / cm^2
   p1.PhoIndex  thawed            1           -3           10           
   p1.norm      thawed            1            0        1e+24   

sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 1.2744e+11
Final fit statistic   = 35.9809 at function evaluation 22
Data points           = 42
Degrees of freedom    = 40
Probability [Q-value] = 0.651765
Reduced statistic     = 0.899521
Change in statistic   = 1.2744e+11
   p1.PhoIndex    2.09882      +/- 0.061386    
   p1.norm        0.000214303  +/- 1.15866e-05 

Note that it is possible to do freeze(abs1.nh) instead of freeze(abs1), which is useful when you only want to freeze a subset of the parameters of the model component.

Estimating Uncertainties by Simulating the Flux Distribution

We now simulate the energy flux between 0.5–7.0 keV one-hundred times, using sample_energy_flux, which returns one or more samples of the model flux, accounting for errors in the model parameters. The simulations choose parameters assuming their Gaussian distributions, with the mean set to the best-fit values and the variance given by the covariance matrix for all thawed parameters; the flux is then calculated for these parameters. The results will contain one-hundred realizations of the model and one-hundred values of the energy flux.

sherpa> flux100 = sample_energy_flux(0.5, 7, num=100)

Similarly, we may do this with the photon flux distribution using sample_photon_flux. Examining the structure of the array produced by sample_energy_flux:

sherpa> flux100.shape
        (100, 4)

sherpa> print(flux100[0:5,:])
[[7.42967316e-13 2.06347054e+00 2.05165034e-04 0.00000000e+00]
 [8.23185464e-13 2.03875795e+00 2.23434956e-04 0.00000000e+00]
 [7.44270428e-13 2.13751647e+00 2.15983885e-04 0.00000000e+00]
 [7.24188050e-13 2.14987214e+00 2.11842903e-04 0.00000000e+00]
 [7.11194555e-13 2.14271208e+00 2.07082478e-04 0.00000000e+00]]

where we have used NumPy indexing and slicing to view a subset of the data.

Determining Statistical Values from the Sample of Flux Values

The attributes of the array are the values of flux, photon spectral index, normalization factor, and a flag column indicating if the parameters were clipped (1) or not (0), respectively. The clipping flag was introduced in CIAO 4.13 to indicate whether or not the parameter value falls on the parameter limits. We can operate on these arrays to determine additional statistical values using the built-in NumPy functions, numpy.mean, numpy.median, and numpy.std from the distribution determined by the flux simulation.

Python arrays are called on as [row,column], so we can examine only the flux elements of the flux100 array by doing:

sherpa> print(flux100[:,0])
[7.42967316e-13 8.23185464e-13 7.44270428e-13 7.24188050e-13
 7.11194555e-13 7.78059151e-13 7.62872852e-13 7.42051895e-13
 8.57759373e-13 7.84995423e-13 7.56465425e-13 7.67538618e-13
 …
 7.43596551e-13 7.69685425e-13 7.08308244e-13 7.71557158e-13
 8.53540703e-13 8.74070366e-13 8.17605660e-13 7.68909104e-13
 8.12366418e-13 7.59807314e-13 7.27155694e-13 7.27796573e-13]

and using the NumPy functions to determine the mean, median, and standard deviation of the sample of fluxes; in units of ergs/cm²/sec, we find:

sherpa> fluxes = flux100[:,0]

sherpa> np.mean(fluxes)
        7.692973996980474e-13

sherpa> np.median(fluxes)
        7.594420125103322e-13

sherpa> np.std(fluxes)
        4.7560360108049994e-14

---

Determining Quantile Values of the Flux Sample

By estimating the number of points in the array—which is known from the number of simulations performed using the information from the len and NumPy sort functions—the quantiles for this sample may be obtained.

First, we sort the fluxes into a 1D array in ascending order:

sherpa> sf = np.sort(fluxes)
  
sherpa> print(sf[:10]) # check sorting by printing the first 10 elements
[7.01192991e-13 7.06895554e-13 7.07119996e-13 7.08308244e-13
 7.10231422e-13 7.11194555e-13 7.12932613e-13 7.13031638e-13
 7.13871859e-13 7.15247757e-13]

and then get the value of the flux for the 95% quantile and 50% quantile, which corresponds to the median. Since we ran one-hundred simulations, we know that we are looking for the 95^th and 50^th elements of the array. But more generally, the quantiles can be determined using the length function, len, to determine the number of elements in the column.

sherpa> sf[int(0.95*len(fluxes)-1)]
        8.662378127900178e-13

sherpa> sf[int(0.5*len(fluxes)-1)]
        7.590767107395887e-13

Remembering that the array indicies begin at zero, we subtract one from the quantile value.

---

Probability Density and Cumulative Distributions

The plotting functions plot_pdf and plot_cdf provide a simple way to display the probability density and cumulative distribution for the flux arrays.

We start with generating a large number of simulations and plotting the PDF and CDF functions of the flux arrays stored in the first element of the array, numbered 0. We also use get_cdf_plot to obtain the characteristic values for the distribution.

sherpa> flux10000 = sample_energy_flux(0.5, 7, num=10000)

sherpa> flux10000.shape
        (10000, 4)

sherpa> f = flux10000[:,0]

sherpa> plot_pdf(f, bins=40)

sherpa> plot_cdf(f)

sherpa> cplot = get_cdf_plot()

sherpa> cplot.median
        7.578842996703238e-13

sherpa> cplot.lower
        7.083522098815013e-13

sherpa> cplot.upper
        8.109447930504374e-13

which produces the PDF (Figure 1) and CDF (Figure 2) plots shown below.

Figure 1: Flux Probability Density

Figure 1: Flux Probability Density

Figure 2: Flux Cumulative Distribution

Figure 2: Flux Cumulative Distribution

Creating a Histogram of a Flux Probability Distribution

We can visualize the flux distribution with a histogram created by the get_energy_flux_hist function and check whether the distribution is Gaussian or another shape. We simulate the flux distribution between 0.5–7.0 keV ten-thousand times, to produce good statistics. We use the function get_energy_flux_hist, which produces a data object that contains all the information about the simulated sample of parameters and a histogram, normalized to unity, representing the flux probability distribution. By default, get_energy_flux_hist creates a histogram with 75 bins; however, the optional parameter, bins, may be included to change the binning.

sherpa> hist10000 = get_energy_flux_hist(0.5, 7, num=10000)

Similarly, this may be done in photon units with get_photon_flux_hist. To learn more about the data product produced by get_photon_flux_hist, enter 'ahelp "get_photon_flux_hist"' on the Sherpa command-line.

We may display the flux probability distribution histogram determined by get_energy_flux_hist using plot_energy_flux, which plots the calculated energy flux probability distribution—which is the flux distribution for the model component accounting for the errors on the model parameters (Figure 3).

sherpa> plot_energy_flux(0.5, 7, recalc=False)

Figure 3: Histogram of the Flux Distribution

Figure 3: Histogram of the Flux Distribution

This is a histogram of the energy flux probability distribution, plotted using plot_energy_flux, which is determined by get_energy_flux_hist.

Determining the Uncertainties from the Flux Distributions

Taking the histogram data array, we replot the histogram as a scatter plot to check the probability distribution with a Gaussian function, fitting the Gaussian to the histogram with least-squares statistics.

sherpa> load_arrays(2, hist10000.xlo, hist10000.xhi, hist10000.y, Data1DInt)
sherpa> plot_data(2) 
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with chi2datavar
zeros or negative values found

sherpa> set_source(2, gauss1d.g0)
sherpa> g0.integrate = False
sherpa> guess(2, g0)
sherpa> set_par(g0.ampl, min=0, max=10, val=1)
sherpa> set_stat("leastsq")

sherpa> fit(2)
Dataset               = 2
Method                = levmar
Statistic             = leastsq
Initial fit statistic = 1.00662
Final fit statistic   = 0.0601916 at function evaluation 21
Data points           = 76
Degrees of freedom    = 73
Change in statistic   = 0.946432
   g0.fwhm        1.21002e-13  +/- 3.63177e-14 
   g0.pos         7.53403e-13  +/- 1.88867e-14 
   g0.ampl        0.965293     +/- 0           

We use load_arrays to load the histogram flux probability distribution into dataset ID=2, where hist10000.xlo and hist10000.xhi are the low and high bins in the histogram and hist10000.y is the probability of getting the flux within the bin.

Plotting this fit, we can see the good agreement between the Gaussian distribution and flux probability distribution (Figure 4).

sherpa> plot_fit(2, yerrorbars=False)
sherpa> limits(Y_AXIS, -0.05, 1.05)
sherpa> plt.ylabel("Frequency")
sherpa> plt.xlabel(r"Energy Flux (ergs cm$^{-2}$ sec$^{-1}$)")

Figure 4: Gaussian Fitted Flux Distribution

Figure 4: Gaussian Fitted Flux Distribution

The energy flux probability distribution replotted as a scatter-plot with a Gaussian function fitted to the data (red).

We can now determine the uncertainty in the modeled flux by calculating the standard deviation, instead of the FWHM of the Gaussian, using the relationship: \( \mathrm{FWHM} = \sigma \sqrt{8 \ln{2}} \approx 2.35482 \sigma \).

sherpa> fac0 = np.sqrt(8 * np.log(2))
sherpa> sig = g0.fwhm.val/fac0

sherpa> sig
        5.138478634298466e-14

Since we have generated two distributions of fluxes, we can examine the results of the mean of the flux100 sample and compare it with the mean position of the Gaussian fit to the histogram, g0.pos. We see that the two mean values, 〈flux₁₀₀〉 ≅ 7.693×10^-13 ergs cm^-2 s^-1 and 〈hist₁₀₀₀₀〉 ≅ 7.592×10^-13 ergs cm^-2 s^-1, are approximately the same. The standard deviations, similarly, are equivalent to each other, σ_flux100 ≅ 4.756×10^-14 ergs cm^-2 s^-1 and σ_hist10000 ≅ 5.159×10^-14 ergs cm^-2 s^-1. This is expected for normal distributions, as was simulated; however, if the flux distribution were different (e.g. skewed-, bi-modal-, or κ-distributions), then the calculated mean and σ using NumPy may be very different. Visualization of the distributions using the histogram is a useful sanity check.

We can write the array to an ASCII file using the save_arrays function, including column headers:

sherpa> save_arrays("hist10000.dat", [hist10000.xlo,hist10000.xhi,hist10000.y], fields=["xlo","xhi","y"])

and reload the simulation results into another session with the load_ascii function:

sherpa> load_ascii("hist10000.dat", ncols=3, dstype=Data1DInt)

for future analysis.

Scripting It

The file, fit.py , 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.)