Sherpa Template Models

Sherpa Threads (CIAO 4.5 Sherpa v1)

Overview

Synopsis:

The Sherpa template model is available for comparing a source data spectrum against a set of user-provided template models, in order to find the single template model which best matches the source data. A special grid-searching fit optimization method, gridsearch, is used for fitting a template model to a data set in Sherpa.

Starting with a directory of template model files conforming to a specific format, and a single index file listing the file contents of that directory, the models are loaded into the Sherpa session using the load_template_model function, an extension of load_table_model. After fitting the template model to the source data using the gridsearch fit optimization method, the parameters of the best matched template model are returned. Fit results may be examined in the usual way with the appropriate plotting commands.

The Sherpa template model supports linear, nearest-neighbor, and polynomial interpolation. Interpolation is used by the template model to match the data grid to the model grid - which must match before the fit statistics can be calculated for fitting.

Last Update: 17 Jan 2012 - The Sherpa template model is new in CIAO 4.4.

Getting Started: Setting up the Template Model Files

Model files

A collection of template models may be read into Sherpa from a directory full of template files using the load_template_model function, in order to compare a data set to all of the templates in that collection. The Sherpa model library accommodates ASCII template model files containing separated columns of x and y coordinates, such as wavelength, energy, or time, versus flux, counts, or intensity. In this thread we consider a set of accretion disk models in the Kerr geometry (Siemiginowska et al., ApJ, 454, p.77, 1995), where there is a separate model file for each combination of the black hole mass, accretion rate, and inclination angle parameter values.

The contents of the directory of Kerr model files is shown below, followed by the contents of a selected file, k10_001_05.dat.

% ls /data/Kerr/

k10_001_01.dat   k6_01_05.dat   k7_01_1.dat    k8_03_025.dat   k9_01_075.dat
k10_001_025.dat  k6_01_075.dat  k7_03_01.dat   k8_03_05.dat    k9_01_1.dat
k10_001_05.dat   k6_01_1.dat    k7_03_025.dat  k8_03_075.dat   k9_03_01.dat
k10_001_075.dat  k6_03_01.dat   k7_03_05.dat   k8_03_1.dat     k9_03_025.dat
k10_001_1.dat    k6_03_025.dat  k7_03_075.dat  k8_08_01.dat    k9_03_05.dat
k10_01_01.dat    k6_03_05.dat   k7_03_1.dat    k8_08_025.dat   k9_03_075.dat
k10_01_025.dat   k6_03_075.dat  k7_08_01.dat   k8_08_05.dat    k9_03_1.dat
k10_01_05.dat    k6_03_1.dat    k7_08_025.dat  k8_08_075.dat   k9_08_01.dat
k10_01_075.dat   k6_08_01.dat   k7_08_05.dat   k8_08_1.dat     k9_08_025.dat
k10_01_1.dat     k6_08_025.dat  k7_08_075.dat  k8_1_01.dat     k9_08_05.dat
k10_03_01.dat    k6_08_05.dat   k7_08_1.dat    k8_1_025.dat    k9_08_075.dat
k10_03_025.dat   k6_08_075.dat  k7_1_01.dat    k8_1_05.dat     k9_08_1.dat
k10_03_05.dat    k6_08_1.dat    k7_1_025.dat   k8_1_075.dat    k9_1_01.dat
k10_03_075.dat   k6_1_01.dat    k7_1_05.dat    k8_1_1.dat      k9_1_025.dat
k10_03_1.dat     k6_1_025.dat   k7_1_075.dat   k9_001_01.dat   k9_1_05.dat
k10_08_01.dat    k6_1_05.dat    k7_1_1.dat     k9_001_025.dat  k9_1_075.dat
k10_08_025.dat   k6_1_075.dat   k8_01_01.dat   k9_001_05.dat   k9_1_1.dat
k10_08_05.dat    k6_1_1.dat     k8_01_025.dat  k9_001_075.dat  templ_index.dat
k10_08_075.dat   k7_01_01.dat   k8_01_05.dat   k9_001_1.dat   
k10_08_1.dat     k7_01_025.dat  k8_01_075.dat  k9_01_01.dat
k6_01_01.dat     k7_01_05.dat   k8_01_1.dat    k9_01_025.dat
k6_01_025.dat    k7_01_075.dat  k8_03_01.dat   k9_01_05.dat




