Fitting Models to Spectra

Whew! This was a little tougher than expected, and I struggled with it quite a bit. The next phase of my research will be focusing on determining the relative strengths of emission and absorption lines in order to characterize the classification parameters of observed galaxies. For this, I will eventually be using a “Baldwin-Phillips-Terlevich” Diagram (BPT).  By measuring the ratios of the observed selected emission/absorption lines, we can empirically calculate indicators of the primary excitation mechanism in the observed galaxies. A great, and easy-to-read rundown of the application of BPT diagrams can be found here.  But, first, you have to determine the location of the targeted emission/absorption lines and measure their intensity.  This is where fitting models to spectra becomes handy.

Fig 1: Line Fitting For Hydrogen Alpha Line. The observed flux, on the right, is nearly covered by the model generated, visually indicating an accurate fit. The continuum spectrum uses a median-smoothing algorithm to produce a generic continuum. An unfortunate drawback to this method results in a less accurate fit when approaching the bounds of the spectrum. By confining the measured Hydrogen Alpha region to within 2-sigma, a more accurate calculation of the equivalent width can be performed.
Fig 2: Line Fitting for [OIII] Line. Again, this turned out to be a really good fit.
Fig 3: Line Fitting for Hydrogen Beta. Here, you can see that there is some deviation in the fit of the model, but is pretty consistent within 2-Sigma.

I thought about writing a tutorial, but a great one already exists. I found this “Quick Fit” modeling tutorial to be extremely useful when working with model parameters, especially when performing error analysis. Once you feel comfortable with the basics, I recommend moving on to the “User-defined Model” tutorial. In this tutorial, you’ll learn the tools necessary to fit Gaussian lines to emission/absorption lines using more complex models. Using this method, you can fit a Gaussian line to the emission/absorption lines, AND fit a one-degree polynomial to estimate the local continuum. Super useful! One thing to mention that wasn’t included in the tutorial: the model fitters give you the ability to have certain information returned to you, including a covariance matrix for the error in your fit parameters. Let’s look at In[22] from the User-defined Model tutorial.

Screen Grab of In[22] from the tutorial

Using this as a starting point, we can add a few lines of code to extract more information. Add the following lines to your code.

print(fitter.fit_info)

The output from ‘fitter.fit_info’.

There are an awful lot of numbers there. Let’s just focus on the parameter covariance matrix. We can call the keyword ‘param_cov’ from the fitter to retrieve the matrix. Next, we can “diagonalize” the matrix, and take the square-root of the values to determine the uncertainty in the associated parameters and then print the parameters associated with those errors..

print(np.sqrt(np.diag(fitter.fit_info['param_cov']))
print(compound_fit_bounded.param_names)
Errors in the fit are show in Out[9] and their associated parameters are shown in Out[10].

Armed with this information, you can ensure that your model is statistically sound.