Modeling and Fitting SED Data in Iris

After one or more SED data segments and/or photometric points have been read into Iris, and data display preferences set, the data may be fit with a customized model expression using the Iris Fitting Tool. Clicking the Fitting Tool icon on the Iris desktop opens a new window in which individual model components may be added in order to define a custom model expression to fit to the data, and where initial model parameter values and the spectral fitting range are set. In Iris, model amplitudes are in units of photon flux density by default, and model spectral coordinates are in Angstroms.

When the fitting window opens, the default model that is provided, powerlaw.c1, is automatically plotted as a red line, and is calculated using the default power-law parameters. (At this point, the model has not been fit to the data, therefore it appears offset from the data curve.)

The Component section of the fitting window is the place to list the model components which will be used to construct the full model expression for fitting; these components can then be arbitrarily combined to form the full model expression, in the Model Expression field. The Component field is populated with a power-law model component by default, with model identifier "c1". The model component naming convention uses an ordered component list from top to bottom, starting with "c1". Removing various model components from the list will cause remaining components to potentially be renamed according to their position in the list, e.g., "c2", "c3", and so on.

Initial model parameter values and ranges may be set by selecting a model in the Component list, and then "Edit". This will open a new window containing a list of the parameters for the selected model component, with associated fields for editing.

Entering a value into a model parameter field, and then pressing the "Return" key, will update that model parameter in the Component section of the fitting window. The "Fit" checkbox next to each parameter in the editing window is for specifying whether or not to freeze that parameter during the fit; if checked, the parameter will be allowed to vary. Selecting the button labeled with a blue arrow beneath the 'Fit' checkbox opens a separate window in which parameter units, minimum and maximum values, and links to other parameters may be specified.

The "Fix all" and "Default Fit" options in the main fitting window may be used (while the model editor window remains open) to freeze all model parameters at current values, or thaw all model parameters (even those initially fixed), respectively.

The initial model parameter values are where Sherpa starts the fit in parameter space; the parameter ranges are the bounds of parameter space. It is not possible to fit using parameter values outside the specified parameter ranges.

The next step in the preparation for fitting involves choosing a fit statistic and optimization method appropriate for your analysis. Selecting "Fit" in the main fitting window (where model components are listed) launches a new model fitting window (labeled "Fitted") containing the statistic and method options. These are listed below with brief descriptions (refer to the Statistics and Optimization pages of the Sherpa website for a detailed explanation of each of the Sherpa fit statistics and methods available in Iris).

Fit Optimization Methods

• Levenberg-Marquardt - a method to find the best-fit model parameter values, by finding the local minimum in parameter space; the functions to be minimized are nonlinear least-squares functions of the model parameters.
• Nelder-Mead - a method to find the best-fit model parameter values, by finding the local minimum in parameter space; a direct search method is used to continually move "downhill" in parameter space, until settling in a local minimum.
• Monte Carlo - a method using the differential evolution algorithm, to find a global minimum in parameter space; the search is seeded with randomly selected starting points, and can continue the search for the true minimum where the other methods could get stuck in local minima that do not represent the actual best fit. (Most useful for complex models and fits that could explore complicated regions of parameter space.)

Fit Statistics

• least-squares - Sum of the squares of the differences between data and model values.
• Cash - A maximum likelihood function based on Poisson statistics. More suitable for counts data than for fluxes.
• C-stat - A maximum likelihood function similar to Cash, but with a chi-squared-like probability distribution. More suitable for counts data than for fluxes.
• chi2datavar - Chi-squared statistic with variance computed from the data. If measured errors are provided, the variance is taken from these errors; else, the variance is computed from the y-values of the data points. This is the preferred chi-squared option to be used in Iris.
• chi2gehrels - Chi-squared statistic, where the variance is computed with a function from Gehrels et al. Suitable for low counts data (e.g., X-ray data) to correct for bias in using chi-squared.
• chi2modvar - Chi-squared statistic, where the variance is computed from the model values instead of data.
• chi2constvar - Chi-squared statistic, where the variance is set to be a constant value. That constant is the average of the y-values of the data points.
• chi2xspecvar - Chi-squared statistic, where the variance is computed as the X-ray spectral fitting program XSPEC would compute the variance (i.e., where the variance would be less than one, reset it to one). More suitable for low counts data (e.g., X-ray data).