% more /data/Kerr/k10_001_05.dat

    13.3830    43.6970
    13.4330    43.8500
    13.4830    43.9320
    13.5330    44.0030
    13.5830    44.0720
    13.6330    44.1400
    13.6830    44.2080
    13.7330    44.2760
    .
    .
    .
    16.4830    33.9740
    16.5330    32.2860
    16.5830    30.3930
    16.6330    28.2710
    16.6830    25.8880
    16.7330    23.2100
    16.7830    20.1980
    16.8330    16.8090
    16.8830    12.9930

Each of the example ASCII-format template files contains columns of frequency and luminosity values defining the predicted spectra emitted by an accretion disk surrounding a supermassive black hole (see section 4.1 of the referenced paper), where local modified blackbody emission is assumed. The first column has units of Log frequency in Hz, and the second is the luminosity in log of the rest frame integrated luminosity Log (nu Lnu) in erg/s.

Index File

The template model index file should contain a table with one line per template data file, with three groups of columns in the following order:

Model parameter columns, an arbitrary number of columns

The MODELFLAG column which separates the parameter list from the filenames/model arrays, and marks lines which are to be used or not: MODELFLAG = 1 - use the file; MODELFLAG = 0 - do not use the file

The FILENAME column which points to the data file for that instance, including the full directory path to the file.

Sherpa reads the index file in order to set up the template model with the parameters specified in the first line, and the arrays from the columns given by the data files.

Our example index file appears below, with model flags (all having value 1) and filenames listed in the two right-most columns, and the black hole mass (Log M_bh/M_sun), accretion rate (Eddington), and inclination angle (cosine theta) model parameters listed from the left. The combination of the three parameter values sets the spectral shape of the model.

% more templ_index.dat

# mass rate angle modelflag filename 

6.0 0.1 0.1    1   /data/Kerr/k6_01_01.dat 
6.0 0.1 0.25   1   /data/Kerr/k6_01_025.dat 
6.0 0.1 0.5    1   /data/Kerr/k6_01_05.dat 
6.0 0.1 0.75   1   /data/Kerr/k6_01_075.dat 
6.0 0.1 1.0    1   /data/Kerr/k6_01_1.dat 

6.0 0.3 0.1    1   /data/Kerr/k6_03_01.dat 
6.0 0.3 0.25   1   /data/Kerr/k6_03_025.dat 
6.0 0.3 0.5    1   /data/Kerr/k6_03_05.dat 
6.0 0.3 0.75   1   /data/Kerr/k6_03_075.dat 
6.0 0.3 1.0    1   /data/Kerr/k6_03_1.dat 

.
.
.

10.0 0.3 0.1    1   /data/Kerr/k10_03_01.dat 
10.0 0.3 0.25   1   /data/Kerr/k10_03_025.dat 
10.0 0.3 0.5    1   /data/Kerr/k10_03_05.dat 
10.0 0.3 0.75   1   /data/Kerr/k10_03_075.dat 
10.0 0.3 1.0    1   /data/Kerr/k10_03_1.dat 

10.0 0.8 0.1    1   /data/Kerr/k10_08_01.dat 
10.0 0.8 0.25   1   /data/Kerr/k10_08_025.dat 
10.0 0.8 0.5    1   /data/Kerr/k10_08_05.dat 
10.0 0.8 0.75   1   /data/Kerr/k10_08_075.dat 
1.0  0.8 1.0    1   /data/Kerr/k10_08_1.dat

Loading the Source Data

