Simple Templates Fitting of a Quasar SED in Sherpa

Sherpa can use the template models and combined them with the other models. Here we show a simple template fitting to the SED of a quasar. A set of accretion disk spectral models with the standard parameters (mass, accretion rate, inclination anlge) has been stored in the subdirectory Templates. table.txt ascii file index these spectra as required for Sherpa.

First import Sherpa packages

In [1]:
from sherpa.astro.ui import *
WARNING: imaging routines will not be available, 
failed to import sherpa.image.ds9_backend due to 
'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'

Define the optimization method and load the ascii data file ( $\log {\nu}, \log {\nu F_{\nu}}$, $1\sigma$ errors), plot the data and set the data filter. The SED covers broad-band and we filter the data to include only the optical-UV part for fitting with the disk models.

In [2]:
set_method('simplex')
load_ascii('3098_errors.dat', ncols=3, dstype=Data1D)
plot_data()
notice(13.5,16.)

Input template model defined via an ascii index file table.txt. The model name for Sherpa session is required and given as 'tbl' in this initial step. set_model() assignes the template model to the SED data. get_model() brings the information about the initial model parameters, which are based on the first entry in the table.txt index file.

In [3]:
load_template_model('tbl', 'table.txt')
set_model(tbl)
print get_model()
template.tbl
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   tbl.mass     thawed            6            6           10           
   tbl.rate     thawed          0.1         0.01            1           
   tbl.angle    thawed          0.1          0.1            1           

In the next step the SED is fit with the templates using simplex neldermead algorithm. The chi2 statistics will be used with the measurement errors entered together with the data. The fit returns the screen output with the information about the best fit parameters and related statistical values.

In [4]:
fit()
plot_fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2
Initial fit statistic = 426180
Final fit statistic   = 91.4691 at function evaluation 311
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 6.4175e-19
Reduced statistic     = 22.8673
Change in statistic   = 426089
   tbl.mass       8.9899      
   tbl.rate       0.306754    
   tbl.angle      0.995986    

plot_fit() generates a figure showing the data points overplotted with the line defined by the best fit mode parameters given above. The fit goes over the data points, but the reduced chi2 statistics is large, mainly due to the deviations at lower frequencies. We add a second model component to account for these deviations.

In [5]:
ignore()
notice(14.5,15.4)
fit()
plot_fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2
Initial fit statistic = 26.2259
Final fit statistic   = 19.3396 at function evaluation 258
Data points           = 5
Degrees of freedom    = 2
Probability [Q-value] = 6.31633e-05
Reduced statistic     = 9.66979
Change in statistic   = 6.88635
   tbl.mass       9.00932     
   tbl.rate       0.294398    
   tbl.angle      0.974253    

The reduced chi2 is still high. In order to obtain the confidence bounds for the best fit parameters we will set up the options for confidence before running it. 'sigma' defines the confidence level, 'max_rstat' defines the maximum value of reduced statistics allowed by confidence, if rstat is higher than the max_rstat conf() will fail.

In [6]:
set_conf_opt('sigma',1)
set_conf_opt('max_rstat',10)
conf()
Dataset               = 1
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = neldermead
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   tbl.mass          9.00932   -0.0457346    0.0271102
   tbl.rate         0.294398    -0.024382    0.0318352
   tbl.angle        0.974253   -0.0106672        -----
tbl.angle lower bound:	-0.0106672
tbl.mass lower bound:	-0.0457346
tbl.angle upper bound:	-----
tbl.mass upper bound:	0.0271102
tbl.rate lower bound:	-0.024382
tbl.rate +: WARNING: The confidence level lies within (3.262462e-01, 3.262202e-01)
tbl.rate upper bound:	0.0318352

The screen output contains the confidence levels as well as the print out of the runs. These values can be accessed later as well using get_conf_results()

In [10]:
print get_conf_results()
datasets    = (1,)
methodname  = confidence
iterfitname = none
fitname     = neldermead
statname    = chi2gehrels
sigma       = 1
percent     = 68.2689492137
parnames    = ('tbl.mass', 'tbl.rate', 'tbl.angle')
parvals     = (9.0093186041805478, 0.29439800604025079, 0.9742525366691086)
parmins     = (-0.045734558207220388, -0.024381957835282408, -0.010667221038111796)
parmaxes    = (0.027110207851803736, 0.031835194148028267, None)
nfits       = 103
In [12]:
print get_conf_results().parmins
(-0.045734558207220388, -0.024381957835282408, -0.010667221038111796)
In [13]:
print get_model()
plot_fit()
template.tbl
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   tbl.mass     thawed      9.00932            6           10           
   tbl.rate     thawed     0.294398         0.01            1           
   tbl.angle    thawed     0.974253          0.1            1           
In [14]:
ignore()
notice(14.5,15.4)
plot_fit()