The Iris default fitting method and statistic are "neldermead" and "leastsq", respectively, which represent good choices for a robust, quick, initial fit of a relatively simple model to a data set covering potentially many orders of magnitude in flux and/or wavelength. The fit can also be done with a chi-squared statistic, with various methods for estimating variance, or with either of two maximum likelihood statistics that are useful when the data have low numbers of counts. The fitting method can be changed to "levmar" or "moncar", but note that switching the method is less important than switching the statistic, e.g., from least-squares to chi-squared.

Nelder-Mead is the optimal fitting method to start with because it does not depend on taking derivatives of the model function. As for the fit statistic, least-squares is preferred because it does not use measured errors and thus essentially weights all data points equally (it seeks to minimize the sum of the squares of the differences between data and model values, without taking errors into account). Chi-squared fitting, in contrast, can be decidedly biased towards the data points with the lowest fluxes and smallest measured error bars (when there are data points with measured errors that are orders of magnitude smaller than the error bars on most of the data points, then these points become the biggest contributors to the chi-squared value).

In subsequent fits - e.g., after selecting a different data range to fit, adding in new model components as needed, and so on - the statistic can be switched to a chi-squared option to use measured errors (preferably chi2datavar, which is chi-squared with data variance). Since measured errors are provided, this means the variance is taken directly from the errors provided with the data, when the data file is read in.

Note: When one or multiple SED segments is fit in Iris, any data points with associated zero-value errors are ignored in the fit. This design choice is intended as a safeguard against yielding potentially misleading fit results in the analysis of your SED data. Iris interprets a zero-value error not to mean that the uncertainty on the associated data value is actually zero, but that a measure of the uncertainty is not available for that particular photometric point.

Before initiating a fit of the defined model to a SED in Iris, the specific subset of the SED data to be included in the fit (if not the entire SED) should be specified using the "Define Range" option in the main fitting window. Selecting this option will prompt you to click on the data display twice, once to define the lower endpoint of the data range to be fit, and once for the upper endpoint.

Having built a model expression and set initial parameter values, chosen an appropriate fit statistic and optimization method for your analysis, and defining the range of data to be fit, you are ready to initiate the fitting process by selecting "Start" in the Fitted window. When the fit is complete, a red model curve appears overplotted on the SED data in the Iris main display (beyond the specified data range to be fit; this extraneous portion of the fitted model may be ignored).

When the fit has finished, the model parameter values will appear updated in the Component field of the fitting window, and fit statistics will be displayed in the fitting tab of the Fitted window, including the number of data points used in the fit, the degrees of freedom, and the probability that the fit is consistent with the data (q-value), where applicable (when fitting with the least-squares or Cash statistic, the "reduced statistic" value and q-value are not available, as they cannot be computed).

When fitting with one of the chi-squared statisics, or the chi-squared-like cstat statistic, the Confidence tab returns the requested sigma/percent confidence intervals on the best-fit model parameters, calculated using the statistic and method chosen for the fit. For example, entering "1.6" in the "sigma" field will return th 90% confidence limits on the model parameter values..

If the statistic is left at the default "leastsq", Iris will issue an error stating that least-squares cannot be used with the confidence limit function (this is because the only way calculate to confidence limits is to know the probability distribution for the fit statistic, which is known for the chi-squared and chi-squared-lke cstat statistic, but not least squares).