After starting Sherpa, we read in our source data spectrum to be fitted by the template model in the usual way, using the appropriate load command. In this example, we use the load_ascii command to read the first three columns of ASCII data from file "source.dat", namely the x (frequency), y (integrated luminsity), and y error columns.

% ciao

% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.4 Sherpa version 1 Friday, December 2, 2011


sherpa> load_ascii('source.dat', ncols=3, dstype=Data1D)

The show_data and get_indep/get_dep commands may be used to check that the data was properly read:

sherpa> show_data()
Data Set: 1
Filter: 17.6983-7.4149 x
name      = source.dat
x         = Float64[46]
y         = Float64[46]
staterror = Float64[46]
syserror  = None


sherpa> get_indep()

(array([ 17.69831459,  15.22778313,  15.14764634,  15.04849851,
        14.69471047,  14.57931705,  14.4512198 ,  14.2154123 ,
        10.4876855 ,  10.34388301,  10.01346923,  10.01346923,
        10.00024097,  10.00024097,   9.74586299,   9.74586299,
         9.46062726,   9.46062726,   9.46062726,   9.46062726,
         9.46062726,   9.46062726,   9.18956049,   9.18956049,
         9.18956049,   8.92515939,   8.87679209,   8.56491923,
         8.56491923,   8.56491923,   8.51861921,   8.49634282,
         8.49634282,   8.24899768,   8.21758921,   8.09265048,
         7.89428282,   7.89428282,   7.73445498,   7.73445498,
         7.71243924,   7.6608522 ,   7.61552922,   7.5372157 ,
         7.48181656,   7.41486977]),)


sherpa> get_dep()

array([ 44.83025145,  46.20275955,  46.03082425,  46.04305871,
        45.63597532,  45.47709506,  45.46993128,  45.41821408,
        43.73912372,  43.74154589,  44.06194405,  44.0665398 ,
        43.96097247,  43.95510354,  43.98644329,  43.9809108 ,
        44.02583231,  44.03575947,  44.05023729,  44.05023729,
        44.08004227,  44.04785759,  44.07782081,  44.10165342,
        44.10165342,  44.14189464,  44.16799744,  44.18745992,
        44.22553677,  44.15139818,  44.18572621,  44.16069821,
        44.16436316,  44.24415224,  44.19602617,  44.18918674,
        44.05729915,  44.13020888,  43.94007455,  43.94007455,
        44.07334351,  43.86696745,  43.87414603,  44.02331465,
        44.05495793,  43.96676315])

We may also visualize the data before fitting with plot_data, and utilize several ChIPS commands to customize the data plot.

sherpa> plot_data()

sherpa> linear_scale()
sherpa> set_plot_xlabel(r"log \nu [Hz]")
sherpa> set_plot_ylabel(r"log \nu L_{\nu} [erg/s]")

Figure 1: Plot of source data to be fitted with the template model

[Print media version: Fit residuals plot]

Figure 1: Plot of source data to be fitted with the template model

Loading the Template Models

We may now use the load_template_model function to read in our collection of templates and assign them to a Sherpa model instance, as demonstrated below.

sherpa> load_template_model("kerr_templ", "templ_index.txt")

sherpa> set_model(kerr_templ)

Here, we have loaded the Kerr disk model templates listed in the index file "templ_index.txt" into the Sherpa model instance "kerr_templ", and assigned this model to our source data set with set_model. We can examine the model parameters using show_model, where the parameter names contained in our index file are shown, along with the parameter values matching the first template in our list.

sherpa> show_model()

Model: 1
templatemodel.kerr_templ
   Param               Type          Value          Min          Max      Units
   -----               ----          -----          ---          ---      -----
   kerr_templ.mass     thawed            6            6           10
   kerr_templ.rate     thawed          0.1         0.01            1
   kerr_templ.angle    thawed          0.1          0.1            1

Using the Gridsearch Optimization Method

In order to fit a Sherpa template model to a data set, 'gridsearch' must be chosen as the fit optimization method. We do this using set_method, after checking the current method setting with show_method:

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> set_method("gridsearch")
sherpa> set_method_opt('sequence', kerr_templ.parvals)