Also of note is that blank or "NaN" values are returned in the confidence results when a parameter bound is found to lie outside the hard limit boundary for a model parameter (which ought not to be changed by the user, so there is not an option to do so in Iris). This could result from an issue with the signal-to-noise of the data, the applicability of the model to the data, systematic errors in the data, among others things. A parameter hard limit represents either a hard physical limit (e.g., temperature is not allowed to go below zero), a mathematical limit (e.g., prevent a number from going to zero or below, when the log of that number will be taken), or the limit of what a float or double can hold (the fit should not be driven above or below the maximum or minimum values a variable can hold).

You can iterate through the fitting process in Iris as many times as necessary - adding or deleting models from the model expression; changing parameter values or ranges; including or excluding new points from the SED - until you find a satisfactory model that best describes your SED data.

Once a spectral model is satisfactorily fitted, the custom fitted model parameters can be saved to a local file before exiting the fitting session, by making an appropriate selection in the File menu of the fitting window.

Selecting the File -> Write to text file option saves the spectral model to file in a human-readable text format. The "Write active components" option in this submenu writes only the model components which were used in the fit; example output is shown below:

$ more 3c273_atten_bpl_active_comps.txt

File: 3C 273

Thu Feb 26 15:46:31 EDT 2015  Iris 2

TARGNAME: 3C 273

Model Expression:
        c1*c2

Fit parameters:
        Final fit statistic:      0.010714114984176791
        Reduced statistic:        nan
        Probability [Q-value]:    nan
        Degrees of freedom:       413.0
        Data points:              419
        Last function evaluation: 394

Component 1:    atten.c1 
       hcol      = 1.0E20     (NaN NaN)                                   
       heiRatio  = 0.04177875 (NaN NaN)                                   
       heiiRatio = 0.74654233 (NaN NaN)                                   

Component 2:    brokenpowerlaw.c2 
   F   refer  = 300000.0                            angstroms             
       ampl   = 0.0062282784 (0.1065081 NaN)                              
       index1 = -0.039390083 (NaN       -4.82435)                         
       index2 = -0.6063001   (2.4033892 NaN)                              

Selecting File -> Write to file, instead, saves the model to a file in a format (CDB) that can be read back by the tool in a future data analysis session. (CDB is an XML format for reading model files in and out of Iris.)

This customized Iris fitting session may be restored by selecting the File->Read from file menu option within the Iris fitting window.

After fitting a SED, you can estimate fluxes through user-defined passbands and photometric filters by integrating under the model components. This is done from the "Calculate Flux" tab in the Science tools, which is opened by clicking the "Shift, Interpolate, Integrate" icon on the desktop. The user may integrate under the full model or an arbitrary combination of the model components. The Fitting Tool must remain open for model integration.

After opening the "Calculate Flux" tab, make sure Model Integration is selected (it should "Integrate Model (YES)"). By default, the full model will be used to evaluate the fluxes. For our fit of 3c273, the model expression is "c1*c2". A quick-look view of the model and the parameter values can be displayed by clicking the "Show Model" button.

You can change the number of spectral bins used to calculate the integral by adjusting the value in the "Model evaluation bin #" box.

By checking-off "Full model," the user can arbitrarily combine the model components, then integrate the expression in the "Model Expression" field. This means that users can integrate under individual model components. Multiplication and addition of other model components and numeric values are acceptable. For example, the expressions "c2 * 3.678", "c2 + c1", and "c1 + c2 * 1.2" are allowed.

Say we want to estimate the total infrared flux and the flux through the Herschel PACS and SPIRE bands . Going back to the "Calculate Flux" in the Science tool, we add the Hershcel photometry filters (select "Photometry Filter" -> "Choose" and search for "Herschel") and add a user-defined pass band from 0.8 to 1000 microns (8000 - 1E7 Angstroms). To integrate under the fitted model, we turn "Model Integration" on (YES), and click "Calculate."

The results can be exported as a new SED with the "Create SED" button, or may be saved to a text file with "Save." The text file can be loaded back into Iris as an ASCII Table.

NOTE: If the user exits the main Fitting Tool window, the user will no longer be able to integrate under the fitted model, as the models will be lost from closing the fitting session. Iris will ask the user if they really want to leave the fitting session if they try to close the Fitting Tool.