sherpa> show_method()

Optimization Method: GridSearch
name     = gridsearch
num      = 16
sequence = Float64[315]
numcores = 1
maxfev   = None
ftol     = 1.19209289551e-07
method   = None
verbose  = 0

The gridsearch method evaluates the fit statistic for each point in the parameter space grid provided by the user, which is specified as a list in the 'sequence' parameter (the list which the method should iterate through to evaluate the model function). In this example, we have set the 'sequence' parameter to use the full list of mass, rate, and angle parameter values associated with the template model files in our collection. The list of parameter values may be accessed as follows:

sherpa> print kerr_templ.parvals

[[  6.     0.1    0.1 ]
 [  6.     0.1    0.25]
 [  6.     0.1    0.5 ]
 [  6.     0.1    0.75]
 [  6.     0.1    1.  ]
 [  6.     0.3    0.1 ]
 [  6.     0.3    0.25]
 [  6.     0.3    0.5 ]
 [  6.     0.3    0.75]
 [  6.     0.3    1.  ]

    .       .      .   
    .       .      .
    .       .      .

 [ 10.     0.3    0.1 ]
 [ 10.     0.3    0.25]
 [ 10.     0.3    0.5 ]
 [ 10.     0.3    0.75]
 [ 10.     0.3    1.  ]
 [ 10.     0.8    0.1 ]
 [ 10.     0.8    0.25]
 [ 10.     0.8    0.5 ]
 [ 10.     0.8    0.75]
 [ 10.     0.8    1.  ]]

The best match to the source spectrum in this sequence is the grid point with the lowest value of the fit statistic; gridsearch reports back the parameter values associated with this point. (If the 'sequence' parameter is kept at the default value of 'None', the gridsearch 'num' parameter is used to generate a sequence of parameters to evaluate. A uniform grid of size npar^num elements is generated.)

Fitting the Template Model to the Spectrum

Once the source data has been loaded, the template model properly set, and the fit optimization method switched to 'gridsearch' and customized, we may fit the template model to the source spectrum in the usual way, using the fit() command.

In this example, we switch the fit statistic from the default 'chi2gehrels' to 'chi2datavar' before running the fit (keeping chi-squared statistics in order to include the supplied data errors in the fit), and filter the source data so that only the 13.5-17 Hz range is fit.

sherpa> show_stat()

Statistic: Chi2Gehrels
Chi Squared with Gehrels variance

sherpa> set_stat("chi2datavar")

sherpa> notice(13.,17.)

sherpa> fit()
Dataset               = 1
Method                = gridsearch
Statistic             = chi2datavar
Initial fit statistic = 97.4647
Final fit statistic   = 97.4647 at function evaluation 106
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 3.40753e-20
Reduced statistic     = 24.3662
Change in statistic   = 0
   kerr_templ.mass   9
   kerr_templ.rate   0.3
   kerr_templ.angle   1
WARNING: parameter value kerr_templ.angle is at its maximum boundary 1.0

The initial fit results indicate that, of all the template models in our collection, the one contained in template model file k9_03_1.dat represents the best match to the source data.

We can further examine the fit results using one of the plotting commands designed for this puropse, such as plot_fit_delchi or plot_fit_resid, to visualize the quality of the fit.

sherpa> plot_fit_resid()

sherpa> set_plot_xlabel(r"log \nu [Hz]")
sherpa> set_current_plot("plot1")
sherpa> set_plot_ylabel(r"log \nu L_{\nu} [erg/s]")
sherpa> set_current_plot("plot2")
sherpa> set_plot_ylabel("Residuals")

Figure 2: Plot of fit residuals resulting from template model fit

Note that the confidence and projection commands available for estimating errors on model parameters, such as confidence and region projection, are not compatible with the gridsearch optimization method used for fitting template models to source data.

History

17 Jan 2012

The Sherpa template model is new in CIAO 4.4.